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, uint *output_array );
  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         printf( "SEQ_NAME: %s\n", entry->seq_name );
 178
 179         rescan_seq( entry->seq, entry->seq_len, count_array, cutoff );
 180
 181         fprintf( stderr, "done.\n" );
 182     }
 183
 184     close_stream( fp );
 185 }
 186
 187
 188 void scan_seq( char *seq, size_t seq_len, uint *count_array )
 189 {
 190     /* Martin A. Hansen, September 2008 */
 191
 192     /* Run a sliding window over a given sequence. The window */
 193     /* consists of a list where new blocks of 4 nucleotides */
 194     /* are pushed onto one end while at the same time old */
 195     /* blocks are popped from the other end. The number of */
 196     /* in the list is determined by the maximum seperator. */
 197     /* Everytime we have a full window, the window is scanned */
 198     /* for motifs. */
 199
 200     bitblock *block      = NULL;
 201     size_t    b_count    = 0;
 202     ushort    n_count    = 0;
 203     size_t    i          = 0;
 204     uchar     bin        = 0;
 205     bool      first_node = TRUE;
 206     node_sl *new_node    = NULL;
 207     node_sl *old_node    = NULL;
 208     list_sl *list        = list_sl_new();
 209
 210     for ( i = 0; seq[ i ]; i++ )
 211     {
 212         bin <<= BITS_IN_NT;
 213
 214         switch( seq[ i ] )
 215         {
 216             case 'A': case 'a': add_A( bin ); break;
 217             case 'T': case 't': add_T( bin ); break;
 218             case 'C': case 'c': add_C( bin ); break;
 219             case 'G': case 'g': add_G( bin ); break;
 220             default: n_count = BLOCK_SIZE_NT; break;
 221         }
 222
 223         if ( i > BLOCK_SIZE_NT - 2 )
 224         {
 225             b_count++;
 226
 227             block      = bitblock_new();
 228             block->bin = bin;
 229
 230             if ( n_count > 0 )
 231             {
 232                  block->hasN = TRUE;
 233                  n_count--;
 234             }
 235
 236             new_node      = node_sl_new();
 237             new_node->val = block;
 238
 239             if ( first_node ) {
 240                 list_sl_add_beg( &list, &new_node );
 241             } else {
 242                 list_sl_add_after( &old_node, &new_node );
 243             }
 244
 245             old_node = new_node;
 246
 247             first_node = FALSE;
 248
 249             if ( b_count > BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
 250             {
 251                 // bitblock_list_print( list );   /* DEBUG */
 252
 253                 scan_list( list, count_array );
 254
 255                 mem_free( &list->first->val );
 256
 257                 list_sl_remove_beg( &list );
 258             }
 259         }
 260     }
 261
 262     /* if the list is shorter than BLOCK_SPACE_MAX + BLOCK_SIZE_NT */
 263     if ( b_count <= BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
 264     {
 265         // bitblock_list_print( list );  /* DEBUG */
 266
 267         scan_list( list, count_array );
 268     }
 269
 270     list_sl_destroy( &list );
 271 }
 272
 273
 274 void rescan_seq( char *seq, size_t seq_len, uint *count_array, size_t cutoff )
 275 {
 276     /* Martin A. Hansen, September 2008 */
 277
 278     /* Run a sliding window over a given sequence. The window */
 279     /* consists of a list where new blocks of 4 nucleotides */
 280     /* are pushed onto one end while at the same time old */
 281     /* blocks are popped from the other end. The number of */
 282     /* in the list is determined by the maximum seperator. */
 283     /* Everytime we have a full window, the window is scanned */
 284     /* for motifs. */
 285
 286     bitblock *block      = NULL;
 287     size_t    b_count    = 0;
 288     ushort    n_count    = 0;
 289     size_t    i          = 0;
 290     uchar     bin        = 0;
 291     bool      first_node = TRUE;
 292     node_sl *new_node    = NULL;
 293     node_sl *old_node    = NULL;
 294     list_sl *list        = list_sl_new();
 295     uint *output_array   = NULL;
 296
 297     output_array         = mem_get_zero( sizeof( uint ) * ( seq_len + 1 ) );
 298
 299     for ( i = 0; seq[ i ]; i++ )
 300     {
 301         bin <<= BITS_IN_NT;
 302
 303         switch( seq[ i ] )
 304         {
 305             case 'A': case 'a': add_A( bin ); break;
 306             case 'T': case 't': add_T( bin ); break;
 307             case 'C': case 'c': add_C( bin ); break;
 308             case 'G': case 'g': add_G( bin ); break;
 309             default: n_count = BLOCK_SIZE_NT; break;
 310         }
 311
 312         if ( i > BLOCK_SIZE_NT - 2 )
 313         {
 314             b_count++;
 315
 316             block      = bitblock_new();
 317             block->bin = bin;
 318
 319             if ( n_count > 0 )
 320             {
 321                  block->hasN = TRUE;
 322                  n_count--;
 323             }
 324
 325             new_node      = node_sl_new();
 326             new_node->val = block;
 327
 328             if ( first_node ) {
 329                 list_sl_add_beg( &list, &new_node );
 330             } else {
 331                 list_sl_add_after( &old_node, &new_node );
 332             }
 333
 334             old_node = new_node;
 335
 336             first_node = FALSE;
 337
 338             if ( b_count > BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
 339             {
 340                 // bitblock_list_print( list );   /* DEBUG */
 341
 342                 rescan_list( list, count_array, i, cutoff, output_array );
 343
 344                 mem_free( &list->first->val );
 345
 346                 list_sl_remove_beg( &list );
 347             }
 348         }
 349     }
 350
 351     /* if the list is shorter than BLOCK_SPACE_MAX + BLOCK_SIZE_NT */
 352     if ( b_count <= BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
 353     {
 354         // bitblock_list_print( list );  /* DEBUG */
 355
 356         rescan_list( list, count_array, i, cutoff, output_array );
 357     }
 358
 359     list_sl_destroy( &list );
 360
 361     for ( i = 0; i < seq_len; i++ ) {
 362         printf( "%zu\t%u\n", i, output_array[ i ] );
 363     }
 364
 365     free( output_array );
 366 }
 367
 368
 369 void scan_list( list_sl *list, uint *count_array )
 370 {
 371     /* Martin A. Hansen, September 2008 */
 372
 373     /* Scan a list of blocks for biparite motifs by creating */
 374     /* a binary motif consisting of two blocks of 4 nucleotides */
 375     /* along with the distance separating them. Motifs containing */
 376     /* N's are skipped. */
 377
 378     node_sl  *first_node = NULL;
 379     node_sl  *next_node  = NULL;
 380     bitblock *block1     = NULL;
 381     bitblock *block2     = NULL;
 382     int       i          = 0;
 383     ushort    dist       = 0;
 384     uint      motif_bin  = 0;
 385
 386 //    bitblock_list_print( list );
 387
 388     first_node = list->first;
 389
 390     block1 = ( bitblock * ) first_node->val;
 391
 392     if ( ! block1->hasN )
 393     {
 394         next_node = first_node->next;
 395
 396         for ( i = 0; i < BLOCK_SIZE_NT - 1; i++ ) {
 397             next_node = next_node->next;
 398         }
 399
 400         for ( next_node = next_node; next_node != NULL; next_node = next_node->next )
 401         {
 402             block2 = ( bitblock * ) next_node->val;
 403
 404             if ( ! block2->hasN )
 405             {
 406                 motif_bin = blocks2motif( block1->bin, block2->bin, dist );
 407
 408                 // bitblock_list_print( list ); /* DEBUG */
 409
 410                 count_array[ motif_bin ]++;
 411             }
 412
 413             dist++;
 414         }
 415     }
 416 }
 417
 418
 419 void rescan_list( list_sl *list, uint *count_array, size_t pos, size_t cutoff, uint *output_array )
 420 {
 421     /* Martin A. Hansen, September 2008 */
 422
 423     /* Scan a list of blocks for biparite motifs by creating */
 424     /* a binary motif consisting of two blocks of 4 nucleotides */
 425     /* along with the distance separating them. Motifs containing */
 426     /* N's are skipped. */
 427
 428     node_sl  *first_node = NULL;
 429     node_sl  *next_node  = NULL;
 430     bitblock *block1     = NULL;
 431     bitblock *block2     = NULL;
 432     int       i          = 0;
 433     int       k          = 0;
 434     ushort    dist       = 0;
 435     uint      motif_bin  = 0;
 436     uint      j          = 0;
 437     uint      count      = 0;
 438
 439 //    bitblock_list_print( list );
 440
 441     first_node = list->first;
 442
 443     block1 = ( bitblock * ) first_node->val;
 444
 445     if ( ! block1->hasN )
 446     {
 447         next_node = first_node->next;
 448
 449         for ( i = 0; i < BLOCK_SIZE_NT - 1; i++ )
 450         {
 451             next_node = next_node->next;
 452
 453             j++;
 454         }
 455
 456         for ( next_node = next_node; next_node != NULL; next_node = next_node->next )
 457         {
 458             block2 = ( bitblock * ) next_node->val;
 459
 460             if ( ! block2->hasN )
 461             {
 462                 motif_bin = blocks2motif( block1->bin, block2->bin, dist );
 463
 464                 // bitblock_list_print( list ); /* DEBUG */
 465
 466                 count = count_array[ motif_bin ];
 467
 468                 if ( count > cutoff )
 469                 {
 470                     // printf( "%zu\t%u\t%u\n", pos + j, motif_bin, count );
 471
 472                     for ( k = 0; k < BLOCK_SIZE_NT - 1; k++ )
 473                     {
 474                         output_array[ pos - j + k - BLOCK_SIZE_NT ]    += count;
 475                         output_array[ pos + k - dist - BLOCK_SIZE_NT ] += count;
 476                     }
 477                 }
 478             }
 479
 480             dist++;
 481
 482             j++;
 483         }
 484     }
 485 }
 486
 487
 488 bitblock *bitblock_new()
 489 {
 490     /* Martin A. Hansen, September 2008 */
 491
 492     /* Initializes a new empty bitblock. */
 493
 494     bitblock *new_block = NULL;
 495
 496     new_block = mem_get( sizeof( bitblock ) );
 497
 498     new_block->bin  = 0;
 499     new_block->hasN = FALSE;
 500
 501     return new_block;
 502 }
 503
 504
 505 uint blocks2motif( uchar bin1, uchar bin2, ushort dist )
 506 {
 507     /* Martin A. Hansen, September 2008 */
 508
 509     /* Given two binary encoded tetra nuceotide blocks, */
 510     /* and the distance separating this, create a binary */
 511     /* bipartite motif. */
 512
 513     uint motif = 0;
 514
 515     motif |= bin1;
 516
 517     motif <<= sizeof( uchar ) * BITS_IN_BYTE;
 518
 519     motif |= bin2;
 520
 521     motif <<= sizeof( uchar ) * BITS_IN_BYTE;
 522
 523     motif |= dist;
 524
 525     return motif;
 526 }
 527
 528
 529 void count_array_print( uint *count_array, size_t nmemb, size_t cutoff )
 530 {
 531     /* Martin A. Hansen, Seqptember 2008. */
 532
 533     /* Print all bipartite motifs in count_array as */
 534     /* tabular output. */
 535
 536     uint i     = 0;
 537     uint motif = 0;
 538     uint count = 0;
 539
 540     for ( i = 0; i < nmemb; i++ )
 541     {
 542         motif = i;
 543         count = count_array[ i ];
 544
 545         if ( count >= cutoff ) {
 546             printf( "%u\t%u\n", motif, count );
 547         }
 548     }
 549 }
 550
 551
 552 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> UNIT TESTS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 553
 554
 555 void run_tests()
 556 {
 557     fprintf( stderr, "Running tests\n" );
 558
 559     test_count_array_new();
 560     test_bitblock_new();
 561     test_bitblock_print();
 562     test_bitblock_list_print();
 563     test_scan_seq();
 564     test_blocks2motif();
 565
 566     fprintf( stderr, "All tests OK\n" );
 567 }
 568
 569
 570 void test_count_array_new()
 571 {
 572     fprintf( stderr, "   Running test_count_array_new ... " );
 573
 574     size_t  i     = 0;
 575     size_t  nmemb = 128;
 576     uint   *array = NULL;
 577
 578     array = count_array_new( nmemb );
 579
 580     for ( i = 0; i < nmemb; i++ ) {
 581         assert( array[ i ] == 0 );
 582     }
 583
 584     mem_free( &array );
 585
 586     fprintf( stderr, "done.\n" );
 587 }
 588
 589
 590 void test_bitblock_new()
 591 {
 592     fprintf( stderr, "   Running test_bitblock_new ... " );
 593
 594     bitblock *new_block = bitblock_new();
 595
 596     assert( new_block->bin  == 0 );
 597     assert( new_block->hasN == FALSE );
 598
 599     fprintf( stderr, "done.\n" );
 600 }
 601
 602
 603 void test_bitblock_print()
 604 {
 605     fprintf( stderr, "   Running test_bitblock_print ... " );
 606
 607     bitblock *new_block = bitblock_new();
 608
 609     new_block->bin  = 7;
 610     new_block->hasN = TRUE;
 611
 612 //    bitblock_print( new_block );
 613
 614     fprintf( stderr, "done.\n");
 615 }
 616
 617
 618 void test_bitblock_list_print()
 619 {
 620     fprintf( stderr, "   Running test_bitblock_list_print ... " );
 621
 622     list_sl  *list   = list_sl_new();
 623     node_sl  *node1  = node_sl_new();
 624     node_sl  *node2  = node_sl_new();
 625     node_sl  *node3  = node_sl_new();
 626     bitblock *block1 = bitblock_new();
 627     bitblock *block2 = bitblock_new();
 628     bitblock *block3 = bitblock_new();
 629
 630     block1->bin  = 0;
 631     block1->hasN = TRUE;
 632     block2->bin  = 1;
 633     block2->hasN = TRUE;
 634     block3->bin  = 2;
 635     block3->hasN = TRUE;
 636
 637     node1->val  = block1;
 638     node2->val  = block2;
 639     node3->val  = block3;
 640
 641     list_sl_add_beg( &list, &node1 );
 642     list_sl_add_beg( &list, &node2 );
 643     list_sl_add_beg( &list, &node3 );
 644
 645     // bitblock_list_print( list );
 646
 647     fprintf( stderr, "done.\n" );
 648 }
 649
 650
 651 void test_scan_seq()
 652 {
 653     fprintf( stderr, "   Running test_scan_seq ... " );
 654
 655     //char   *seq       = "AAAANTCGGCTNGGGG";
 656     //char   *seq         = "AAAATCGGCTGGGG";
 657     char   *seq         = "AAAAAAAAAAAAAAG";
 658     size_t  seq_len     = strlen( seq );
 659     size_t  nmemb       = 1 << 5;
 660
 661     uint   *count_array = count_array_new( nmemb );
 662
 663     scan_seq( seq, seq_len, count_array );
 664
 665     // count_array_print( count_array, nmemb, 1 );   /* DEBUG */
 666
 667     fprintf( stderr, "done.\n" );
 668 }
 669
 670
 671 static void test_blocks2motif()
 672 {
 673     fprintf( stderr, "   Running test_blocks2motif ... " );
 674
 675     uchar  bin1  = 4;
 676     uchar  bin2  = 3;
 677     ushort dist  = 256;
 678     uint   motif = 0;
 679
 680     motif = blocks2motif( bin1, bin2, dist );
 681
 682 //    printf( "motif: %d\n", motif );
 683
 684     fprintf( stderr, "done.\n");
 685 }
 686
 687
 688 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 689
 690