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

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