code_c/Maasha/src/bipartite_scan.c

   1 #include "common.h"
   2 #include "mem.h"
   3 #include "filesys.h"
   4 #include "seq.h"
   5 #include "fasta.h"
   6 #include "list.h"
   7
   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 COUNT_ARRAY_SIZE ( sizeof( uint ) * ( 1 << 30 ) )  /* size of the unsigned int count array. */
  13
  14 #define add_A( c )              /* add 00 to the rightmost two bits of bin (i.e. do nothing). */
  15 #define add_T( c ) ( c |= 3 )   /* add 11 on the rightmost two bits of c. */
  16 #define add_C( c ) ( c |= 1 )   /* add 01 on the rightmost two bits of c. */
  17 #define add_G( c ) ( c |= 2 )   /* add 10 on the rightmost two bits of c. */
  18
  19 /* Structure that will hold one tetra nucleotide block. */
  20 struct _bitblock
  21 {
  22     uchar bin;   /* Tetra nucleotide encoded binary. */
  23     bool hasN;   /* Flag indicating any N's in the block. */
  24 };
  25
  26 typedef struct _bitblock bitblock;
  27
  28 /* Byte array for fast convertion of binary blocks back to DNA. */
  29 char *bin2dna[256] = {
  30     "AAAA", "AAAC", "AAAG", "AAAT", "AACA", "AACC", "AACG", "AACT",
  31     "AAGA", "AAGC", "AAGG", "AAGT", "AATA", "AATC", "AATG", "AATT",
  32     "ACAA", "ACAC", "ACAG", "ACAT", "ACCA", "ACCC", "ACCG", "ACCT",
  33     "ACGA", "ACGC", "ACGG", "ACGT", "ACTA", "ACTC", "ACTG", "ACTT",
  34     "AGAA", "AGAC", "AGAG", "AGAT", "AGCA", "AGCC", "AGCG", "AGCT",
  35     "AGGA", "AGGC", "AGGG", "AGGT", "AGTA", "AGTC", "AGTG", "AGTT",
  36     "ATAA", "ATAC", "ATAG", "ATAT", "ATCA", "ATCC", "ATCG", "ATCT",
  37     "ATGA", "ATGC", "ATGG", "ATGT", "ATTA", "ATTC", "ATTG", "ATTT",
  38     "CAAA", "CAAC", "CAAG", "CAAT", "CACA", "CACC", "CACG", "CACT",
  39     "CAGA", "CAGC", "CAGG", "CAGT", "CATA", "CATC", "CATG", "CATT",
  40     "CCAA", "CCAC", "CCAG", "CCAT", "CCCA", "CCCC", "CCCG", "CCCT",
  41     "CCGA", "CCGC", "CCGG", "CCGT", "CCTA", "CCTC", "CCTG", "CCTT",
  42     "CGAA", "CGAC", "CGAG", "CGAT", "CGCA", "CGCC", "CGCG", "CGCT",
  43     "CGGA", "CGGC", "CGGG", "CGGT", "CGTA", "CGTC", "CGTG", "CGTT",
  44     "CTAA", "CTAC", "CTAG", "CTAT", "CTCA", "CTCC", "CTCG", "CTCT",
  45     "CTGA", "CTGC", "CTGG", "CTGT", "CTTA", "CTTC", "CTTG", "CTTT",
  46     "GAAA", "GAAC", "GAAG", "GAAT", "GACA", "GACC", "GACG", "GACT",
  47     "GAGA", "GAGC", "GAGG", "GAGT", "GATA", "GATC", "GATG", "GATT",
  48     "GCAA", "GCAC", "GCAG", "GCAT", "GCCA", "GCCC", "GCCG", "GCCT",
  49     "GCGA", "GCGC", "GCGG", "GCGT", "GCTA", "GCTC", "GCTG", "GCTT",
  50     "GGAA", "GGAC", "GGAG", "GGAT", "GGCA", "GGCC", "GGCG", "GGCT",
  51     "GGGA", "GGGC", "GGGG", "GGGT", "GGTA", "GGTC", "GGTG", "GGTT",
  52     "GTAA", "GTAC", "GTAG", "GTAT", "GTCA", "GTCC", "GTCG", "GTCT",
  53     "GTGA", "GTGC", "GTGG", "GTGT", "GTTA", "GTTC", "GTTG", "GTTT",
  54     "TAAA", "TAAC", "TAAG", "TAAT", "TACA", "TACC", "TACG", "TACT",
  55     "TAGA", "TAGC", "TAGG", "TAGT", "TATA", "TATC", "TATG", "TATT",
  56     "TCAA", "TCAC", "TCAG", "TCAT", "TCCA", "TCCC", "TCCG", "TCCT",
  57     "TCGA", "TCGC", "TCGG", "TCGT", "TCTA", "TCTC", "TCTG", "TCTT",
  58     "TGAA", "TGAC", "TGAG", "TGAT", "TGCA", "TGCC", "TGCG", "TGCT",
  59     "TGGA", "TGGC", "TGGG", "TGGT", "TGTA", "TGTC", "TGTG", "TGTT",
  60     "TTAA", "TTAC", "TTAG", "TTAT", "TTCA", "TTCC", "TTCG", "TTCT",
  61     "TTGA", "TTGC", "TTGG", "TTGT", "TTTA", "TTTC", "TTTG", "TTTT"
  62 };
  63
  64 /* Function declarations. */
  65 void      run_scan( int argc, char *argv[] );
  66 void      print_usage();
  67 void      scan_file( char *file, seq_entry *entry, uint *count_array );
  68 uint     *count_array_new( size_t size );
  69 void      scan_seq( char *seq, size_t seq_len, uint *count_array );
  70 void      scan_list( list_sl *list, uint *count_array );
  71 bitblock *bitblock_new();
  72 uint      blocks2motif( uchar bin1, uchar bin2, ushort dist );
  73 void      count_array_print( uint *count_array );
  74 void      motif_print( uint motif, uint count );
  75 void      bitblock_list_print( list_sl *list );
  76 void      bitblock_print( bitblock *out );
  77
  78 /* Unit test declarations. */
  79 static void run_tests();
  80 static void test_count_array_new();
  81 static void test_bitblock_new();
  82 static void test_bitblock_print();
  83 static void test_bitblock_list_print();
  84 static void test_scan_seq();
  85 static void test_blocks2motif();
  86
  87
  88 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MAIN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
  89
  90
  91 int main( int argc, char *argv[] )
  92 {
  93     if ( argc == 1 ) {
  94         print_usage();
  95     }
  96
  97     run_tests();
  98     run_scan( argc, argv );
  99
 100     return EXIT_SUCCESS;
 101 }
 102
 103
 104 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FUNCTIONS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 105
 106
 107 void print_usage()
 108 {
 109     /* Martin A. Hansen, September 2008 */
 110
 111     /* Print usage and exit if no files in argument. */
 112
 113     fprintf( stderr,
 114         "Usage: bipartite_scam <FASTA file(s)> > result.csv\n"
 115     );
 116
 117     exit( EXIT_SUCCESS );
 118 }
 119
 120
 121 void run_scan( int argc, char *argv[] )
 122 {
 123     /* Martin A. Hansen, September 2008 */
 124
 125     /* For each file in argv scan the file for */
 126     /* bipartite motifs and output the motifs */
 127     /* and their count. */
 128
 129     char      *file        = NULL;
 130     int        i           = 0;
 131     seq_entry *entry       = NULL;
 132     uint      *count_array = NULL;
 133
 134     count_array = count_array_new( COUNT_ARRAY_SIZE );
 135
 136     entry = seq_new( MAX_SEQ_NAME, MAX_SEQ );
 137
 138     for ( i = 1; i < argc; i++ )
 139     {
 140         file = argv[ i ];
 141
 142         fprintf( stderr, "Scanning file: %s\n", file );
 143
 144         scan_file( file, entry, count_array );
 145
 146         fprintf( stderr, "done.\n" );
 147     }
 148
 149     fprintf( stderr, "Printing motifs: ... " );
 150
 151     count_array_print( count_array );
 152
 153     fprintf( stderr, "done.\n" );
 154
 155     seq_destroy( entry );
 156
 157     mem_free( &count_array );
 158 }
 159
 160
 161 uint *count_array_new( size_t size )
 162 {
 163     /* Martin A. Hansen, September 2008 */
 164
 165     /* Initialize a new zeroed uint array of a given byte size. */
 166
 167     uint *array = NULL;
 168
 169     array = mem_get_zero( size );
 170
 171     return array;
 172 }
 173
 174
 175 void scan_file( char *file, seq_entry *entry, uint *count_array )
 176 {
 177     /* Martin A. Hansen, September 2008 */
 178
 179     /* Scan all FASTA entries of a file in both */
 180     /* sense and anti-sense directions. */
 181
 182     FILE *fp = read_open( file );
 183
 184     while ( fasta_get_entry( fp, &entry ) == TRUE )
 185     {
 186         fprintf( stderr, "   Scanning: %s (%zu nt) ... ", entry->seq_name, entry->seq_len );
 187
 188         scan_seq( entry->seq, entry->seq_len, count_array );
 189
 190         fprintf( stderr, "done.\n" );
 191     }
 192
 193     close_stream( fp );
 194 }
 195
 196
 197 void scan_seq( char *seq, size_t seq_len, uint *count_array )
 198 {
 199     /* Martin A. Hansen, September 2008 */
 200
 201     /* Run a sliding window over a given sequence. The window */
 202     /* consists of a list where new blocks of 4 nucleotides */
 203     /* are pushed onto one end while at the same time old */
 204     /* blocks are popped from the other end. The number of */
 205     /* in the list is determined by the maximum seperator. */
 206     /* Everytime we have a full window, the window is scanned */
 207     /* for motifs. */
 208
 209     bitblock *block      = NULL;
 210     size_t    b_count    = 0;
 211     ushort    n_count    = 0;
 212     size_t    i          = 0;
 213     uchar     bin        = 0;
 214     bool      first_node = TRUE;
 215     node_sl *new_node    = NULL;
 216     node_sl *old_node    = NULL;
 217     list_sl *list        = list_sl_new();
 218
 219     for ( i = 0; seq[ i ]; i++ )
 220     {
 221         bin <<= BITS_IN_NT;
 222
 223         switch( seq[ i ] )
 224         {
 225             case 'A': case 'a': add_A( bin ); break;
 226             case 'T': case 't': add_T( bin ); break;
 227             case 'C': case 'c': add_C( bin ); break;
 228             case 'G': case 'g': add_G( bin ); break;
 229             default: n_count = BLOCK_SIZE_NT; break;
 230         }
 231
 232         if ( i > BLOCK_SIZE_NT - 2 )
 233         {
 234             b_count++;
 235
 236             block      = bitblock_new();
 237             block->bin = bin;
 238
 239             if ( n_count > 0 )
 240             {
 241                  block->hasN = TRUE;
 242                  n_count--;
 243             }
 244
 245             new_node      = node_sl_new();
 246             new_node->val = block;
 247
 248             if ( first_node ) {
 249                 list_sl_add_beg( &list, &new_node );
 250             } else {
 251                 list_sl_add_after( &old_node, &new_node );
 252             }
 253
 254             old_node = new_node;
 255
 256             first_node = FALSE;
 257
 258             if ( b_count > BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
 259             {
 260                 // bitblock_list_print( list );
 261
 262                 scan_list( list, count_array );
 263
 264                 list_sl_remove_beg( &list );
 265             }
 266         }
 267     }
 268
 269     list_sl_destroy( &list );
 270 }
 271
 272
 273 void scan_list( list_sl *list, uint *count_array )
 274 {
 275     /* Martin A. Hansen, September 2008 */
 276
 277     /* Scan a list of blocks for biparite motifs by creating */
 278     /* a binary motif consisting of two blocks of 4 nucleotides */
 279     /* along with the distance separating them. Motifs containing */
 280     /* N's are skipped. */
 281
 282     node_sl  *first_node = NULL;
 283     node_sl  *next_node  = NULL;
 284     bitblock *block1     = NULL;
 285     bitblock *block2     = NULL;
 286     int       i          = 0;
 287     ushort    dist       = 0;
 288     uint      motif_bin  = 0;
 289
 290 //    bitblock_list_print( list );
 291
 292     first_node = list->first;
 293
 294     block1 = ( bitblock * ) first_node->val;
 295
 296     if ( ! block1->hasN )
 297     {
 298         next_node = first_node->next;
 299
 300         for ( i = 0; i < BLOCK_SIZE_NT - 1; i++ ) {
 301             next_node = next_node->next;
 302         }
 303
 304         for ( next_node = next_node; next_node != NULL; next_node = next_node->next )
 305         {
 306             block2 = ( bitblock * ) next_node->val;
 307
 308             // printf( "block1: %s  block2: %s   dist: %d\n", bin2dna[ block1->bin ], bin2dna[ block2->bin ], dist ); /* DEBUG */
 309
 310             if ( ! block2->hasN )
 311             {
 312                 motif_bin = blocks2motif( block1->bin, block2->bin, dist );
 313
 314 //                motif_print( motif_bin, 0 ); /* DEBUG */
 315
 316                 count_array[ motif_bin ]++;
 317             }
 318
 319             dist++;
 320         }
 321     }
 322 }
 323
 324
 325 bitblock *bitblock_new()
 326 {
 327     /* Martin A. Hansen, September 2008 */
 328
 329     /* Initializes a new empty bitblock. */
 330
 331     bitblock *new_block = NULL;
 332
 333     new_block = mem_get( sizeof( bitblock ) );
 334
 335     new_block->bin  = 0;
 336     new_block->hasN = FALSE;
 337
 338     return new_block;
 339 }
 340
 341
 342 uint blocks2motif( uchar bin1, uchar bin2, ushort dist )
 343 {
 344     /* Martin A. Hansen, September 2008 */
 345
 346     /* Given two binary encoded tetra nuceotide blocks, */
 347     /* and the distance separating this, create a binary */
 348     /* bipartite motif. */
 349
 350     uint motif = 0;
 351
 352     motif |= bin1;
 353
 354     motif <<= sizeof( uchar ) * BITS_IN_BYTE;
 355
 356     motif |= bin2;
 357
 358     motif <<= sizeof( uchar ) * BITS_IN_BYTE;
 359
 360     motif |= dist;
 361
 362     return motif;
 363 }
 364
 365
 366 void count_array_print( uint *count_array )
 367 {
 368     /* Martin A. Hansen, Seqptember 2008. */
 369
 370     /* Print all bipartite motifs in count_array as */
 371     /* tabular output. */
 372
 373     uint i     = 0;
 374     uint motif = 0;
 375     uint count = 0;
 376
 377     for ( i = 0; i < COUNT_ARRAY_SIZE / 4; i++ )
 378     {
 379         motif = i;
 380         count = count_array[ i ];
 381
 382         if ( count > 0 ) {
 383             motif_print( motif, count );
 384         }
 385     }
 386 }
 387
 388
 389 void motif_print( uint motif, uint count )
 390 {
 391     /* Martin A. Hansen, September 2008 */
 392
 393     /* Converts a binary encoded bipartite motif */
 394     /* into DNA and output the motif, distance and */
 395     /* count seperated by tabs: */
 396     /* BLOCK1 \t BLOCK2 \t DIST \t COUNT */
 397
 398     uchar  bin1 = 0;
 399     uchar  bin2 = 0;
 400     ushort dist = 0;
 401
 402     dist = ( ushort ) motif & 255;
 403
 404     motif >>= sizeof( uchar ) * BITS_IN_BYTE;
 405
 406     bin2 = ( uchar ) motif;
 407
 408     motif >>= sizeof( uchar ) * BITS_IN_BYTE;
 409
 410     bin1 = ( uchar ) motif;
 411
 412     printf( "%s\t%s\t%d\t%d\n", bin2dna[ bin1 ], bin2dna[ bin2 ], dist, count );
 413 }
 414
 415
 416 void bitblock_list_print( list_sl *list )
 417 {
 418     /* Martin A. Hansen, September 2008 */
 419
 420     /* Debug function to print all blocks in a list. */
 421
 422     node_sl *node = NULL;
 423
 424     printf( "\nbitblock_list_print:\n" );
 425
 426     for ( node = list->first; node != NULL; node = node->next ) {
 427         bitblock_print( ( bitblock * ) node->val );
 428     }
 429 }
 430
 431
 432 void bitblock_print( bitblock *out )
 433 {
 434     /* Martin A. Hansen, September 2008 */
 435
 436     /* Debug function to print a given block. */
 437
 438     printf( "bin: %d  dna: %s   hasN: %d\n", out->bin, bin2dna[ ( int ) out->bin ], out->hasN );
 439 }
 440
 441
 442 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> UNIT TESTS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 443
 444
 445 void run_tests()
 446 {
 447     fprintf( stderr, "Running tests\n" );
 448
 449     test_count_array_new();
 450     test_bitblock_new();
 451     test_bitblock_print();
 452     test_bitblock_list_print();
 453     test_scan_seq();
 454     test_blocks2motif();
 455
 456     fprintf( stderr, "All tests OK\n" );
 457 }
 458
 459
 460 void test_count_array_new()
 461 {
 462     fprintf( stderr, "   Running test_count_array_new ... " );
 463
 464     size_t  i     = 0;
 465     size_t  size  = 128;
 466     uint   *array = NULL;
 467
 468     array = count_array_new( 4 * size );
 469
 470     for ( i = 0; i < size; i++ ) {
 471         assert( array[ i ] == 0 );
 472     }
 473
 474     mem_free( &array );
 475
 476     fprintf( stderr, "done.\n" );
 477 }
 478
 479
 480 void test_bitblock_new()
 481 {
 482     fprintf( stderr, "   Running test_bitblock_new ... " );
 483
 484     bitblock *new_block = bitblock_new();
 485
 486     assert( new_block->bin  == 0 );
 487     assert( new_block->hasN == FALSE );
 488
 489     fprintf( stderr, "done.\n" );
 490 }
 491
 492
 493 void test_bitblock_print()
 494 {
 495     fprintf( stderr, "   Running test_bitblock_print ... " );
 496
 497     bitblock *new_block = bitblock_new();
 498
 499     new_block->bin  = 7;
 500     new_block->hasN = TRUE;
 501
 502 //    bitblock_print( new_block );
 503
 504     fprintf( stderr, "done.\n");
 505 }
 506
 507
 508 void test_bitblock_list_print()
 509 {
 510     fprintf( stderr, "   Running test_bitblock_list_print ... " );
 511
 512     list_sl  *list   = list_sl_new();
 513     node_sl  *node1  = node_sl_new();
 514     node_sl  *node2  = node_sl_new();
 515     node_sl  *node3  = node_sl_new();
 516     bitblock *block1 = bitblock_new();
 517     bitblock *block2 = bitblock_new();
 518     bitblock *block3 = bitblock_new();
 519
 520     block1->bin  = 0;
 521     block1->hasN = TRUE;
 522     block2->bin  = 1;
 523     block2->hasN = TRUE;
 524     block3->bin  = 2;
 525     block3->hasN = TRUE;
 526
 527     node1->val  = block1;
 528     node2->val  = block2;
 529     node3->val  = block3;
 530
 531     list_sl_add_beg( &list, &node1 );
 532     list_sl_add_beg( &list, &node2 );
 533     list_sl_add_beg( &list, &node3 );
 534
 535     // bitblock_list_print( list );
 536
 537     fprintf( stderr, "done.\n" );
 538 }
 539
 540
 541 void test_scan_seq()
 542 {
 543     fprintf( stderr, "   Running test_scan_seq ... " );
 544
 545     //char   *seq       = "AAAANTCGGCTNGGGG";
 546     //char   *seq         = "AAAATCGGCTGGGG";
 547     char   *seq         = "AAAAAAAAAAAAAAA";
 548     size_t  seq_len     = strlen( seq );
 549     uint   *count_array = mem_get_zero( sizeof( uint ) * ( 1 << 5 ) );
 550
 551     scan_seq( seq, seq_len, count_array );
 552
 553     fprintf( stderr, "done.\n" );
 554 }
 555
 556
 557 static void test_blocks2motif()
 558 {
 559     fprintf( stderr, "   Running test_blocks2motif ... " );
 560
 561     uchar  bin1  = 4;
 562     uchar  bin2  = 3;
 563     ushort dist  = 256;
 564     uint   motif = 0;
 565
 566     motif = blocks2motif( bin1, bin2, dist );
 567
 568 //    printf( "motif: %d\n", motif );
 569
 570     fprintf( stderr, "done.\n");
 571 }
 572
 573
 574 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 575
 576