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

   1 #include <getopt.h>
   2 #include "common.h"
   3 #include "mem.h"
   4 #include "filesys.h"
   5 #include "seq.h"
   6 #include "fasta.h"
   7
   8 #define add_A( c )            /* add 00 to the rightmost two bits of bin (i.e. do nothing). */
   9 #define add_T( c ) ( c |= 3 ) /* add 11 on the rightmost two bits of c. */
  10 #define add_C( c ) ( c |= 1 ) /* add 01 on the rightmost two bits of c. */
  11 #define add_G( c ) ( c |= 2 ) /* add 10 on the rightmost two bits of c. */
  12
  13 static uint mask_create( int oligo_size );
  14 static void oligo_count( char *path, uint **array_ppt, uint nmer, uint mask );
  15 static void oligo_count_output( char *path, uint *array, uint nmer, uint mask, bool log10_flag );
  16 static void fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size, bool log10_flag );
  17
  18
  19 static void usage()
  20 {
  21     fprintf( stderr,
  22         "\n"
  23         "repeat-O-matic determines the repetiveness of a genome by determining\n"
  24         "the number of identical n-mers for each position in the genome.\n"
  25         "\n"
  26         "The output is a fixedStep file ala the phastCons files from the UCSC\n"
  27         "Genome browser.\n"
  28         "\n"
  29         "Usage: repeat-O-matic [options] <FASTA file>\n"
  30         "\n"
  31         "Options:\n"
  32         "   [-n <int> | --nmer <int>]   # nmer size between 1 and 15 (Default 15).\n"
  33         "   [-1       | --log10]        # output log10 (Default no).\n"
  34         "\n"
  35         "Examples:\n"
  36         "   repeat-O-matic -n 14 -l hg18.fna > hg18.fixedStep\n"
  37         "\n"
  38         "Copyright (C) 2008, Martin A. Hansen\n"
  39         "\n"
  40         );
  41
  42     exit( EXIT_FAILURE );
  43 }
  44
  45
  46 int main( int argc, char *argv[] )
  47 {
  48     int   opt        = 0;
  49     uint  nmer       = 15;
  50     bool  log10_flag = FALSE;
  51     char *path       = NULL;
  52     uint  mask       = 0;
  53     uint *array      = NULL;
  54
  55     static struct option longopts[] = {
  56         { "nmer",  required_argument, NULL, 'n' },
  57         { "log10", no_argument,       NULL, 'l' },
  58         { NULL,    0,                 NULL,  0 }
  59     };
  60
  61     while ( ( opt = getopt_long( argc, argv, "n:l", longopts, NULL ) ) != -1 )
  62     {
  63         switch ( opt ) {
  64             case 'n': nmer  = strtol( optarg, NULL, 0 ); break;
  65             case 'l': log10_flag = TRUE;                 break;
  66             default:                                     break;
  67         }
  68     }
  69
  70     argc -= optind;
  71     argv += optind;
  72
  73     if ( nmer < 1 || nmer > 15 )
  74     {
  75         fprintf( stderr, "ERROR: nmer must be between 1 and 15 inclusive - not %d\n", nmer );
  76         abort();
  77     }
  78
  79     if ( argc < 1 ) {
  80         usage();
  81     }
  82
  83     path = argv[ argc - 1 ];
  84
  85     mask = mask_create( nmer );
  86
  87     oligo_count( path, &array, nmer, mask );
  88
  89     oligo_count_output( path, array, nmer, mask, log10_flag );
  90
  91     return EXIT_SUCCESS;
  92 }
  93
  94
  95 uint mask_create( int oligo_size )
  96 {
  97     /* Martin A. Hansen, June 2008 */
  98
  99     /* Create a bit mask for binary encoded oligos less than 32 bits. */
 100
 101     uint i    = 0;
 102     uint mask = 0;
 103
 104     for ( i = 0; i < oligo_size; i++ )
 105     {
 106         mask <<= 2;
 107
 108         mask |= 3; /* add 11 to mask */
 109     }
 110
 111     return mask;
 112 }
 113
 114
 115 void oligo_count( char *path, uint **array_ppt, uint nmer, uint mask )
 116 {
 117     /* Martin A. Hansen, June 2008 */
 118
 119     /* Count the occurence of all oligos of a fixed size in a FASTA file. */
 120
 121     uint       *array      = *array_ppt;
 122     uint        array_size = 0;
 123     uint        i          = 0;
 124     uint        j          = 0;
 125     uint        bin        = 0;
 126     seq_entry  *entry      = NULL;
 127     FILE       *fp         = NULL;
 128
 129     array_size = ( 1 << ( nmer * 2 ) );
 130     array = mem_get_zero( sizeof( uint ) * array_size );
 131
 132     fp = read_open( path );
 133
 134     seq_new( &entry, MAX_SEQ_NAME, MAX_SEQ );
 135
 136     while ( ( fasta_get_entry( fp, &entry ) ) != 0 )
 137     {
 138         fprintf( stderr, "Counting oligos in: %s ... ", entry->seq_name );
 139
 140         /* ---- Sense strand ---- */
 141
 142         bin = 0;
 143         j   = 0;
 144
 145         for ( i = 0; entry->seq[ i ]; i++ )
 146         {
 147             bin <<= 2;
 148
 149             switch( entry->seq[ i ] )
 150             {
 151                 case 'A': case 'a': add_A( bin ); j++; break;
 152                 case 'T': case 't': add_T( bin ); j++; break;
 153                 case 'C': case 'c': add_C( bin ); j++; break;
 154                 case 'G': case 'g': add_G( bin ); j++; break;
 155                 default: bin = 0; j = 0; break;
 156             }
 157
 158             if ( j >= nmer )
 159             {
 160                 array[ ( bin & mask ) ]++;
 161 /*
 162                 printf( "\n" );
 163                 printf( "mask          : %s\n", bits2string( mask ) );
 164                 printf( "bin           : %s\n", bits2string( bin ) );
 165                 printf( "bin     & mask: %s\n", bits2string( bin & mask ) );
 166 */
 167             }
 168         }
 169
 170         /* ---- Anti-sense strand ---- */
 171
 172         revcomp_dna( entry->seq );
 173
 174         bin = 0;
 175         j   = 0;
 176
 177         for ( i = 0; entry->seq[ i ]; i++ )
 178         {
 179             bin <<= 2;
 180
 181             switch( entry->seq[ i ] )
 182             {
 183                 case 'A': case 'a': add_A( bin ); j++; break;
 184                 case 'T': case 't': add_T( bin ); j++; break;
 185                 case 'C': case 'c': add_C( bin ); j++; break;
 186                 case 'G': case 'g': add_G( bin ); j++; break;
 187                 default: bin = 0; j = 0; break;
 188             }
 189
 190             if ( j >= nmer )
 191             {
 192                 array[ ( bin & mask ) ]++;
 193 /*
 194                 printf( "\n" );
 195                 printf( "mask          : %s\n", bits2string( mask ) );
 196                 printf( "bin           : %s\n", bits2string( bin ) );
 197                 printf( "bin     & mask: %s\n", bits2string( bin & mask ) );
 198 */
 199             }
 200         }
 201
 202         fprintf( stderr, "done.\n" );
 203     }
 204
 205     close_stream( fp );
 206
 207     free( entry->seq_name );
 208     free( entry->seq );
 209     entry = NULL;
 210
 211     *array_ppt = array;
 212 }
 213
 214
 215 void oligo_count_output( char *path, uint *array, uint nmer, uint mask, bool log10_flag )
 216 {
 217     /* Martin A. Hansen, June 2008 */
 218
 219     /* Output oligo count for each sequence position. */
 220
 221     uint       i;
 222     uint       j;
 223     uint       bin;
 224     int        count;
 225     uint      *block;
 226     uint       block_pos;
 227     uint       block_beg;
 228     uint       block_size;
 229     uint       chr_pos;
 230     seq_entry *entry;
 231     FILE      *fp;
 232
 233     seq_new( &entry, MAX_SEQ_NAME, MAX_SEQ );
 234
 235     fp = read_open( path );
 236
 237     while ( ( fasta_get_entry( fp, &entry ) ) != 0 )
 238     {
 239         fprintf( stderr, "Writing results for: %s ... ", entry->seq_name );
 240
 241         bin        = 0;
 242         j          = 0;
 243         block_pos  = 0;
 244         block_size = sizeof( uint ) * ( entry->seq_len + nmer );
 245         block      = mem_get_zero( block_size );
 246
 247         for ( i = 0; entry->seq[ i ]; i++ )
 248         {
 249             bin <<= 2;
 250
 251             switch( entry->seq[ i ] )
 252             {
 253                 case 'A': case 'a': add_A( bin ); j++; break;
 254                 case 'T': case 't': add_T( bin ); j++; break;
 255                 case 'C': case 'c': add_C( bin ); j++; break;
 256                 case 'G': case 'g': add_G( bin ); j++; break;
 257                 default: bin = 0; j = 0; break;
 258             }
 259
 260             if ( j >= nmer )
 261             {
 262                 count = array[ ( bin & mask ) ];
 263
 264                 if ( count > 1 )
 265                 {
 266                     chr_pos = i - nmer + 1;
 267
 268                     if ( block_pos == 0 )
 269                     {
 270                         block_beg = chr_pos;
 271
 272                         block[ block_pos ] = count;
 273
 274                         block_pos++;
 275                     }
 276                     else
 277                     {
 278                         if ( chr_pos > block_beg + block_pos )
 279                         {
 280                             fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos, log10_flag );
 281
 282                             block_pos = 0;
 283
 284                             memset( block, '\0', block_pos );
 285                         }
 286                         else
 287                         {
 288                             block[ block_pos ] = count;
 289
 290                             block_pos++;
 291                         }
 292                     }
 293                 }
 294             }
 295         }
 296
 297         if ( block_pos > 0 )
 298         {
 299             fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos, log10_flag );
 300
 301             free( block );
 302             block = NULL;
 303         }
 304
 305         fprintf( stderr, "done.\n" );
 306     }
 307
 308     free( entry->seq_name );
 309     free( entry->seq );
 310     entry = NULL;
 311
 312     close_stream( fp );
 313 }
 314
 315
 316 void fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size, bool log10_flag )
 317 {
 318     /* Martin A. Hansen, June 2008 */
 319
 320     /* Outputs a block of fixedStep values. */
 321
 322     int i;
 323
 324     if ( log10_flag )
 325     {
 326         if ( block_size > 0 )
 327         {
 328             beg += 1; /* fixedStep format is 1 based. */
 329
 330             printf( "fixedStep chrom=%s start=%d step=%d\n", chr, beg, step_size );
 331
 332             for ( i = 0; i < block_size; i++ ) {
 333                 printf( "%lf\n", log10( block_array[ i ] ) );
 334             }
 335         }
 336     }
 337     else
 338     {
 339         if ( block_size > 0 )
 340         {
 341             beg += 1; /* fixedStep format is 1 based. */
 342
 343             printf( "fixedStep chrom=%s start=%d step=%d\n", chr, beg, step_size );
 344
 345             for ( i = 0; i < block_size; i++ ) {
 346                 printf( "%i\n", block_array[ i ] );
 347             }
 348         }
 349     }
 350 }
 351
 352