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

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