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            10000                            /* minimum number of motifs in output. */
  17
  18 #define add_A( c )              /* add 00 to the rightmost two bits of bin (i.e. do nothing). */
  19 #define add_T( c ) ( c |= 3 )   /* add 11 on the rightmost two bits of c. */
  20 #define add_C( c ) ( c |= 1 )   /* add 01 on the rightmost two bits of c. */
  21 #define add_G( c ) ( c |= 2 )   /* add 10 on the rightmost two bits of c. */
  22
  23 /* Structure that will hold one tetra nucleotide block. */
  24 struct _bitblock
  25 {
  26     uchar bin;   /* Tetra nucleotide binary encoded. */
  27     bool hasN;   /* Flag indicating any N's in the block. */
  28 };
  29
  30 typedef struct _bitblock bitblock;
  31
  32 /* Function declarations. */
  33 void      run_scan( int argc, char *argv[] );
  34 void      print_usage();
  35 void      scan_file( char *file, seq_entry *entry, uint *count_array );
  36 uint     *count_array_new( size_t nmemb );
  37 void      scan_seq( char *seq, size_t seq_len, uint *count_array );
  38 void      scan_list( list_sl *list, uint *count_array );
  39 bitblock *bitblock_new();
  40 uint      blocks2motif( uchar bin1, uchar bin2, ushort dist );
  41 void      count_array_print( uint *count_array, size_t nmemb, size_t cutoff );
  42 void      bitblock_list_print( list_sl *list );
  43 void      bitblock_print( bitblock *out );
  44
  45 /* Unit test declarations. */
  46 static void run_tests();
  47 static void test_count_array_new();
  48 static void test_bitblock_new();
  49 static void test_bitblock_print();
  50 static void test_bitblock_list_print();
  51 static void test_scan_seq();
  52 static void test_blocks2motif();
  53
  54
  55 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MAIN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
  56
  57
  58 int main( int argc, char *argv[] )
  59 {
  60     run_tests();
  61
  62     if ( argc == 1 ) {
  63         print_usage();
  64     }
  65
  66     run_scan( argc, argv );
  67
  68     return EXIT_SUCCESS;
  69 }
  70
  71
  72 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FUNCTIONS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
  73
  74
  75 void print_usage()
  76 {
  77     /* Martin A. Hansen, September 2008 */
  78
  79     /* Print usage and exit if no files in argument. */
  80
  81     fprintf( stderr,
  82         "Usage: bipartite_scam <FASTA file(s)> > result.csv\n"
  83     );
  84
  85     exit( EXIT_SUCCESS );
  86 }
  87
  88
  89 void run_scan( int argc, char *argv[] )
  90 {
  91     /* Martin A. Hansen, September 2008 */
  92
  93     /* For each file in argv scan the file for */
  94     /* bipartite motifs and output the motifs */
  95     /* and their count. */
  96
  97     char      *file        = NULL;
  98     int        i           = 0;
  99     seq_entry *entry       = NULL;
 100     uint      *count_array = NULL;
 101 //    size_t    new_nmemb    = 0;
 102
 103     count_array = count_array_new( COUNT_ARRAY_NMEMB );
 104
 105     entry = seq_new( MAX_SEQ_NAME, MAX_SEQ );
 106
 107     for ( i = 1; i < argc; i++ )
 108     {
 109         file = argv[ i ];
 110
 111         fprintf( stderr, "Scanning file: %s\n", file );
 112         scan_file( file, entry, count_array );
 113         fprintf( stderr, "done.\n" );
 114     }
 115
 116     fprintf( stderr, "Printing motifs: ... " );
 117     count_array_print( count_array, COUNT_ARRAY_NMEMB, CUTOFF );
 118     fprintf( stderr, "done.\n" );
 119
 120     seq_destroy( entry );
 121
 122     mem_free( &count_array );
 123 }
 124
 125
 126 uint *count_array_new( size_t nmemb )
 127 {
 128     /* Martin A. Hansen, September 2008 */
 129
 130     /* Initialize a new zeroed uint array of nmemb objects. */
 131
 132     uint *array = NULL;
 133
 134     assert( nmemb > 0 );
 135
 136     array = mem_get_zero( nmemb * sizeof( uint ) );
 137
 138     return array;
 139 }
 140
 141
 142 void scan_file( char *file, seq_entry *entry, uint *count_array )
 143 {
 144     /* Martin A. Hansen, September 2008 */
 145
 146     /* Scan all FASTA entries of a file in both */
 147     /* sense and anti-sense directions. */
 148
 149     FILE *fp = read_open( file );
 150
 151     while ( fasta_get_entry( fp, &entry ) == TRUE )
 152     {
 153         fprintf( stderr, "   Scanning: %s (%zu nt) ... ", entry->seq_name, entry->seq_len );
 154
 155         scan_seq( entry->seq, entry->seq_len, count_array );
 156
 157         fprintf( stderr, "done.\n" );
 158     }
 159
 160     close_stream( fp );
 161 }
 162
 163
 164 void scan_seq( char *seq, size_t seq_len, uint *count_array )
 165 {
 166     /* Martin A. Hansen, September 2008 */
 167
 168     /* Run a sliding window over a given sequence. The window */
 169     /* consists of a list where new blocks of 4 nucleotides */
 170     /* are pushed onto one end while at the same time old */
 171     /* blocks are popped from the other end. The number of */
 172     /* in the list is determined by the maximum seperator. */
 173     /* Everytime we have a full window, the window is scanned */
 174     /* for motifs. */
 175
 176     bitblock *block      = NULL;
 177     size_t    b_count    = 0;
 178     ushort    n_count    = 0;
 179     size_t    i          = 0;
 180     uchar     bin        = 0;
 181     bool      first_node = TRUE;
 182     node_sl *new_node    = NULL;
 183     node_sl *old_node    = NULL;
 184     list_sl *list        = list_sl_new();
 185
 186     for ( i = 0; seq[ i ]; i++ )
 187     {
 188         bin <<= BITS_IN_NT;
 189
 190         switch( seq[ i ] )
 191         {
 192             case 'A': case 'a': add_A( bin ); break;
 193             case 'T': case 't': add_T( bin ); break;
 194             case 'C': case 'c': add_C( bin ); break;
 195             case 'G': case 'g': add_G( bin ); break;
 196             default: n_count = BLOCK_SIZE_NT; break;
 197         }
 198
 199         if ( i > BLOCK_SIZE_NT - 2 )
 200         {
 201             b_count++;
 202
 203             block      = bitblock_new();
 204             block->bin = bin;
 205
 206             if ( n_count > 0 )
 207             {
 208                  block->hasN = TRUE;
 209                  n_count--;
 210             }
 211
 212             new_node      = node_sl_new();
 213             new_node->val = block;
 214
 215             if ( first_node ) {
 216                 list_sl_add_beg( &list, &new_node );
 217             } else {
 218                 list_sl_add_after( &old_node, &new_node );
 219             }
 220
 221             old_node = new_node;
 222
 223             first_node = FALSE;
 224
 225             if ( b_count > BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
 226             {
 227                 // bitblock_list_print( list );   /* DEBUG */
 228
 229                 scan_list( list, count_array );
 230
 231                 list_sl_remove_beg( &list );
 232             }
 233         }
 234     }
 235
 236     /* if the list is shorter than BLOCK_SPACE_MAX + BLOCK_SIZE_NT */
 237     if ( b_count <= BLOCK_SPACE_MAX + BLOCK_SIZE_NT )
 238     {
 239         // bitblock_list_print( list );  /* DEBUG */
 240
 241         scan_list( list, count_array );
 242     }
 243
 244     list_sl_destroy( &list );
 245 }
 246
 247
 248 void scan_list( list_sl *list, uint *count_array )
 249 {
 250     /* Martin A. Hansen, September 2008 */
 251
 252     /* Scan a list of blocks for biparite motifs by creating */
 253     /* a binary motif consisting of two blocks of 4 nucleotides */
 254     /* along with the distance separating them. Motifs containing */
 255     /* N's are skipped. */
 256
 257     node_sl  *first_node = NULL;
 258     node_sl  *next_node  = NULL;
 259     bitblock *block1     = NULL;
 260     bitblock *block2     = NULL;
 261     int       i          = 0;
 262     ushort    dist       = 0;
 263     uint      motif_bin  = 0;
 264
 265 //    bitblock_list_print( list );
 266
 267     first_node = list->first;
 268
 269     block1 = ( bitblock * ) first_node->val;
 270
 271     if ( ! block1->hasN )
 272     {
 273         next_node = first_node->next;
 274
 275         for ( i = 0; i < BLOCK_SIZE_NT - 1; i++ ) {
 276             next_node = next_node->next;
 277         }
 278
 279         for ( next_node = next_node; next_node != NULL; next_node = next_node->next )
 280         {
 281             block2 = ( bitblock * ) next_node->val;
 282
 283             if ( ! block2->hasN )
 284             {
 285                 motif_bin = blocks2motif( block1->bin, block2->bin, dist );
 286
 287                 // bitblock_list_print( list ); /* DEBUG */
 288
 289                 count_array[ motif_bin ]++;
 290             }
 291
 292             dist++;
 293         }
 294     }
 295 }
 296
 297
 298 bitblock *bitblock_new()
 299 {
 300     /* Martin A. Hansen, September 2008 */
 301
 302     /* Initializes a new empty bitblock. */
 303
 304     bitblock *new_block = NULL;
 305
 306     new_block = mem_get( sizeof( bitblock ) );
 307
 308     new_block->bin  = 0;
 309     new_block->hasN = FALSE;
 310
 311     return new_block;
 312 }
 313
 314
 315 uint blocks2motif( uchar bin1, uchar bin2, ushort dist )
 316 {
 317     /* Martin A. Hansen, September 2008 */
 318
 319     /* Given two binary encoded tetra nuceotide blocks, */
 320     /* and the distance separating this, create a binary */
 321     /* bipartite motif. */
 322
 323     uint motif = 0;
 324
 325     motif |= bin1;
 326
 327     motif <<= sizeof( uchar ) * BITS_IN_BYTE;
 328
 329     motif |= bin2;
 330
 331     motif <<= sizeof( uchar ) * BITS_IN_BYTE;
 332
 333     motif |= dist;
 334
 335     return motif;
 336 }
 337
 338
 339 void count_array_print( uint *count_array, size_t nmemb, size_t cutoff )
 340 {
 341     /* Martin A. Hansen, Seqptember 2008. */
 342
 343     /* Print all bipartite motifs in count_array as */
 344     /* tabular output. */
 345
 346     uint i     = 0;
 347     uint motif = 0;
 348     uint count = 0;
 349
 350     for ( i = 0; i < nmemb; i++ )
 351     {
 352         motif = i;
 353         count = count_array[ i ];
 354
 355         if ( count >= cutoff ) {
 356             printf( "%u\t%u\n", motif, count );
 357         }
 358     }
 359 }
 360
 361
 362 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> UNIT TESTS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 363
 364
 365 void run_tests()
 366 {
 367     fprintf( stderr, "Running tests\n" );
 368
 369     test_count_array_new();
 370     test_bitblock_new();
 371     test_bitblock_print();
 372     test_bitblock_list_print();
 373     test_scan_seq();
 374     test_blocks2motif();
 375
 376     fprintf( stderr, "All tests OK\n" );
 377 }
 378
 379
 380 void test_count_array_new()
 381 {
 382     fprintf( stderr, "   Running test_count_array_new ... " );
 383
 384     size_t  i     = 0;
 385     size_t  nmemb = 128;
 386     uint   *array = NULL;
 387
 388     array = count_array_new( nmemb );
 389
 390     for ( i = 0; i < nmemb; i++ ) {
 391         assert( array[ i ] == 0 );
 392     }
 393
 394     mem_free( &array );
 395
 396     fprintf( stderr, "done.\n" );
 397 }
 398
 399
 400 void test_bitblock_new()
 401 {
 402     fprintf( stderr, "   Running test_bitblock_new ... " );
 403
 404     bitblock *new_block = bitblock_new();
 405
 406     assert( new_block->bin  == 0 );
 407     assert( new_block->hasN == FALSE );
 408
 409     fprintf( stderr, "done.\n" );
 410 }
 411
 412
 413 void test_bitblock_print()
 414 {
 415     fprintf( stderr, "   Running test_bitblock_print ... " );
 416
 417     bitblock *new_block = bitblock_new();
 418
 419     new_block->bin  = 7;
 420     new_block->hasN = TRUE;
 421
 422 //    bitblock_print( new_block );
 423
 424     fprintf( stderr, "done.\n");
 425 }
 426
 427
 428 void test_bitblock_list_print()
 429 {
 430     fprintf( stderr, "   Running test_bitblock_list_print ... " );
 431
 432     list_sl  *list   = list_sl_new();
 433     node_sl  *node1  = node_sl_new();
 434     node_sl  *node2  = node_sl_new();
 435     node_sl  *node3  = node_sl_new();
 436     bitblock *block1 = bitblock_new();
 437     bitblock *block2 = bitblock_new();
 438     bitblock *block3 = bitblock_new();
 439
 440     block1->bin  = 0;
 441     block1->hasN = TRUE;
 442     block2->bin  = 1;
 443     block2->hasN = TRUE;
 444     block3->bin  = 2;
 445     block3->hasN = TRUE;
 446
 447     node1->val  = block1;
 448     node2->val  = block2;
 449     node3->val  = block3;
 450
 451     list_sl_add_beg( &list, &node1 );
 452     list_sl_add_beg( &list, &node2 );
 453     list_sl_add_beg( &list, &node3 );
 454
 455     // bitblock_list_print( list );
 456
 457     fprintf( stderr, "done.\n" );
 458 }
 459
 460
 461 void test_scan_seq()
 462 {
 463     fprintf( stderr, "   Running test_scan_seq ... " );
 464
 465     //char   *seq       = "AAAANTCGGCTNGGGG";
 466     //char   *seq         = "AAAATCGGCTGGGG";
 467     char   *seq         = "AAAAAAAAAAAAAAG";
 468     size_t  seq_len     = strlen( seq );
 469     size_t  nmemb       = 1 << 5;
 470
 471     uint   *count_array = count_array_new( nmemb );
 472
 473     scan_seq( seq, seq_len, count_array );
 474
 475     // count_array_print( count_array, nmemb, 1 );   /* DEBUG */
 476
 477     fprintf( stderr, "done.\n" );
 478 }
 479
 480
 481 static void test_blocks2motif()
 482 {
 483     fprintf( stderr, "   Running test_blocks2motif ... " );
 484
 485     uchar  bin1  = 4;
 486     uchar  bin2  = 3;
 487     ushort dist  = 256;
 488     uint   motif = 0;
 489
 490     motif = blocks2motif( bin1, bin2, dist );
 491
 492 //    printf( "motif: %d\n", motif );
 493
 494     fprintf( stderr, "done.\n");
 495 }
 496
 497
 498 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
 499
 500