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