code_c/Maasha/src/align_two_seq.c

   1 #include "common.h"
   2 #include "mem.h"
   3 #include "hash.h"
   4
   5
   6 /*********************************** DECLARATIONS ************************************/
   7
   8
   9 struct _search_space
  10 {
  11     char         *q_seq;
  12     char         *s_seq;
  13     unsigned int  q_beg;
  14     unsigned int  s_beg;
  15     unsigned int  q_end;
  16     unsigned int  s_end;
  17 };
  18
  19 typedef struct _search_space search_space;
  20
  21
  22 struct _match
  23 {
  24     struct _match *next;
  25     unsigned int   q_beg;
  26     unsigned int   s_beg;
  27     unsigned int   len;
  28     float          score;
  29 };
  30
  31 typedef struct _match match;
  32
  33
  34 struct _list
  35 {
  36     struct _list *next;
  37     unsigned int  val;
  38 };
  39
  40 typedef struct _list list;
  41
  42
  43 /*********************************** FUNCTIONS ************************************/
  44
  45
  46 search_space *search_space_new( char *q_seq, char *s_seq, unsigned int q_beg, unsigned int s_beg, unsigned int q_end, unsigned int s_end )
  47 {
  48     search_space *new = NULL;
  49
  50     new = mem_get( sizeof( search_space ) );
  51
  52     assert( q_seq != NULL );
  53     assert( s_seq != NULL );
  54     assert( q_beg >= 0 );
  55     assert( s_beg >= 0 );
  56     assert( q_end > q_beg );
  57     assert( s_end > s_beg );
  58
  59     new->q_seq = q_seq;
  60     new->s_seq = s_seq;
  61     new->q_beg = q_beg;
  62     new->s_beg = s_beg;
  63     new->q_end = q_end;
  64     new->s_end = s_end;
  65
  66     return new;
  67 }
  68
  69
  70 match *match_new( unsigned int q_beg, unsigned int s_beg, unsigned int len )
  71 {
  72     match *new = NULL;
  73
  74     new = mem_get( sizeof( match ) );
  75
  76     new->q_beg = q_beg;
  77     new->s_beg = s_beg;
  78     new->len   = len;
  79     new->next  = NULL;
  80     new->score = 0.0;
  81
  82     return new;
  83 }
  84
  85
  86 void match_add( match **old_ppt, match **new_ppt )
  87 {
  88     match *old = *old_ppt;
  89     match *new = *new_ppt;
  90
  91     new->next = old->next;
  92     old->next = new;
  93 }
  94
  95
  96 void match_print( match *match_pt )
  97 {
  98     printf( "q_beg: %d\ts_beg: %d\tlen: %d\tscore: %f\n", match_pt->q_beg, match_pt->s_beg, match_pt->len, match_pt->score );
  99 }
 100
 101
 102 void matches_print( match *matches )
 103 {
 104     match *m = NULL;
 105
 106     for ( m = matches; m != NULL; m = m->next ) {
 107         match_print( m );
 108     }
 109 }
 110
 111
 112 void match_expand_forward( match *match_pt, search_space *ss )
 113 {
 114     while ( match_pt->q_beg + match_pt->len <= ss->q_end &&
 115             match_pt->s_beg + match_pt->len <= ss->s_end &&
 116             toupper( ss->q_seq[ match_pt->q_beg + match_pt->len ] ) == toupper( ss->s_seq[ match_pt->s_beg + match_pt->len ] ) )
 117     {
 118         match_pt->len++;
 119     }
 120 }
 121
 122
 123 void match_expand_backward( match *match_pt, search_space *ss )
 124 {
 125     while ( match_pt->q_beg > ss->q_beg &&
 126             match_pt->s_beg > ss->q_beg &&
 127             toupper( ss->q_seq[ match_pt->q_beg ] ) == toupper( ss->s_seq[ match_pt->s_beg ] ) )
 128     {
 129         match_pt->q_beg--;
 130         match_pt->s_beg--;
 131         match_pt->len++;
 132     }
 133 }
 134
 135
 136 void match_expand( match *match_pt, search_space *ss )
 137 {
 138     match_expand_forward( match_pt, ss );
 139     match_expand_backward( match_pt, ss );
 140 }
 141
 142
 143 bool match_redundant( match *matches, match *new )
 144 {
 145     // FIXME replace with binary search scheme. (test if matches always are sorted).
 146
 147     match *m = NULL;
 148
 149     for ( m = matches; m != NULL; m = m->next )
 150     {
 151         if ( new->q_beg >= m->q_beg &&
 152              new->s_beg >= m->s_beg &&
 153              new->q_beg + new->len - 1 <= m->q_beg + m->len - 1 &&
 154              new->s_beg + new->len - 1 <= m->s_beg + m->len - 1
 155            )
 156         {
 157             return TRUE;
 158         }
 159     }
 160
 161     return FALSE;
 162 }
 163
 164
 165 list *list_new()
 166 {
 167     list *new = NULL;
 168
 169     new = mem_get( sizeof( list ) );
 170
 171     new->next = NULL;
 172     new->val  = 0;
 173
 174     return new;
 175 }
 176
 177
 178 void list_add( list **old_ppt, list **new_ppt )
 179 {
 180     list *old = *old_ppt;
 181     list *new = *new_ppt;
 182
 183     new->next = old->next;
 184     old->next = new;
 185 }
 186
 187
 188 void list_print( list *list_pt )
 189 {
 190     list *node = NULL;
 191
 192     for ( node = list_pt; node != NULL; node = node->next ) {
 193         printf( "node->val: %d\n", node->val );
 194     }
 195 }
 196
 197
 198 void seq_index( hash **index_ppt, char *seq, unsigned int seq_beg, unsigned int seq_end, unsigned int word_size )
 199 {
 200     hash         *index = *index_ppt;
 201     unsigned int  i     = 0;
 202     char         *word  = NULL;
 203     list         *old   = NULL;
 204     list         *new   = NULL;
 205
 206     assert( 0 <= seq_beg );
 207     assert( seq_beg < seq_end );
 208     assert( word_size > 0 );
 209
 210     word = mem_get_zero( sizeof( char ) * word_size + 1 );
 211
 212     for ( i = seq_beg; i < seq_end - word_size + 2; i++ )
 213     {
 214         memcpy( word, &seq[ i ], word_size );
 215
 216         new = list_new();
 217
 218         new->val = i;
 219
 220         if ( ( old = hash_get( index, word ) ) != NULL ) {
 221             list_add( &old, &new );
 222         } else {
 223             hash_add( index, word, new );
 224         }
 225     }
 226
 227     free( word );
 228 }
 229
 230
 231 void index_print( hash *index )
 232 {
 233     hash_elem *bucket   = NULL;
 234     list      *pos_list = NULL;
 235     list      *node     = NULL;
 236     size_t     i        = 0;
 237
 238     printf( "\ntable_size: %zu mask: %zu elem_count: %zu\n", index->table_size, index->mask, index->nmemb );
 239
 240     for ( i = 0; i < index->table_size; i++ )
 241     {
 242         bucket = index->table[ i ];
 243
 244         while ( bucket != NULL  )
 245         {
 246             printf( "i: %zu   key: %s   val: ", i, bucket->key );
 247
 248             pos_list = bucket->val;
 249
 250             for ( node = pos_list; node != NULL; node = node->next ) {
 251                 printf( "%d,", node->val );
 252             }
 253
 254             printf( "\n" );
 255
 256             bucket = bucket->next;
 257         }
 258     }
 259 }
 260
 261
 262 unsigned int matches_find_s( match **matches_ppt, hash *index, search_space *ss, unsigned int word_size )
 263 {
 264     match        *matches     = *matches_ppt;
 265     unsigned int  s_pos       = 0;
 266     unsigned int  match_count = 0;
 267     char         *word        = NULL;
 268     list         *pos_list    = NULL;
 269     list         *node        = NULL;
 270     match        *new         = NULL;
 271
 272     word = mem_get_zero( sizeof( char ) * word_size + 1 );
 273
 274     for ( s_pos = ss->s_beg; s_pos < ss->s_end - word_size + 2; s_pos++ )
 275     {
 276         memcpy( word, &ss->s_seq[ s_pos ], word_size );
 277
 278         if ( ( pos_list = hash_get( index, word ) ) != NULL )
 279         {
 280             for ( node = pos_list; node != NULL; node = node->next )
 281             {
 282                 new = match_new( node->val, s_pos, word_size );
 283
 284                 if ( ! match_redundant( matches, new ) )
 285                 {
 286                     match_expand( new, ss );
 287
 288                     match_add( &matches, &new );
 289
 290                     match_count++;
 291                 }
 292             }
 293         }
 294     }
 295
 296     printf( "her\n" ); // DEBUG
 297     matches_print( matches );
 298     index_print( index );
 299
 300
 301     free( word );
 302
 303     return match_count;
 304 }
 305
 306
 307 unsigned int matches_find( match **matches_ppt, search_space *ss, unsigned int word_size )
 308 {
 309     match        *matches     = *matches_ppt;
 310     hash         *index       = NULL;
 311     unsigned int  match_count = 0;
 312
 313     assert( word_size > 0 );
 314
 315     if ( ss->q_end - ss->q_beg < ss->s_end - ss->s_beg )
 316     {
 317         seq_index( &index, ss->q_seq, ss->q_beg, ss->q_end, word_size );
 318
 319         match_count = matches_find_s( &matches, index, ss, word_size );
 320     }
 321     else
 322     {
 323         seq_index( &index, ss->s_seq, ss->s_beg, ss->s_end, word_size );
 324
 325         // match_count = matches_find_q( &matches, index, q_seq, s_seq, q_beg, q_end, word_size );
 326     }
 327
 328     // index_destroy( index );
 329
 330     return match_count;
 331 }
 332
 333
 334 /*********************************** UNIT TESTS ************************************/
 335
 336
 337 void test_search_space_new()
 338 {
 339     fprintf( stderr, "   Testing search_space_new ... " );
 340
 341     search_space *new = NULL;
 342
 343     new = search_space_new( "ATCG", "GCTA", 0, 1, 2, 3 );
 344
 345     assert( new->q_seq == "ATCG" );
 346     assert( new->s_seq == "GCTA" );
 347     assert( new->q_beg == 0 );
 348     assert( new->s_beg == 1 );
 349     assert( new->q_end == 2 );
 350     assert( new->s_end == 3 );
 351
 352     fprintf( stderr, "done.\n" );
 353 }
 354
 355
 356 void test_match_new()
 357 {
 358     fprintf( stderr, "   Testing match_new ... " );
 359
 360     unsigned int q_beg = 1;
 361     unsigned int s_beg = 1;
 362     unsigned int len   = 2;
 363     match *new = NULL;
 364
 365     new = match_new( q_beg, s_beg, len );
 366
 367     assert( new->next  == NULL );
 368     assert( new->q_beg == 1 );
 369     assert( new->s_beg == 1 );
 370     assert( new->len   == 2 );
 371     assert( new->score == 0.0 );
 372
 373     fprintf( stderr, "done.\n" );
 374 }
 375
 376
 377 void test_match_expand_forward()
 378 {
 379     fprintf( stderr, "   Testing match_expand_forward ... " );
 380
 381     search_space *new_ss    = NULL;
 382     match        *new_match = NULL;
 383
 384     new_ss    = search_space_new( "ATCGG", "atcg", 0, 0, 4, 3 );
 385     new_match = match_new( 1, 1, 2 );
 386
 387     match_expand_forward( new_match, new_ss );
 388
 389     assert( new_match->q_beg == 1 );
 390     assert( new_match->s_beg == 1 );
 391     assert( new_match->len   == 3 );
 392
 393     fprintf( stderr, "done.\n" );
 394 }
 395
 396
 397 void test_match_expand_backward()
 398 {
 399     fprintf( stderr, "   Testing match_expand_backward ... " );
 400
 401     search_space *new_ss    = NULL;
 402     match        *new_match = NULL;
 403
 404     new_ss    = search_space_new( "ATCGG", "atcg", 0, 0, 4, 3 );
 405     new_match = match_new( 1, 1, 2 );
 406
 407     match_expand_backward( new_match, new_ss );
 408
 409     assert( new_match->q_beg == 0 );
 410     assert( new_match->s_beg == 0 );
 411     assert( new_match->len   == 3 );
 412
 413     fprintf( stderr, "done.\n" );
 414 }
 415
 416
 417 void test_match_expand()
 418 {
 419     fprintf( stderr, "   Testing match_expand ... " );
 420
 421     search_space *new_ss    = NULL;
 422     match        *new_match = NULL;
 423
 424     new_ss    = search_space_new( "ATCGG", "atcg", 0, 0, 4, 3 );
 425     new_match = match_new( 1, 1, 2 );
 426
 427     match_expand( new_match, new_ss );
 428
 429     assert( new_match->q_beg == 0 );
 430     assert( new_match->s_beg == 0 );
 431     assert( new_match->len   == 4 );
 432
 433     fprintf( stderr, "done.\n" );
 434 }
 435
 436
 437 void test_match_redundant_duplicate()
 438 {
 439     fprintf( stderr, "   Testing match_redundant_duplicate ... " );
 440
 441     match *matches = NULL;
 442     match *new     = NULL;
 443
 444     matches = match_new( 0, 1, 2 );
 445     new     = match_new( 0, 1, 2 );
 446
 447     assert( match_redundant( matches, new ) == TRUE );
 448
 449     fprintf( stderr, "done.\n" );
 450 }
 451
 452
 453 void test_match_redundant_upstream()
 454 {
 455     fprintf( stderr, "   Testing match_redundant_upstream ... " );
 456
 457     match *matches = NULL;
 458     match *new     = NULL;
 459
 460     matches = match_new( 5, 5, 2 );
 461     new     = match_new( 4, 4, 3 );
 462
 463     assert( match_redundant( matches, new ) == FALSE );
 464
 465     fprintf( stderr, "done.\n" );
 466 }
 467
 468
 469 void test_match_redundant_downstream()
 470 {
 471     fprintf( stderr, "   Testing match_redundant_downstream ... " );
 472
 473     match *matches = NULL;
 474     match *new     = NULL;
 475
 476     matches = match_new( 5, 5, 2 );
 477     new     = match_new( 6, 6, 2 );
 478
 479     assert( match_redundant( matches, new ) == FALSE );
 480
 481     fprintf( stderr, "done.\n" );
 482 }
 483
 484
 485 void test_list_new()
 486 {
 487     fprintf( stderr, "   Testing list_new ... " );
 488
 489     list *new = NULL;
 490
 491     new = list_new();
 492
 493     assert( new->next == NULL );
 494     assert( new->val  == 0 );
 495
 496     fprintf( stderr, "done.\n" );
 497 }
 498
 499
 500 void test_list_add()
 501 {
 502     fprintf( stderr, "   Testing list_add ... " );
 503
 504     list *old = NULL;
 505     list *new = NULL;
 506
 507     old = list_new();
 508     new = list_new();
 509
 510     old->val = 10;
 511     new->val = 20;
 512
 513     list_add( &old, &new );
 514
 515     assert( old->val       == 10 );
 516     assert( old->next->val == 20 );
 517
 518     fprintf( stderr, "done.\n" );
 519 }
 520
 521
 522 void test_seq_index()
 523 {
 524     fprintf( stderr, "   Testing seq_index ... " );
 525
 526     char         *seq       = "ATCGATCG";
 527     unsigned int  seq_beg   = 0;
 528     unsigned int  seq_end   = strlen( seq ) - 1;
 529     unsigned int  word_size = 2;
 530     hash         *index     = NULL;
 531     unsigned int  size      = 16;
 532     list         *node      = NULL;
 533
 534     index = hash_new( size );
 535
 536     seq_index( &index, seq, seq_beg, seq_end, word_size );
 537
 538 //    index_print( index );
 539
 540     node = hash_get( index, "GA" );
 541
 542     assert( node->val == 3 );
 543
 544     fprintf( stderr, "done.\n" );
 545 }
 546
 547
 548 void test_matches_find_s()
 549 {
 550     fprintf( stderr, "   Testing matches_find_s ... " );
 551
 552     search_space *ss          = NULL;
 553     unsigned int  word_size   = 2;
 554     unsigned int  match_count = 0;
 555     match        *matches     = NULL;
 556     hash         *index       = NULL;
 557     unsigned int  size        = 16;
 558
 559     ss      = search_space_new( "ATCG", "ATCG", 0, 0, 3, 3 );
 560 //    matches = match_new( 0, 0, 0 ); // FIXME: dummy match
 561     matches = mem_get( sizeof( match ) );
 562
 563     index = hash_new( size );
 564
 565     seq_index( &index, ss->q_seq, ss->q_beg, ss->q_end, word_size );
 566
 567     match_count = matches_find_s( &matches, index, ss, word_size );
 568
 569     fprintf( stderr, "done.\n" );
 570 }
 571
 572
 573 void test_matches_find()
 574 {
 575     fprintf( stderr, "   Testing matches_find ... " );
 576
 577     search_space *ss          = NULL;
 578     unsigned int  word_size   = 2;
 579     unsigned int  match_count = 0;
 580     match        *matches     = NULL;
 581
 582     ss = search_space_new( "ATCG", "ATCG", 0, 0, 3, 3 );
 583
 584     match_count = matches_find( &matches, ss, word_size );
 585
 586     fprintf( stderr, "done.\n" );
 587 }
 588
 589
 590 void test_all()
 591 {
 592     fprintf( stderr, "Running all tests:\n" );
 593
 594     test_search_space_new();
 595
 596     test_match_new();
 597     test_match_expand_forward();
 598     test_match_expand_backward();
 599     test_match_expand();
 600     test_match_redundant_duplicate();
 601     test_match_redundant_upstream();
 602     test_match_redundant_downstream();
 603
 604     test_list_new();
 605     test_list_add();
 606
 607     test_seq_index();
 608
 609     // test_matches_find_q();
 610     test_matches_find_s();
 611     test_matches_find();
 612
 613     fprintf( stderr, "Done.\n\n" );
 614 }
 615
 616
 617 /*********************************** MAIN ************************************/
 618
 619
 620 int main()
 621 {
 622     test_all();
 623
 624     return EXIT_SUCCESS;
 625 }
 626
 627
 628 /***********************************************************************/
 629
 630