code_c/Maasha/src/bipartite_scan.c

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