8 #define BLOCK_SIZE_NT 4 /* one block holds 4 nucleotides. */
9 #define BITS_IN_NT 2 /* two bits holds 1 nucleotide. */
10 #define BITS_IN_BYTE 8 /* number of bits in one byte. */
11 #define BLOCK_SPACE_MAX 64 /* maximum space between two blocks. */
12 #define BLOCK_MASK ( ( BLOCK_SPACE_MAX << 1 ) - 1 ) /* mask for printing block space. */
13 #define COUNT_ARRAY_NMEMB ( 1 << 30 ) /* number of objects in the unsigned int count array. */
14 #define CUTOFF 10000 /* minimum number of motifs in output. */
16 #define add_A( c ) /* add 00 to the rightmost two bits of bin (i.e. do nothing). */
17 #define add_T( c ) ( c |= 3 ) /* add 11 on the rightmost two bits of c. */
18 #define add_C( c ) ( c |= 1 ) /* add 01 on the rightmost two bits of c. */
19 #define add_G( c ) ( c |= 2 ) /* add 10 on the rightmost two bits of c. */
21 /* Structure that will hold one tetra nucleotide block. */
24 uchar bin; /* Tetra nucleotide binary encoded. */
25 bool hasN; /* Flag indicating any N's in the block. */
28 typedef struct _bitblock bitblock;
30 /* Byte array for fast convertion of binary blocks back to DNA. */
31 char *bin2dna[256] = {
32 "AAAA", "AAAC", "AAAG", "AAAT", "AACA", "AACC", "AACG", "AACT",
33 "AAGA", "AAGC", "AAGG", "AAGT", "AATA", "AATC", "AATG", "AATT",
34 "ACAA", "ACAC", "ACAG", "ACAT", "ACCA", "ACCC", "ACCG", "ACCT",
35 "ACGA", "ACGC", "ACGG", "ACGT", "ACTA", "ACTC", "ACTG", "ACTT",
36 "AGAA", "AGAC", "AGAG", "AGAT", "AGCA", "AGCC", "AGCG", "AGCT",
37 "AGGA", "AGGC", "AGGG", "AGGT", "AGTA", "AGTC", "AGTG", "AGTT",
38 "ATAA", "ATAC", "ATAG", "ATAT", "ATCA", "ATCC", "ATCG", "ATCT",
39 "ATGA", "ATGC", "ATGG", "ATGT", "ATTA", "ATTC", "ATTG", "ATTT",
40 "CAAA", "CAAC", "CAAG", "CAAT", "CACA", "CACC", "CACG", "CACT",
41 "CAGA", "CAGC", "CAGG", "CAGT", "CATA", "CATC", "CATG", "CATT",
42 "CCAA", "CCAC", "CCAG", "CCAT", "CCCA", "CCCC", "CCCG", "CCCT",
43 "CCGA", "CCGC", "CCGG", "CCGT", "CCTA", "CCTC", "CCTG", "CCTT",
44 "CGAA", "CGAC", "CGAG", "CGAT", "CGCA", "CGCC", "CGCG", "CGCT",
45 "CGGA", "CGGC", "CGGG", "CGGT", "CGTA", "CGTC", "CGTG", "CGTT",
46 "CTAA", "CTAC", "CTAG", "CTAT", "CTCA", "CTCC", "CTCG", "CTCT",
47 "CTGA", "CTGC", "CTGG", "CTGT", "CTTA", "CTTC", "CTTG", "CTTT",
48 "GAAA", "GAAC", "GAAG", "GAAT", "GACA", "GACC", "GACG", "GACT",
49 "GAGA", "GAGC", "GAGG", "GAGT", "GATA", "GATC", "GATG", "GATT",
50 "GCAA", "GCAC", "GCAG", "GCAT", "GCCA", "GCCC", "GCCG", "GCCT",
51 "GCGA", "GCGC", "GCGG", "GCGT", "GCTA", "GCTC", "GCTG", "GCTT",
52 "GGAA", "GGAC", "GGAG", "GGAT", "GGCA", "GGCC", "GGCG", "GGCT",
53 "GGGA", "GGGC", "GGGG", "GGGT", "GGTA", "GGTC", "GGTG", "GGTT",
54 "GTAA", "GTAC", "GTAG", "GTAT", "GTCA", "GTCC", "GTCG", "GTCT",
55 "GTGA", "GTGC", "GTGG", "GTGT", "GTTA", "GTTC", "GTTG", "GTTT",
56 "TAAA", "TAAC", "TAAG", "TAAT", "TACA", "TACC", "TACG", "TACT",
57 "TAGA", "TAGC", "TAGG", "TAGT", "TATA", "TATC", "TATG", "TATT",
58 "TCAA", "TCAC", "TCAG", "TCAT", "TCCA", "TCCC", "TCCG", "TCCT",
59 "TCGA", "TCGC", "TCGG", "TCGT", "TCTA", "TCTC", "TCTG", "TCTT",
60 "TGAA", "TGAC", "TGAG", "TGAT", "TGCA", "TGCC", "TGCG", "TGCT",
61 "TGGA", "TGGC", "TGGG", "TGGT", "TGTA", "TGTC", "TGTG", "TGTT",
62 "TTAA", "TTAC", "TTAG", "TTAT", "TTCA", "TTCC", "TTCG", "TTCT",
63 "TTGA", "TTGC", "TTGG", "TTGT", "TTTA", "TTTC", "TTTG", "TTTT"
66 /* Function declarations. */
67 void run_scan( int argc, char *argv[] );
69 void scan_file( char *file, seq_entry *entry, uint *count_array );
70 uint *count_array_new( size_t nmemb );
71 void scan_seq( char *seq, size_t seq_len, uint *count_array );
72 void scan_list( list_sl *list, uint *count_array );
73 bitblock *bitblock_new();
74 uint blocks2motif( uchar bin1, uchar bin2, ushort dist );
75 int cmp_uint_desc( const void *a, const void *b );
76 void count_array_print( uint *count_array, size_t nmemb, size_t cutoff );
77 void motif_print( uint motif, uint count );
78 void bitblock_list_print( list_sl *list );
79 void bitblock_print( bitblock *out );
81 /* Unit test declarations. */
82 static void run_tests();
83 static void test_count_array_new();
84 static void test_bitblock_new();
85 static void test_bitblock_print();
86 static void test_bitblock_list_print();
87 static void test_scan_seq();
88 static void test_blocks2motif();
91 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MAIN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
94 int main( int argc, char *argv[] )
102 run_scan( argc, argv );
108 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FUNCTIONS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
113 /* Martin A. Hansen, September 2008 */
115 /* Print usage and exit if no files in argument. */
118 "Usage: bipartite_scam <FASTA file(s)> > result.csv\n"
121 exit( EXIT_SUCCESS );
125 void run_scan( int argc, char *argv[] )
127 /* Martin A. Hansen, September 2008 */
129 /* For each file in argv scan the file for */
130 /* bipartite motifs and output the motifs */
131 /* and their count. */
135 seq_entry *entry = NULL;
136 uint *count_array = NULL;
137 // size_t new_nmemb = 0;
139 count_array = count_array_new( COUNT_ARRAY_NMEMB );
141 entry = seq_new( MAX_SEQ_NAME, MAX_SEQ );
143 for ( i = 1; i < argc; i++ )
147 fprintf( stderr, "Scanning file: %s\n", file );
148 scan_file( file, entry, count_array );
149 fprintf( stderr, "done.\n" );
152 fprintf( stderr, "Printing motifs: ... " );
153 count_array_print( count_array, COUNT_ARRAY_NMEMB, CUTOFF );
154 fprintf( stderr, "done.\n" );
156 seq_destroy( entry );
158 mem_free( &count_array );
162 uint *count_array_new( size_t nmemb )
164 /* Martin A. Hansen, September 2008 */
166 /* Initialize a new zeroed uint array of nmemb objects. */
172 array = mem_get_zero( nmemb * sizeof( uint ) );
178 void scan_file( char *file, seq_entry *entry, uint *count_array )
180 /* Martin A. Hansen, September 2008 */
182 /* Scan all FASTA entries of a file in both */
183 /* sense and anti-sense directions. */
185 FILE *fp = read_open( file );
187 while ( fasta_get_entry( fp, &entry ) == TRUE )
189 fprintf( stderr, " Scanning: %s (%zu nt) ... ", entry->seq_name, entry->seq_len );
191 scan_seq( entry->seq, entry->seq_len, count_array );
193 fprintf( stderr, "done.\n" );
200 void scan_seq( char *seq, size_t seq_len, uint *count_array )
202 /* Martin A. Hansen, September 2008 */
204 /* Run a sliding window over a given sequence. The window */
205 /* consists of a list where new blocks of 4 nucleotides */
206 /* are pushed onto one end while at the same time old */
207 /* blocks are popped from the other end. The number of */
208 /* in the list is determined by the maximum seperator. */
209 /* Everytime we have a full window, the window is scanned */
212 bitblock *block = NULL;
217 bool first_node = TRUE;
218 node_sl *new_node = NULL;
219 node_sl *old_node = NULL;
220 list_sl *list = list_sl_new();
222 for ( i = 0; seq[ i ]; i++ )
228 case 'A': case 'a': add_A( bin ); break;
229 case 'T': case 't': add_T( bin ); break;
230 case 'C': case 'c': add_C( bin ); break;
231 case 'G': case 'g': add_G( bin ); break;
232 default: n_count = BLOCK_SIZE_NT; break;
235 if ( i > BLOCK_SIZE_NT - 2 )
239 block = bitblock_new();
248 new_node = node_sl_new();
249 new_node->val = block;
252 list_sl_add_beg( &list, &new_node );
254 list_sl_add_after( &old_node, &new_node );
261 if ( b_count > BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
263 // bitblock_list_print( list ); /* DEBUG */
265 scan_list( list, count_array );
267 list_sl_remove_beg( &list );
272 /* if the list is shorter than BLOCK_SPACE_MAX + BLOCK_SIZE_NT */
273 if ( b_count <= BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
275 // bitblock_list_print( list ); /* DEBUG */
277 scan_list( list, count_array );
280 list_sl_destroy( &list );
284 void scan_list( list_sl *list, uint *count_array )
286 /* Martin A. Hansen, September 2008 */
288 /* Scan a list of blocks for biparite motifs by creating */
289 /* a binary motif consisting of two blocks of 4 nucleotides */
290 /* along with the distance separating them. Motifs containing */
291 /* N's are skipped. */
293 node_sl *first_node = NULL;
294 node_sl *next_node = NULL;
295 bitblock *block1 = NULL;
296 bitblock *block2 = NULL;
301 // bitblock_list_print( list );
303 first_node = list->first;
305 block1 = ( bitblock * ) first_node->val;
307 if ( ! block1->hasN )
309 next_node = first_node->next;
311 for ( i = 0; i < BLOCK_SIZE_NT - 1; i++ ) {
312 next_node = next_node->next;
315 for ( next_node = next_node; next_node != NULL; next_node = next_node->next )
317 block2 = ( bitblock * ) next_node->val;
319 // printf( "block1: %s block2: %s dist: %d\n", bin2dna[ block1->bin ], bin2dna[ block2->bin ], dist ); /* DEBUG */
321 if ( ! block2->hasN )
323 motif_bin = blocks2motif( block1->bin, block2->bin, dist );
325 // motif_print( motif_bin, 0 ); /* DEBUG */
326 // bitblock_list_print( list ); /* DEBUG */
328 count_array[ motif_bin ]++;
337 bitblock *bitblock_new()
339 /* Martin A. Hansen, September 2008 */
341 /* Initializes a new empty bitblock. */
343 bitblock *new_block = NULL;
345 new_block = mem_get( sizeof( bitblock ) );
348 new_block->hasN = FALSE;
354 uint blocks2motif( uchar bin1, uchar bin2, ushort dist )
356 /* Martin A. Hansen, September 2008 */
358 /* Given two binary encoded tetra nuceotide blocks, */
359 /* and the distance separating this, create a binary */
360 /* bipartite motif. */
366 motif <<= sizeof( uchar ) * BITS_IN_BYTE;
370 motif <<= sizeof( uchar ) * BITS_IN_BYTE;
374 // motif_print( motif, 0 ); /* DEBUG */
380 void count_array_print( uint *count_array, size_t nmemb, size_t cutoff )
382 /* Martin A. Hansen, Seqptember 2008. */
384 /* Print all bipartite motifs in count_array as */
385 /* tabular output. */
391 for ( i = 0; i < nmemb; i++ )
394 count = count_array[ i ];
396 if ( count >= cutoff ) {
397 motif_print( motif, count );
403 void motif_print( uint motif, uint count )
405 /* Martin A. Hansen, September 2008 */
407 /* Converts a binary encoded bipartite motif */
408 /* into DNA and output the motif, distance and */
409 /* count seperated by tabs: */
410 /* BLOCK1 \t BLOCK2 \t DIST \t COUNT */
416 // printf( "%d\t", motif ); /* DEBUG */
418 dist = ( ushort ) motif & BLOCK_MASK;
420 motif >>= sizeof( uchar ) * BITS_IN_BYTE;
422 bin2 = ( uchar ) motif;
424 motif >>= sizeof( uchar ) * BITS_IN_BYTE;
426 bin1 = ( uchar ) motif;
428 printf( "%s\t%s\t%d\t%d\n", bin2dna[ bin1 ], bin2dna[ bin2 ], dist, count );
432 void bitblock_list_print( list_sl *list )
434 /* Martin A. Hansen, September 2008 */
436 /* Debug function to print all blocks in a list. */
438 node_sl *node = NULL;
440 printf( "\nbitblock_list_print:\n" );
442 for ( node = list->first; node != NULL; node = node->next ) {
443 bitblock_print( ( bitblock * ) node->val );
448 void bitblock_print( bitblock *out )
450 /* Martin A. Hansen, September 2008 */
452 /* Debug function to print a given block. */
454 printf( "bin: %d dna: %s hasN: %d\n", out->bin, bin2dna[ ( int ) out->bin ], out->hasN );
458 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> UNIT TESTS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
463 fprintf( stderr, "Running tests\n" );
465 test_count_array_new();
467 test_bitblock_print();
468 test_bitblock_list_print();
472 fprintf( stderr, "All tests OK\n" );
476 void test_count_array_new()
478 fprintf( stderr, " Running test_count_array_new ... " );
484 array = count_array_new( nmemb );
486 for ( i = 0; i < nmemb; i++ ) {
487 assert( array[ i ] == 0 );
492 fprintf( stderr, "done.\n" );
496 void test_bitblock_new()
498 fprintf( stderr, " Running test_bitblock_new ... " );
500 bitblock *new_block = bitblock_new();
502 assert( new_block->bin == 0 );
503 assert( new_block->hasN == FALSE );
505 fprintf( stderr, "done.\n" );
509 void test_bitblock_print()
511 fprintf( stderr, " Running test_bitblock_print ... " );
513 bitblock *new_block = bitblock_new();
516 new_block->hasN = TRUE;
518 // bitblock_print( new_block );
520 fprintf( stderr, "done.\n");
524 void test_bitblock_list_print()
526 fprintf( stderr, " Running test_bitblock_list_print ... " );
528 list_sl *list = list_sl_new();
529 node_sl *node1 = node_sl_new();
530 node_sl *node2 = node_sl_new();
531 node_sl *node3 = node_sl_new();
532 bitblock *block1 = bitblock_new();
533 bitblock *block2 = bitblock_new();
534 bitblock *block3 = bitblock_new();
547 list_sl_add_beg( &list, &node1 );
548 list_sl_add_beg( &list, &node2 );
549 list_sl_add_beg( &list, &node3 );
551 // bitblock_list_print( list );
553 fprintf( stderr, "done.\n" );
559 fprintf( stderr, " Running test_scan_seq ... " );
561 //char *seq = "AAAANTCGGCTNGGGG";
562 //char *seq = "AAAATCGGCTGGGG";
563 char *seq = "AAAAAAAAAAAAAAG";
564 size_t seq_len = strlen( seq );
565 size_t nmemb = 1 << 5;
567 uint *count_array = count_array_new( nmemb );
569 scan_seq( seq, seq_len, count_array );
571 // count_array_print( count_array, nmemb, 1 ); /* DEBUG */
573 fprintf( stderr, "done.\n" );
577 static void test_blocks2motif()
579 fprintf( stderr, " Running test_blocks2motif ... " );
586 motif = blocks2motif( bin1, bin2, dist );
588 // printf( "motif: %d\n", motif );
590 fprintf( stderr, "done.\n");
594 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */