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. */
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 );
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"
26 "The output is a fixedStep file ala the phastCons files from the UCSC\n"
29 "Usage: repeat-O-matic [options] <FASTA file>\n"
32 " [-n <int> | --nmer <int>] # nmer size between 1 and 15 (Default 15).\n"
33 " [-1 | --log10] # output log10 (Default no).\n"
36 " repeat-O-matic -n 14 -l hg18.fna > hg18.fixedStep\n"
38 "Copyright (C) 2008, Martin A. Hansen\n"
46 int main( int argc, char *argv[] )
50 bool log10_flag = FALSE;
55 static struct option longopts[] = {
56 { "nmer", required_argument, NULL, 'n' },
57 { "log10", no_argument, NULL, 'l' },
61 while ( ( opt = getopt_long( argc, argv, "n:l", longopts, NULL ) ) != -1 )
64 case 'n': nmer = strtol( optarg, NULL, 0 ); break;
65 case 'l': log10_flag = TRUE; break;
73 if ( nmer < 1 || nmer > 15 )
75 fprintf( stderr, "ERROR: nmer must be between 1 and 15 inclusive - not %d\n", nmer );
83 path = argv[ argc - 1 ];
85 mask = mask_create( nmer );
87 oligo_count( path, &array, nmer, mask );
89 oligo_count_output( path, array, nmer, mask, log10_flag );
95 uint mask_create( int oligo_size )
97 /* Martin A. Hansen, June 2008 */
99 /* Create a bit mask for binary encoded oligos less than 32 bits. */
104 for ( i = 0; i < oligo_size; i++ )
108 mask |= 3; /* add 11 to mask */
115 void oligo_count( char *path, uint **array_ppt, uint nmer, uint mask )
117 /* Martin A. Hansen, June 2008 */
119 /* Count the occurence of all oligos of a fixed size in a FASTA file. */
121 uint *array = *array_ppt;
126 seq_entry *entry = NULL;
129 array_size = ( 1 << ( nmer * 2 ) );
130 array = mem_get_zero( sizeof( uint ) * array_size );
132 fp = read_open( path );
134 seq_new( &entry, MAX_SEQ_NAME, MAX_SEQ );
136 while ( ( fasta_get_entry( fp, &entry ) ) != 0 )
138 fprintf( stderr, "Counting oligos in: %s ... ", entry->seq_name );
140 /* ---- Sense strand ---- */
145 for ( i = 0; entry->seq[ i ]; i++ )
149 switch( entry->seq[ i ] )
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;
160 array[ ( bin & mask ) ]++;
163 printf( "mask : %s\n", bits2string( mask ) );
164 printf( "bin : %s\n", bits2string( bin ) );
165 printf( "bin & mask: %s\n", bits2string( bin & mask ) );
170 /* ---- Anti-sense strand ---- */
172 revcomp_dna( entry->seq );
177 for ( i = 0; entry->seq[ i ]; i++ )
181 switch( entry->seq[ i ] )
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;
192 array[ ( bin & mask ) ]++;
195 printf( "mask : %s\n", bits2string( mask ) );
196 printf( "bin : %s\n", bits2string( bin ) );
197 printf( "bin & mask: %s\n", bits2string( bin & mask ) );
202 fprintf( stderr, "done.\n" );
207 free( entry->seq_name );
215 void oligo_count_output( char *path, uint *array, uint nmer, uint mask, bool log10_flag )
217 /* Martin A. Hansen, June 2008 */
219 /* Output oligo count for each sequence position. */
233 seq_new( &entry, MAX_SEQ_NAME, MAX_SEQ );
235 fp = read_open( path );
237 while ( ( fasta_get_entry( fp, &entry ) ) != 0 )
239 fprintf( stderr, "Writing results for: %s ... ", entry->seq_name );
244 block_size = sizeof( uint ) * ( entry->seq_len + nmer );
245 block = mem_get_zero( block_size );
247 for ( i = 0; entry->seq[ i ]; i++ )
251 switch( entry->seq[ i ] )
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;
262 count = array[ ( bin & mask ) ];
266 chr_pos = i - nmer + 1;
268 if ( block_pos == 0 )
270 memset( block, '\0', block_size );
274 block[ block_pos ] = count;
280 if ( chr_pos > block_beg + block_pos )
282 fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos, log10_flag );
288 block[ block_pos ] = count;
299 fixedstep_put_entry( entry->seq_name, block_beg, 1, block, block_pos, log10_flag );
301 mem_free( ( void * ) &block );
304 fprintf( stderr, "done.\n" );
307 free( entry->seq_name );
315 void fixedstep_put_entry( char *chr, int beg, int step_size, uint *block_array, int block_size, bool log10_flag )
317 /* Martin A. Hansen, June 2008 */
319 /* Outputs a block of fixedStep values. */
325 if ( block_size > 0 )
327 beg += 1; /* fixedStep format is 1 based. */
329 printf( "fixedStep chrom=%s start=%d step=%d\n", chr, beg, step_size );
331 for ( i = 0; i < block_size; i++ ) {
332 printf( "%lf\n", log10( block_array[ i ] ) );
338 if ( block_size > 0 )
340 beg += 1; /* fixedStep format is 1 based. */
342 printf( "fixedStep chrom=%s start=%d step=%d\n", chr, beg, step_size );
344 for ( i = 0; i < block_size; i++ ) {
345 printf( "%i\n", block_array[ i ] );