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