binaries/src/mafft/core/Galign11.c

   1 #include "mltaln.h"
   2 #include "dp.h"
   3
   4 #define DEBUG 0
   5 #define XXXXXXX    0
   6 #define USE_PENALTY_EX  1
   7
   8
   9 #if 1
  10 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
  11 {
  12         char *seq2 = s2[0];
  13         int *intptr = amino_dis[(int)s1[0][i1]];
  14
  15         while( lgth2-- )
  16                 *match++ = intptr[(int)*seq2++];
  17 }
  18 static void match_calc_mtx( int mtx[0x80][0x80], float *match, char **s1, char **s2, int i1, int lgth2 )
  19 {
  20         char *seq2 = s2[0];
  21         int *intptr = mtx[(int)s1[0][i1]];
  22
  23         while( lgth2-- )
  24                 *match++ = intptr[(int)*seq2++];
  25 }
  26 #else
  27 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
  28 {
  29         int j;
  30
  31         for( j=0; j<lgth2; j++ )
  32                 match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
  33 }
  34 #endif
  35
  36 static float Atracking( float *lasthorizontalw, float *lastverticalw,
  37                                                 char **seq1, char **seq2,
  38                         char **mseq1, char **mseq2,
  39                         int **ijp )
  40 {
  41         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
  42         char gap[] = "-";
  43         lgth1 = strlen( seq1[0] );
  44         lgth2 = strlen( seq2[0] );
  45
  46
  47 #if 0
  48         for( i=0; i<lgth1; i++ )
  49         {
  50                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
  51         }
  52 #endif
  53
  54     for( i=0; i<lgth1+1; i++ )
  55     {
  56         ijp[i][0] = i + 1;
  57     }
  58     for( j=0; j<lgth2+1; j++ )
  59     {
  60         ijp[0][j] = -( j + 1 );
  61     }
  62
  63
  64         mseq1[0] += lgth1+lgth2;
  65         *mseq1[0] = 0;
  66         mseq2[0] += lgth1+lgth2;
  67         *mseq2[0] = 0;
  68         iin = lgth1; jin = lgth2;
  69         limk = lgth1+lgth2 + 1;
  70         for( k=0; k<limk; k++ )
  71         {
  72                 if( ijp[iin][jin] < 0 )
  73                 {
  74                         ifi = iin-1; jfi = jin+ijp[iin][jin];
  75                 }
  76                 else if( ijp[iin][jin] > 0 )
  77                 {
  78                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
  79                 }
  80                 else
  81                 {
  82                         ifi = iin-1; jfi = jin-1;
  83                 }
  84                 l = iin - ifi;
  85                 while( --l )
  86                 {
  87                         *--mseq1[0] = seq1[0][ifi+l];
  88                         *--mseq2[0] = *gap;
  89                         k++;
  90                 }
  91                 l= jin - jfi;
  92                 while( --l )
  93                 {
  94                         *--mseq1[0] = *gap;
  95                         *--mseq2[0] = seq2[0][jfi+l];
  96                         k++;
  97                 }
  98                 if( iin <= 0 || jin <= 0 ) break;
  99                 *--mseq1[0] = seq1[0][ifi];
 100                 *--mseq2[0] = seq2[0][jfi];
 101                 k++;
 102                 iin = ifi; jin = jfi;
 103         }
 104         return( 0.0 );
 105 }
 106
 107
 108 float G__align11( char **seq1, char **seq2, int alloclen )
 109 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 110 {
 111 //      int k;
 112         register int i, j;
 113         int lasti;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
 114         int lgth1, lgth2;
 115         int resultlen;
 116         float wm;   /* int ?????? */
 117         float g;
 118         float *currentw, *previousw;
 119         float fpenalty = (float)penalty;
 120 #if USE_PENALTY_EX
 121         float fpenalty_ex = (float)penalty_ex;
 122 #endif
 123 #if 1
 124         float *wtmp;
 125         int *ijppt;
 126         float *mjpt, *prept, *curpt;
 127         int *mpjpt;
 128 #endif
 129         static float mi, *m;
 130         static int **ijp;
 131         static int mpi, *mp;
 132         static float *w1, *w2;
 133         static float *match;
 134         static float *initverticalw;    /* kufuu sureba iranai */
 135         static float *lastverticalw;    /* kufuu sureba iranai */
 136         static char **mseq1;
 137         static char **mseq2;
 138         static char **mseq;
 139         static int **intwork;
 140         static float **floatwork;
 141         static int orlgth1 = 0, orlgth2 = 0;
 142
 143         wm = 0.0;
 144
 145         if( orlgth1 == 0 )
 146         {
 147                 mseq1 = AllocateCharMtx( njob, 0 );
 148                 mseq2 = AllocateCharMtx( njob, 0 );
 149         }
 150
 151
 152         lgth1 = strlen( seq1[0] );
 153         lgth2 = strlen( seq2[0] );
 154
 155
 156
 157         if( lgth1 <= 0 || lgth2 <= 0 )
 158         {
 159                 fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
 160         }
 161
 162         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 163         {
 164                 int ll1, ll2;
 165
 166                 if( orlgth1 > 0 && orlgth2 > 0 )
 167                 {
 168                         FreeFloatVec( w1 );
 169                         FreeFloatVec( w2 );
 170                         FreeFloatVec( match );
 171                         FreeFloatVec( initverticalw );
 172                         FreeFloatVec( lastverticalw );
 173
 174                         FreeFloatVec( m );
 175                         FreeIntVec( mp );
 176
 177                         FreeCharMtx( mseq );
 178
 179
 180
 181                         FreeFloatMtx( floatwork );
 182                         FreeIntMtx( intwork );
 183                 }
 184
 185                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 186                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 187
 188 #if DEBUG
 189                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 190 #endif
 191
 192                 w1 = AllocateFloatVec( ll2+2 );
 193                 w2 = AllocateFloatVec( ll2+2 );
 194                 match = AllocateFloatVec( ll2+2 );
 195
 196                 initverticalw = AllocateFloatVec( ll1+2 );
 197                 lastverticalw = AllocateFloatVec( ll1+2 );
 198
 199                 m = AllocateFloatVec( ll2+2 );
 200                 mp = AllocateIntVec( ll2+2 );
 201
 202                 mseq = AllocateCharMtx( njob, ll1+ll2 );
 203
 204
 205                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
 206                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
 207
 208 #if DEBUG
 209                 fprintf( stderr, "succeeded\n" );
 210 #endif
 211
 212                 orlgth1 = ll1 - 100;
 213                 orlgth2 = ll2 - 100;
 214         }
 215
 216
 217         mseq1[0] = mseq[0];
 218         mseq2[0] = mseq[1];
 219
 220
 221         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
 222         {
 223                 int ll1, ll2;
 224
 225                 if( commonAlloc1 && commonAlloc2 )
 226                 {
 227                         FreeIntMtx( commonIP );
 228                 }
 229
 230                 ll1 = MAX( orlgth1, commonAlloc1 );
 231                 ll2 = MAX( orlgth2, commonAlloc2 );
 232
 233 #if DEBUG
 234                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
 235 #endif
 236
 237                 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
 238
 239 #if DEBUG
 240                 fprintf( stderr, "succeeded\n\n" );
 241 #endif
 242
 243                 commonAlloc1 = ll1;
 244                 commonAlloc2 = ll2;
 245         }
 246         ijp = commonIP;
 247
 248
 249 #if 0
 250         for( i=0; i<lgth1; i++ )
 251                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
 252 #endif
 253
 254         currentw = w1;
 255         previousw = w2;
 256
 257
 258         match_calc( initverticalw, seq2, seq1, 0, lgth1 );
 259
 260
 261         match_calc( currentw, seq1, seq2, 0, lgth2 );
 262
 263         if( outgap == 1 )
 264         {
 265                 for( i=1; i<lgth1+1; i++ )
 266                 {
 267                         initverticalw[i] += fpenalty;
 268                 }
 269                 for( j=1; j<lgth2+1; j++ )
 270                 {
 271                         currentw[j] += fpenalty;
 272                 }
 273         }
 274
 275         for( j=1; j<lgth2+1; ++j )
 276         {
 277                 m[j] = currentw[j-1]; mp[j] = 0;
 278         }
 279
 280         if( lgth2 == 0 )
 281                 lastverticalw[0] = 0.0;               // lgth2==0 no toki error
 282         else
 283                 lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
 284
 285         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
 286
 287 #if XXXXXXX
 288 fprintf( stderr, "currentw = \n" );
 289 for( i=0; i<lgth1+1; i++ )
 290 {
 291         fprintf( stderr, "%5.2f ", currentw[i] );
 292 }
 293 fprintf( stderr, "\n" );
 294 fprintf( stderr, "initverticalw = \n" );
 295 for( i=0; i<lgth2+1; i++ )
 296 {
 297         fprintf( stderr, "%5.2f ", initverticalw[i] );
 298 }
 299 fprintf( stderr, "\n" );
 300 #endif
 301
 302         for( i=1; i<lasti; i++ )
 303         {
 304                 wtmp = previousw;
 305                 previousw = currentw;
 306                 currentw = wtmp;
 307
 308                 previousw[0] = initverticalw[i-1];
 309
 310                 match_calc( currentw, seq1, seq2, i, lgth2 );
 311 #if XXXXXXX
 312 fprintf( stderr, "\n" );
 313 fprintf( stderr, "i=%d\n", i );
 314 fprintf( stderr, "currentw = \n" );
 315 for( j=0; j<lgth2; j++ )
 316 {
 317         fprintf( stderr, "%5.2f ", currentw[j] );
 318 }
 319 fprintf( stderr, "\n" );
 320 #endif
 321 #if XXXXXXX
 322 fprintf( stderr, "\n" );
 323 fprintf( stderr, "i=%d\n", i );
 324 fprintf( stderr, "currentw = \n" );
 325 for( j=0; j<lgth2; j++ )
 326 {
 327         fprintf( stderr, "%5.2f ", currentw[j] );
 328 }
 329 fprintf( stderr, "\n" );
 330 #endif
 331                 currentw[0] = initverticalw[i];
 332
 333                 mi = previousw[0]; mpi = 0;
 334
 335                 ijppt = ijp[i] + 1;
 336                 mjpt = m + 1;
 337                 prept = previousw;
 338                 curpt = currentw + 1;
 339                 mpjpt = mp + 1;
 340                 for( j=1; j<lgth2+1; j++ )
 341                 {
 342                         wm = *prept;
 343                         *ijppt = 0;
 344
 345 #if 0
 346                         fprintf( stderr, "%5.0f->", wm );
 347 #endif
 348 #if 0
 349                         fprintf( stderr, "%5.0f?", g );
 350 #endif
 351                         if( (g=mi+fpenalty) > wm )
 352                         {
 353                                 wm = g;
 354                                 *ijppt = -( j - mpi );
 355                         }
 356                         if( (g=*prept) >= mi )
 357                         {
 358                                 mi = g;
 359                                 mpi = j-1;
 360                         }
 361 #if USE_PENALTY_EX
 362                         mi += fpenalty_ex;
 363 #endif
 364
 365 #if 0
 366                         fprintf( stderr, "%5.0f?", g );
 367 #endif
 368                         if( (g=*mjpt + fpenalty) > wm )
 369                         {
 370                                 wm = g;
 371                                 *ijppt = +( i - *mpjpt );
 372                         }
 373                         if( (g=*prept) >= *mjpt )
 374                         {
 375                                 *mjpt = g;
 376                                 *mpjpt = i-1;
 377                         }
 378 #if USE_PENALTY_EX
 379                         m[j] += fpenalty_ex;
 380 #endif
 381
 382 #if 0
 383                         fprintf( stderr, "%5.0f ", wm );
 384 #endif
 385                         *curpt++ += wm;
 386                         ijppt++;
 387                         mjpt++;
 388                         prept++;
 389                         mpjpt++;
 390                 }
 391                 lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
 392         }
 393
 394         Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp );
 395
 396         resultlen = strlen( mseq1[0] );
 397         if( alloclen < resultlen || resultlen > N )
 398         {
 399                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
 400                 ErrorExit( "LENGTH OVER!\n" );
 401         }
 402
 403
 404         strcpy( seq1[0], mseq1[0] );
 405         strcpy( seq2[0], mseq2[0] );
 406 #if 0
 407         fprintf( stderr, "\n" );
 408         fprintf( stderr, ">\n%s\n", mseq1[0] );
 409         fprintf( stderr, ">\n%s\n", mseq2[0] );
 410         fprintf( stderr, "wm = %f\n", wm );
 411 #endif
 412
 413         return( wm );
 414 }
 415
 416 float G__align11_noalign( int scoremtx[0x80][0x80], int penal, int penal_ex, char **seq1, char **seq2, int alloclen )
 417 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 418 {
 419 //      int k;
 420         register int i, j;
 421         int lasti;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
 422         int lgth1, lgth2;
 423 //      int resultlen;
 424         float wm;   /* int ?????? */
 425         float g;
 426         float *currentw, *previousw;
 427         float fpenalty = (float)penal;
 428 #if USE_PENALTY_EX
 429         float fpenalty_ex = (float)penal_ex;
 430 #endif
 431 #if 1
 432         float *wtmp;
 433         float *mjpt, *prept, *curpt;
 434 //      int *mpjpt;
 435 #endif
 436         static float mi, *m;
 437         static float *w1, *w2;
 438         static float *match;
 439         static float *initverticalw;    /* kufuu sureba iranai */
 440         static float *lastverticalw;    /* kufuu sureba iranai */
 441         static int **intwork;
 442         static float **floatwork;
 443         static int orlgth1 = 0, orlgth2 = 0;
 444
 445         wm = 0.0;
 446
 447
 448
 449         lgth1 = strlen( seq1[0] );
 450         lgth2 = strlen( seq2[0] );
 451
 452
 453
 454         if( lgth1 <= 0 || lgth2 <= 0 )
 455         {
 456                 fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
 457         }
 458
 459         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 460         {
 461                 int ll1, ll2;
 462
 463                 if( orlgth1 > 0 && orlgth2 > 0 )
 464                 {
 465                         FreeFloatVec( w1 );
 466                         FreeFloatVec( w2 );
 467                         FreeFloatVec( match );
 468                         FreeFloatVec( initverticalw );
 469                         FreeFloatVec( lastverticalw );
 470
 471                         FreeFloatVec( m );
 472
 473
 474
 475
 476                         FreeFloatMtx( floatwork );
 477                         FreeIntMtx( intwork );
 478                 }
 479
 480                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 481                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 482
 483 #if DEBUG
 484                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 485 #endif
 486
 487                 w1 = AllocateFloatVec( ll2+2 );
 488                 w2 = AllocateFloatVec( ll2+2 );
 489                 match = AllocateFloatVec( ll2+2 );
 490
 491                 initverticalw = AllocateFloatVec( ll1+2 );
 492                 lastverticalw = AllocateFloatVec( ll1+2 );
 493
 494                 m = AllocateFloatVec( ll2+2 );
 495
 496
 497
 498                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
 499                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
 500
 501 #if DEBUG
 502                 fprintf( stderr, "succeeded\n" );
 503 #endif
 504
 505                 orlgth1 = ll1 - 100;
 506                 orlgth2 = ll2 - 100;
 507         }
 508
 509
 510
 511
 512
 513
 514 #if 0
 515         for( i=0; i<lgth1; i++ )
 516                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
 517 #endif
 518
 519         currentw = w1;
 520         previousw = w2;
 521
 522
 523         match_calc_mtx( scoremtx, initverticalw, seq2, seq1, 0, lgth1 );
 524
 525
 526         match_calc_mtx( scoremtx, currentw, seq1, seq2, 0, lgth2 );
 527
 528         if( outgap == 1 )
 529         {
 530                 for( i=1; i<lgth1+1; i++ )
 531                 {
 532                         initverticalw[i] += fpenalty;
 533                 }
 534                 for( j=1; j<lgth2+1; j++ )
 535                 {
 536                         currentw[j] += fpenalty;
 537                 }
 538         }
 539
 540         for( j=1; j<lgth2+1; ++j )
 541         {
 542                 m[j] = currentw[j-1];
 543         }
 544
 545         if( lgth2 == 0 )
 546                 lastverticalw[0] = 0.0;               // lgth2==0 no toki error
 547         else
 548                 lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
 549
 550         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
 551
 552 #if XXXXXXX
 553 fprintf( stderr, "currentw = \n" );
 554 for( i=0; i<lgth1+1; i++ )
 555 {
 556         fprintf( stderr, "%5.2f ", currentw[i] );
 557 }
 558 fprintf( stderr, "\n" );
 559 fprintf( stderr, "initverticalw = \n" );
 560 for( i=0; i<lgth2+1; i++ )
 561 {
 562         fprintf( stderr, "%5.2f ", initverticalw[i] );
 563 }
 564 fprintf( stderr, "\n" );
 565 #endif
 566
 567         for( i=1; i<lasti; i++ )
 568         {
 569                 wtmp = previousw;
 570                 previousw = currentw;
 571                 currentw = wtmp;
 572
 573                 previousw[0] = initverticalw[i-1];
 574
 575                 match_calc_mtx( scoremtx, currentw, seq1, seq2, i, lgth2 );
 576 #if XXXXXXX
 577 fprintf( stderr, "\n" );
 578 fprintf( stderr, "i=%d\n", i );
 579 fprintf( stderr, "currentw = \n" );
 580 for( j=0; j<lgth2; j++ )
 581 {
 582         fprintf( stderr, "%5.2f ", currentw[j] );
 583 }
 584 fprintf( stderr, "\n" );
 585 #endif
 586 #if XXXXXXX
 587 fprintf( stderr, "\n" );
 588 fprintf( stderr, "i=%d\n", i );
 589 fprintf( stderr, "currentw = \n" );
 590 for( j=0; j<lgth2; j++ )
 591 {
 592         fprintf( stderr, "%5.2f ", currentw[j] );
 593 }
 594 fprintf( stderr, "\n" );
 595 #endif
 596                 currentw[0] = initverticalw[i];
 597
 598                 mi = previousw[0];
 599
 600                 mjpt = m + 1;
 601                 prept = previousw;
 602                 curpt = currentw + 1;
 603                 for( j=1; j<lgth2+1; j++ )
 604                 {
 605                         wm = *prept;
 606
 607 #if 0
 608                         fprintf( stderr, "%5.0f->", wm );
 609 #endif
 610 #if 0
 611                         fprintf( stderr, "%5.0f?", g );
 612 #endif
 613                         if( (g=mi+fpenalty) > wm )
 614                         {
 615                                 wm = g;
 616                         }
 617                         if( (g=*prept) >= mi )
 618                         {
 619                                 mi = g;
 620                         }
 621 #if USE_PENALTY_EX
 622                         mi += fpenalty_ex;
 623 #endif
 624
 625 #if 0
 626                         fprintf( stderr, "%5.0f?", g );
 627 #endif
 628                         if( (g=*mjpt + fpenalty) > wm )
 629                         {
 630                                 wm = g;
 631                         }
 632                         if( (g=*prept) >= *mjpt )
 633                         {
 634                                 *mjpt = g;
 635                         }
 636 #if USE_PENALTY_EX
 637                         m[j] += fpenalty_ex;
 638 #endif
 639
 640 #if 0
 641                         fprintf( stderr, "%5.0f ", wm );
 642 #endif
 643                         *curpt++ += wm;
 644                         mjpt++;
 645                         prept++;
 646                 }
 647                 lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
 648         }
 649
 650 #if 0
 651         fprintf( stderr, "\n" );
 652         fprintf( stderr, ">\n%s\n", mseq1[0] );
 653         fprintf( stderr, ">\n%s\n", mseq2[0] );
 654         fprintf( stderr, "wm = %f\n", wm );
 655 #endif
 656
 657         return( wm );
 658 }