code_c/Maasha/src/repeat-O-matic.c

   1 /*
   2     Copyright (C) 2008, Martin A. Hansen
   3
   4     This program determines the repetiveness of a genome by determining
   5     the number of identical 15-mers for each position in the genome.
   6
   7     The output is a fixedStep file ala the phastCons files from the UCSC
   8     Genome browser.
   9
  10     It is very fast and efficient using less than 8 Gb of memory to
  11     complete the human genome in roughly 30 minutes.
  12 */
  13
  14 #include "common.h"
  15 #include "filesys.h"
  16 #include "fasta.h"
  17
  18 // #define OLIGO_SIZE 15
  19 #define OLIGO_SIZE 5
  20 #define SIZE ( 1 << ( OLIGO_SIZE * 2 ) )
  21
  22 #define UINT_BITS 32
  23 #define T    3                            /* 11 on the rightmost two bits of bin. */
  24 #define C    1                            /* 01 on the rightmost two bits of bin. */
  25 #define G    2                            /* 10 on the rightmost two bits of bin. */
  26
  27 uint *oligo_count( char *path );
  28 uint  mask_create( int oligo_size );
  29 void  oligo_count_output( char *path, uint *array );
  30 void  fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size );
  31
  32 int main( int argc, char *argv[] )
  33 {
  34     char *path  = NULL;
  35     uint *array = NULL;
  36
  37     path  = argv[ 1 ];
  38
  39     array = oligo_count( path );
  40
  41     oligo_count_output( path, array );
  42
  43     return 0;
  44 }
  45
  46
  47 uint *oligo_count( char *path )
  48 {
  49     /* Martin A. Hansen, June 2008 */
  50
  51     /* Count the occurence of all oligos of a fixed size in a FASTA file. */
  52
  53     uint             *array   = NULL;
  54     uint              i       = 0;
  55     uint              mask    = 0;
  56     uint              bin     = 0;
  57     uint              bin_rc1 = 0;
  58     uint              bin_rc2 = 0;
  59     uint              j       = 0;
  60     uint              A_rc    = ( 3 << ( UINT_BITS - 2 ) );   /* 11 on the leftmost two bits an uint. */
  61     uint              G_rc    = ( 2 << ( UINT_BITS - 2 ) );   /* 10 on the leftmost two bits an uint. */
  62     uint              C_rc    = ( 1 << ( UINT_BITS - 2 ) );   /* 01 on the leftmost two bits an uint. */
  63     struct seq_entry *entry   = NULL;
  64     FILE             *fp      = NULL;
  65
  66     array = mem_get_zero( sizeof( uint ) * SIZE );
  67
  68     mask  = mask_create( OLIGO_SIZE );
  69
  70     MEM_GET( entry );
  71
  72     fp = read_open( path );
  73
  74     while ( ( fasta_get_entry( fp, entry ) ) )
  75     {
  76         fprintf( stderr, "Counting oligos in: %s ... ", entry->seq_name );
  77
  78         bin     = 0;
  79         bin_rc1 = 0;
  80         j       = 0;
  81
  82         for ( i = 0; entry->seq[ i ]; i++ )
  83         {
  84             bin     <<= 2;
  85             bin_rc1 >>= 2;
  86
  87             switch( entry->seq[ i ] )
  88             {
  89                 case 'A': case 'a':           bin_rc1 |= A_rc; j++; break;
  90                 case 'T': case 't': bin |= T;                  j++; break;
  91                 case 'C': case 'c': bin |= C; bin_rc1 |= G_rc; j++; break;
  92                 case 'G': case 'g': bin |= G; bin_rc1 |= C_rc; j++; break;
  93                 default: bin = 0; bin_rc1 = 0; j = 0; break;
  94             }
  95
  96             if ( j >= OLIGO_SIZE )
  97             {
  98                 array[ ( bin & mask ) ]++;
  99
 100                 bin_rc2 = bin_rc1;
 101
 102                 bin_rc2 >>= ( UINT_BITS - ( OLIGO_SIZE * 2 ) );
 103
 104                 array[ ( bin_rc2 ) ]++;
 105 /*
 106                 printf( "\n" );
 107                 printf( "mask          : %s\n", bits2string( mask ) );
 108                 printf( "bin           : %s\n", bits2string( bin ) );
 109                 printf( "bin     & mask: %s\n", bits2string( bin & mask ) );
 110                 printf( "bin_rc1       : %s\n", bits2string( bin_rc1 ) );
 111                 printf( "bin_rc2       : %s\n", bits2string( bin_rc2 ) );
 112 */
 113             }
 114         }
 115
 116         fprintf( stderr, "done.\n" );
 117     }
 118
 119     close_stream( fp );
 120
 121     fasta_free_entry( entry );
 122
 123     return array;
 124 }
 125
 126
 127 uint mask_create( int oligo_size )
 128 {
 129     /* Martin A. Hansen, June 2008 */
 130
 131     /* Create a bit mask for binary encode oligos less than sizeof( uint ). */
 132
 133     uint i;
 134     uint mask;
 135
 136     mask = 0;
 137
 138     for ( i = 0; i < oligo_size; i++ )
 139     {
 140         mask <<= 2;
 141
 142         mask |= 3;
 143     }
 144
 145     return mask;
 146 }
 147
 148
 149 void oligo_count_output( char *path, uint *array )
 150 {
 151     /* Martin A. Hansen, June 2008 */
 152
 153     /* Output oligo count for each sequence position. */
 154
 155     struct seq_entry *entry;
 156     FILE             *fp;
 157     uint              mask;
 158     uint              i;
 159     uint              j;
 160     uint              bin;
 161     int               count;
 162     uint             *block;
 163     uint              block_pos;
 164     uint              block_beg;
 165     uint              chr_pos;
 166
 167     mask = mask_create( OLIGO_SIZE );
 168
 169     MEM_GET( entry );
 170
 171     fp = read_open( path );
 172
 173     while ( ( fasta_get_entry( fp, entry ) ) )
 174     {
 175         fprintf( stderr, "Writing results for: %s ... ", entry->seq_name );
 176
 177         bin       = 0;
 178         j         = 0;
 179         block_pos = 0;
 180         block     = mem_get_zero( sizeof( uint ) * ( entry->seq_len + OLIGO_SIZE ) );
 181
 182         for ( i = 0; entry->seq[ i ]; i++ )
 183         {
 184             bin <<= 2;
 185
 186             switch( entry->seq[ i ] )
 187             {
 188                 case 'A': case 'a':           j++; break;
 189                 case 'T': case 't': bin |= T; j++; break;
 190                 case 'C': case 'c': bin |= C; j++; break;
 191                 case 'G': case 'g': bin |= G; j++; break;
 192                 default: bin = 0; j = 0; break;
 193             }
 194
 195             if ( j >= OLIGO_SIZE )
 196             {
 197                 count   = array[ ( bin & mask ) ];
 198
 199                 if ( count > 1 )
 200                 {
 201                     chr_pos = i - OLIGO_SIZE + 1;
 202
 203                     if ( block_pos == 0 )
 204                     {
 205                         ZERO( block );
 206
 207                         block_beg = chr_pos;
 208
 209                         block[ block_pos ] = count;
 210
 211                         block_pos++;
 212                     }
 213                     else
 214                     {
 215                         if ( chr_pos > block_beg + block_pos )
 216                         {
 217                             fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos );
 218
 219                             block_pos = 0;
 220                         }
 221                         else
 222                         {
 223                             block[ block_pos ] = count;
 224
 225                             block_pos++;
 226                         }
 227                     }
 228                 }
 229             }
 230         }
 231
 232         if ( block_pos > 0 )
 233         {
 234             fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos );
 235
 236             mem_free( block );
 237         }
 238
 239         fprintf( stderr, "done.\n" );
 240     }
 241
 242     close_stream( fp );
 243
 244     fasta_free_entry( entry );
 245 }
 246
 247
 248 void fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size )
 249 {
 250     /* Martin A. Hansen, June 2008 */
 251
 252     /* Outputs a block of fixedStep values. */
 253
 254     int i;
 255
 256     if ( block_size > 0 )
 257     {
 258         beg += 1; /* fixedStep format is 1 based. */
 259
 260         printf( "fixedStep chrom=%s start=%d step=%d\n", chr, beg, step_size );
 261
 262         for ( i = 0; i < block_size; i++ ) {
 263             printf( "%u\n", block_array[ i ] );
 264         }
 265     }
 266 }