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     MEM_GET( 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              chr_pos;
 167
 168     mask = mask_create( OLIGO_SIZE );
 169
 170     MEM_GET( entry );
 171
 172     fp = read_open( path );
 173
 174     while ( ( fasta_get_entry( fp, entry ) ) )
 175     {
 176         fprintf( stderr, "Writing results for: %s ... ", entry->seq_name );
 177
 178         bin       = 0;
 179         j         = 0;
 180         block_pos = 0;
 181         block     = mem_get_zero( sizeof( uint ) * ( entry->seq_len + OLIGO_SIZE ) );
 182
 183         for ( i = 0; entry->seq[ i ]; i++ )
 184         {
 185             bin <<= 2;
 186
 187             switch( entry->seq[ i ] )
 188             {
 189                 case 'A': case 'a':           j++; break;
 190                 case 'T': case 't': bin |= T; j++; break;
 191                 case 'C': case 'c': bin |= C; j++; break;
 192                 case 'G': case 'g': bin |= G; j++; break;
 193                 default: bin = 0; j = 0; break;
 194             }
 195
 196             if ( j >= OLIGO_SIZE )
 197             {
 198                 count   = array[ ( bin & mask ) ];
 199
 200                 if ( count > 1 )
 201                 {
 202                     chr_pos = i - OLIGO_SIZE + 1;
 203
 204                     if ( block_pos == 0 )
 205                     {
 206                         MEM_ZERO( block );
 207
 208                         block_beg = chr_pos;
 209
 210                         block[ block_pos ] = count;
 211
 212                         block_pos++;
 213                     }
 214                     else
 215                     {
 216                         if ( chr_pos > block_beg + block_pos )
 217                         {
 218                             fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos );
 219
 220                             block_pos = 0;
 221                         }
 222                         else
 223                         {
 224                             block[ block_pos ] = count;
 225
 226                             block_pos++;
 227                         }
 228                     }
 229                 }
 230             }
 231         }
 232
 233         if ( block_pos > 0 )
 234         {
 235             fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos );
 236
 237             mem_free( block );
 238         }
 239
 240         fprintf( stderr, "done.\n" );
 241     }
 242
 243     close_stream( fp );
 244
 245     fasta_free_entry( entry );
 246 }
 247
 248
 249 void fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size )
 250 {
 251     /* Martin A. Hansen, June 2008 */
 252
 253     /* Outputs a block of fixedStep values. */
 254
 255     int i;
 256
 257     if ( block_size > 0 )
 258     {
 259         beg += 1; /* fixedStep format is 1 based. */
 260
 261         printf( "fixedStep chrom=%s start=%d step=%d\n", chr, beg, step_size );
 262
 263         for ( i = 0; i < block_size; i++ ) {
 264             printf( "%u\n", block_array[ i ] );
 265         }
 266     }
 267 }