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 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. */
  15
  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. */
  20
  21 /* Structure that will hold one tetra nucleotide block. */
  22 struct _bitblock
  23 {
  24     uchar bin;   /* Tetra nucleotide binary encoded. */
  25     bool hasN;   /* Flag indicating any N's in the block. */
  26 };
  27
  28 typedef struct _bitblock bitblock;
  29
  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"
  64 };
  65
  66 /* Function declarations. */
  67 void      run_scan( int argc, char *argv[] );
  68 void      print_usage();
  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 void      count_array_print( uint *count_array, size_t nmemb, size_t cutoff );
  76 void      motif_print( uint motif, uint count );
  77 void      bitblock_list_print( list_sl *list );
  78 void      bitblock_print( bitblock *out );
  79
  80 /* Unit test declarations. */
  81 static void run_tests();
  82 static void test_count_array_new();
  83 static void test_bitblock_new();
  84 static void test_bitblock_print();
  85 static void test_bitblock_list_print();
  86 static void test_scan_seq();
  87 static void test_blocks2motif();
  88
  89
  90 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MAIN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
  91
  92
  93 int main( int argc, char *argv[] )
  94 {
  95     run_tests();
  96
  97     if ( argc == 1 ) {
  98         print_usage();
  99     }
 100
 101     run_scan( argc, argv );
 102
 103     return EXIT_SUCCESS;
 104 }
 105
 106
 107 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FUNCTIONS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 108
 109
 110 void print_usage()
 111 {
 112     /* Martin A. Hansen, September 2008 */
 113
 114     /* Print usage and exit if no files in argument. */
 115
 116     fprintf( stderr,
 117         "Usage: bipartite_scam <FASTA file(s)> > result.csv\n"
 118     );
 119
 120     exit( EXIT_SUCCESS );
 121 }
 122
 123
 124 void run_scan( int argc, char *argv[] )
 125 {
 126     /* Martin A. Hansen, September 2008 */
 127
 128     /* For each file in argv scan the file for */
 129     /* bipartite motifs and output the motifs */
 130     /* and their count. */
 131
 132     char      *file        = NULL;
 133     int        i           = 0;
 134     seq_entry *entry       = NULL;
 135     uint      *count_array = NULL;
 136 //    size_t    new_nmemb    = 0;
 137
 138     count_array = count_array_new( COUNT_ARRAY_NMEMB );
 139
 140     entry = seq_new( MAX_SEQ_NAME, MAX_SEQ );
 141
 142     for ( i = 1; i < argc; i++ )
 143     {
 144         file = argv[ i ];
 145
 146         fprintf( stderr, "Scanning file: %s\n", file );
 147         scan_file( file, entry, count_array );
 148         fprintf( stderr, "done.\n" );
 149     }
 150
 151     fprintf( stderr, "Printing motifs: ... " );
 152     count_array_print( count_array, COUNT_ARRAY_NMEMB, CUTOFF );
 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 nmemb )
 162 {
 163     /* Martin A. Hansen, September 2008 */
 164
 165     /* Initialize a new zeroed uint array of nmemb objects. */
 166
 167     uint *array = NULL;
 168
 169     assert( nmemb > 0 );
 170
 171     array = mem_get_zero( nmemb * sizeof( uint ) );
 172
 173     return array;
 174 }
 175
 176
 177 void scan_file( char *file, seq_entry *entry, uint *count_array )
 178 {
 179     /* Martin A. Hansen, September 2008 */
 180
 181     /* Scan all FASTA entries of a file in both */
 182     /* sense and anti-sense directions. */
 183
 184     FILE *fp = read_open( file );
 185
 186     while ( fasta_get_entry( fp, &entry ) == TRUE )
 187     {
 188         fprintf( stderr, "   Scanning: %s (%zu nt) ... ", entry->seq_name, entry->seq_len );
 189
 190         scan_seq( entry->seq, entry->seq_len, count_array );
 191
 192         fprintf( stderr, "done.\n" );
 193     }
 194
 195     close_stream( fp );
 196 }
 197
 198
 199 void scan_seq( char *seq, size_t seq_len, uint *count_array )
 200 {
 201     /* Martin A. Hansen, September 2008 */
 202
 203     /* Run a sliding window over a given sequence. The window */
 204     /* consists of a list where new blocks of 4 nucleotides */
 205     /* are pushed onto one end while at the same time old */
 206     /* blocks are popped from the other end. The number of */
 207     /* in the list is determined by the maximum seperator. */
 208     /* Everytime we have a full window, the window is scanned */
 209     /* for motifs. */
 210
 211     bitblock *block      = NULL;
 212     size_t    b_count    = 0;
 213     ushort    n_count    = 0;
 214     size_t    i          = 0;
 215     uchar     bin        = 0;
 216     bool      first_node = TRUE;
 217     node_sl *new_node    = NULL;
 218     node_sl *old_node    = NULL;
 219     list_sl *list        = list_sl_new();
 220
 221     for ( i = 0; seq[ i ]; i++ )
 222     {
 223         bin <<= BITS_IN_NT;
 224
 225         switch( seq[ i ] )
 226         {
 227             case 'A': case 'a': add_A( bin ); break;
 228             case 'T': case 't': add_T( bin ); break;
 229             case 'C': case 'c': add_C( bin ); break;
 230             case 'G': case 'g': add_G( bin ); break;
 231             default: n_count = BLOCK_SIZE_NT; break;
 232         }
 233
 234         if ( i > BLOCK_SIZE_NT - 2 )
 235         {
 236             b_count++;
 237
 238             block      = bitblock_new();
 239             block->bin = bin;
 240
 241             if ( n_count > 0 )
 242             {
 243                  block->hasN = TRUE;
 244                  n_count--;
 245             }
 246
 247             new_node      = node_sl_new();
 248             new_node->val = block;
 249
 250             if ( first_node ) {
 251                 list_sl_add_beg( &list, &new_node );
 252             } else {
 253                 list_sl_add_after( &old_node, &new_node );
 254             }
 255
 256             old_node = new_node;
 257
 258             first_node = FALSE;
 259
 260             if ( b_count > BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
 261             {
 262                 // bitblock_list_print( list );   /* DEBUG */
 263
 264                 scan_list( list, count_array );
 265
 266                 list_sl_remove_beg( &list );
 267             }
 268         }
 269     }
 270
 271     /* if the list is shorter than BLOCK_SPACE_MAX + BLOCK_SIZE_NT */
 272     if ( b_count <= BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
 273     {
 274         // bitblock_list_print( list );  /* DEBUG */
 275
 276         scan_list( list, count_array );
 277     }
 278
 279     list_sl_destroy( &list );
 280 }
 281
 282
 283 void scan_list( list_sl *list, uint *count_array )
 284 {
 285     /* Martin A. Hansen, September 2008 */
 286
 287     /* Scan a list of blocks for biparite motifs by creating */
 288     /* a binary motif consisting of two blocks of 4 nucleotides */
 289     /* along with the distance separating them. Motifs containing */
 290     /* N's are skipped. */
 291
 292     node_sl  *first_node = NULL;
 293     node_sl  *next_node  = NULL;
 294     bitblock *block1     = NULL;
 295     bitblock *block2     = NULL;
 296     int       i          = 0;
 297     ushort    dist       = 0;
 298     uint      motif_bin  = 0;
 299
 300 //    bitblock_list_print( list );
 301
 302     first_node = list->first;
 303
 304     block1 = ( bitblock * ) first_node->val;
 305
 306     if ( ! block1->hasN )
 307     {
 308         next_node = first_node->next;
 309
 310         for ( i = 0; i < BLOCK_SIZE_NT - 1; i++ ) {
 311             next_node = next_node->next;
 312         }
 313
 314         for ( next_node = next_node; next_node != NULL; next_node = next_node->next )
 315         {
 316             block2 = ( bitblock * ) next_node->val;
 317
 318 //             printf( "block1: %s  block2: %s   dist: %d\n", bin2dna[ block1->bin ], bin2dna[ block2->bin ], dist ); /* DEBUG */
 319
 320             if ( ! block2->hasN )
 321             {
 322                 motif_bin = blocks2motif( block1->bin, block2->bin, dist );
 323
 324                 // motif_print( motif_bin, 0 ); /* DEBUG */
 325                 // bitblock_list_print( list ); /* DEBUG */
 326
 327                 count_array[ motif_bin ]++;
 328             }
 329
 330             dist++;
 331         }
 332     }
 333 }
 334
 335
 336 bitblock *bitblock_new()
 337 {
 338     /* Martin A. Hansen, September 2008 */
 339
 340     /* Initializes a new empty bitblock. */
 341
 342     bitblock *new_block = NULL;
 343
 344     new_block = mem_get( sizeof( bitblock ) );
 345
 346     new_block->bin  = 0;
 347     new_block->hasN = FALSE;
 348
 349     return new_block;
 350 }
 351
 352
 353 uint blocks2motif( uchar bin1, uchar bin2, ushort dist )
 354 {
 355     /* Martin A. Hansen, September 2008 */
 356
 357     /* Given two binary encoded tetra nuceotide blocks, */
 358     /* and the distance separating this, create a binary */
 359     /* bipartite motif. */
 360
 361     uint motif = 0;
 362
 363     motif |= bin1;
 364
 365     motif <<= sizeof( uchar ) * BITS_IN_BYTE;
 366
 367     motif |= bin2;
 368
 369     motif <<= sizeof( uchar ) * BITS_IN_BYTE;
 370
 371     motif |= dist;
 372
 373 //    motif_print( motif, 0 ); /* DEBUG */
 374
 375     return motif;
 376 }
 377
 378
 379 void count_array_print( uint *count_array, size_t nmemb, size_t cutoff )
 380 {
 381     /* Martin A. Hansen, Seqptember 2008. */
 382
 383     /* Print all bipartite motifs in count_array as */
 384     /* tabular output. */
 385
 386     uint i     = 0;
 387     uint motif = 0;
 388     uint count = 0;
 389
 390     for ( i = 0; i < nmemb; i++ )
 391     {
 392         motif = i;
 393         count = count_array[ i ];
 394
 395         if ( count >= cutoff ) {
 396             motif_print( motif, count );
 397         }
 398     }
 399 }
 400
 401
 402 void motif_print( uint motif, uint count )
 403 {
 404     /* Martin A. Hansen, September 2008 */
 405
 406     /* Converts a binary encoded bipartite motif */
 407     /* into DNA and output the motif, distance and */
 408     /* count seperated by tabs: */
 409     /* BLOCK1 \t BLOCK2 \t DIST \t COUNT */
 410
 411     uchar  bin1 = 0;
 412     uchar  bin2 = 0;
 413     ushort dist = 0;
 414
 415     // printf( "%d\t", motif ); /* DEBUG */
 416
 417     dist = ( ushort ) motif & BLOCK_MASK;
 418
 419     motif >>= sizeof( uchar ) * BITS_IN_BYTE;
 420
 421     bin2 = ( uchar ) motif;
 422
 423     motif >>= sizeof( uchar ) * BITS_IN_BYTE;
 424
 425     bin1 = ( uchar ) motif;
 426
 427     printf( "%s\t%s\t%d\t%d\n", bin2dna[ bin1 ], bin2dna[ bin2 ], dist, count );
 428 }
 429
 430
 431 void bitblock_list_print( list_sl *list )
 432 {
 433     /* Martin A. Hansen, September 2008 */
 434
 435     /* Debug function to print all blocks in a list. */
 436
 437     node_sl *node = NULL;
 438
 439     printf( "\nbitblock_list_print:\n" );
 440
 441     for ( node = list->first; node != NULL; node = node->next ) {
 442         bitblock_print( ( bitblock * ) node->val );
 443     }
 444 }
 445
 446
 447 void bitblock_print( bitblock *out )
 448 {
 449     /* Martin A. Hansen, September 2008 */
 450
 451     /* Debug function to print a given block. */
 452
 453     printf( "bin: %d  dna: %s   hasN: %d\n", out->bin, bin2dna[ ( int ) out->bin ], out->hasN );
 454 }
 455
 456
 457 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> UNIT TESTS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 458
 459
 460 void run_tests()
 461 {
 462     fprintf( stderr, "Running tests\n" );
 463
 464     test_count_array_new();
 465     test_bitblock_new();
 466     test_bitblock_print();
 467     test_bitblock_list_print();
 468     test_scan_seq();
 469     test_blocks2motif();
 470
 471     fprintf( stderr, "All tests OK\n" );
 472 }
 473
 474
 475 void test_count_array_new()
 476 {
 477     fprintf( stderr, "   Running test_count_array_new ... " );
 478
 479     size_t  i     = 0;
 480     size_t  nmemb = 128;
 481     uint   *array = NULL;
 482
 483     array = count_array_new( nmemb );
 484
 485     for ( i = 0; i < nmemb; i++ ) {
 486         assert( array[ i ] == 0 );
 487     }
 488
 489     mem_free( &array );
 490
 491     fprintf( stderr, "done.\n" );
 492 }
 493
 494
 495 void test_bitblock_new()
 496 {
 497     fprintf( stderr, "   Running test_bitblock_new ... " );
 498
 499     bitblock *new_block = bitblock_new();
 500
 501     assert( new_block->bin  == 0 );
 502     assert( new_block->hasN == FALSE );
 503
 504     fprintf( stderr, "done.\n" );
 505 }
 506
 507
 508 void test_bitblock_print()
 509 {
 510     fprintf( stderr, "   Running test_bitblock_print ... " );
 511
 512     bitblock *new_block = bitblock_new();
 513
 514     new_block->bin  = 7;
 515     new_block->hasN = TRUE;
 516
 517 //    bitblock_print( new_block );
 518
 519     fprintf( stderr, "done.\n");
 520 }
 521
 522
 523 void test_bitblock_list_print()
 524 {
 525     fprintf( stderr, "   Running test_bitblock_list_print ... " );
 526
 527     list_sl  *list   = list_sl_new();
 528     node_sl  *node1  = node_sl_new();
 529     node_sl  *node2  = node_sl_new();
 530     node_sl  *node3  = node_sl_new();
 531     bitblock *block1 = bitblock_new();
 532     bitblock *block2 = bitblock_new();
 533     bitblock *block3 = bitblock_new();
 534
 535     block1->bin  = 0;
 536     block1->hasN = TRUE;
 537     block2->bin  = 1;
 538     block2->hasN = TRUE;
 539     block3->bin  = 2;
 540     block3->hasN = TRUE;
 541
 542     node1->val  = block1;
 543     node2->val  = block2;
 544     node3->val  = block3;
 545
 546     list_sl_add_beg( &list, &node1 );
 547     list_sl_add_beg( &list, &node2 );
 548     list_sl_add_beg( &list, &node3 );
 549
 550     // bitblock_list_print( list );
 551
 552     fprintf( stderr, "done.\n" );
 553 }
 554
 555
 556 void test_scan_seq()
 557 {
 558     fprintf( stderr, "   Running test_scan_seq ... " );
 559
 560     //char   *seq       = "AAAANTCGGCTNGGGG";
 561     //char   *seq         = "AAAATCGGCTGGGG";
 562     char   *seq         = "AAAAAAAAAAAAAAG";
 563     size_t  seq_len     = strlen( seq );
 564     size_t  nmemb       = 1 << 5;
 565
 566     uint   *count_array = count_array_new( nmemb );
 567
 568     scan_seq( seq, seq_len, count_array );
 569
 570     // count_array_print( count_array, nmemb, 1 );   /* DEBUG */
 571
 572     fprintf( stderr, "done.\n" );
 573 }
 574
 575
 576 static void test_blocks2motif()
 577 {
 578     fprintf( stderr, "   Running test_blocks2motif ... " );
 579
 580     uchar  bin1  = 4;
 581     uchar  bin2  = 3;
 582     ushort dist  = 256;
 583     uint   motif = 0;
 584
 585     motif = blocks2motif( bin1, bin2, dist );
 586
 587 //    printf( "motif: %d\n", motif );
 588
 589     fprintf( stderr, "done.\n");
 590 }
 591
 592
 593 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 594
 595