website/archive/binaries/mac/src/clustalw/src/multipleAlign/MyersMillerProfileAlign.cpp

   1 /**
   2  * Author: Mark Larkin
   3  *
   4  * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
   5  * Changes: Mark 20-6-07, I added a call to calculateMaxlengths
   6  */
   7 #ifdef HAVE_CONFIG_H
   8     #include "config.h"
   9 #endif
  10 #include "MyersMillerProfileAlign.h"
  11 #include "Iteration.h"
  12 #include <math.h>
  13
  14 namespace clustalw
  15 {
  16
  17 /**
  18  *
  19  * @return
  20  */
  21 MyersMillerProfileAlign::MyersMillerProfileAlign()
  22  : _gapPos1(userParameters->getGapPos1()),
  23    _gapPos2(userParameters->getGapPos2())
  24 {
  25
  26 }
  27
  28 /**
  29  *
  30  * @param alnPtr
  31  * @param distMat
  32  * @param group
  33  * @param aligned
  34  * @return
  35  */
  36 int MyersMillerProfileAlign::profileAlign(Alignment* alnPtr, DistMatrix* distMat,
  37                                           vector<int>* group, int* aligned)
  38 {
  39     bool negative = false;
  40     int i = 0, j = 0, count = 0;
  41     int numSeq = 0;
  42     int seqNum = 0;
  43     int len = 0, len1 = 0, len2 = 0, is = 0, minLen = 0;
  44     int se1 = 0, se2 = 0, sb1 = 0, sb2 = 0;
  45     int maxRes = 0;
  46     int c = 0;
  47     int score = 0;
  48     double logmin = 0.0, logdiff = 0.0;
  49     double pcid = 0.0;
  50     int matrix[NUMRES][NUMRES];
  51     int numSeqsProf1 = 0, numSeqsProf2 = 0;
  52     // Initialise the matrix. To get rid of a valgrind error.
  53     for(int i = 0; i < NUMRES; i++)
  54     {
  55         for(int j = 0; j < NUMRES; j++)
  56         {
  57             matrix[i][j] = 0;
  58         }
  59     }
  60
  61     numSeq = alnPtr->getNumSeqs();
  62     seqArray.resize(numSeq);
  63     alnWeight.resize(numSeq);
  64
  65     for (i = 0; i < numSeq; i++)
  66     {
  67         if (aligned[i + 1] == 0)
  68         {
  69             (*group)[i + 1] = 0;
  70         }
  71     }
  72
  73     nseqs1 = nseqs2 = 0;
  74     for (i = 0; i < numSeq; i++)
  75     {
  76         if ((*group)[i + 1] == 1)
  77         {
  78             nseqs1++;
  79         }
  80         else if ((*group)[i + 1] == 2)
  81         {
  82             nseqs2++;
  83         }
  84     }
  85
  86     if ((nseqs1 == 0) || (nseqs2 == 0))
  87     {
  88         return (0);
  89     }
  90     numSeqsProf1 = nseqs1;
  91     numSeqsProf2 = nseqs2;
  92
  93     if (nseqs2 > nseqs1)
  94     {
  95         switchProfiles = true;
  96         for (i = 0; i < numSeq; i++)
  97         {
  98             if ((*group)[i + 1] == 1)
  99             {
 100                 (*group)[i + 1] = 2;
 101             }
 102             else if ((*group)[i + 1] == 2)
 103             {
 104                 (*group)[i + 1] = 1;
 105             }
 106         }
 107     }
 108     else
 109     {
 110         switchProfiles = false;
 111     }
 112
 113     // calculate the mean of the sequence pc identities between the two groups
 114
 115     count = 0;
 116     pcid = 0.0;
 117     negative = userParameters->getUseNegMatrix();
 118
 119     for (i = 0; i < numSeq; i++)
 120     {
 121         if ((*group)[i + 1] == 1)
 122         {
 123             for (j = 0; j < numSeq; j++)
 124             {
 125                 if ((*group)[j + 1] == 2)
 126                 {
 127                     count++;
 128                     pcid += (*distMat)(i + 1, j + 1);
 129                 }
 130             }
 131         }
 132     }
 133
 134     pcid = pcid / (float)count;
 135
 136     // Make the first profile.
 137
 138     prfLength1 = 0;
 139     for (i = 0; i < numSeq; i++)
 140     {
 141         if ((*group)[i + 1] == 1)
 142         {
 143             if (alnPtr->getSeqLength(i + 1) > prfLength1)
 144             {
 145                 prfLength1 = alnPtr->getSeqLength(i + 1);
 146             }
 147         }
 148     }
 149     nseqs1 = 0;
 150
 151     for (i = 0; i < numSeq; i++)
 152     {
 153         if ((*group)[i + 1] == 1)
 154         {
 155             len = alnPtr->getSeqLength(i + 1);
 156
 157             // initialise with the other vector!
 158             seqArray[nseqs1] = vector<int>(alnPtr->getSequence(i + 1)->begin() + 1,
 159                                             alnPtr->getSequence(i + 1)->end());
 160             seqArray[nseqs1].resize(prfLength1 + extraEndElemNum);
 161
 162             for (j = len; j < prfLength1; j++)
 163             {
 164                 seqArray[nseqs1][j] = userParameters->getGapPos1();
 165             }
 166             alnWeight[nseqs1] = alnPtr->getSeqWeight(i);
 167             nseqs1++;
 168         }
 169     }
 170
 171
 172     // Make the second profile.
 173
 174     prfLength2 = 0;
 175     for (i = 0; i < numSeq; i++)
 176     {
 177         if ((*group)[i + 1] == 2)
 178         {
 179             if (alnPtr->getSeqLength(i + 1) > prfLength2)
 180             {
 181                 prfLength2 = alnPtr->getSeqLength(i + 1);
 182             }
 183         }
 184     }
 185     nseqs2 = 0;
 186
 187     for (i = 0; i < numSeq; i++)
 188     {
 189         if ((*group)[i + 1] == 2)
 190         {
 191             len = alnPtr->getSeqLength(i + 1);
 192
 193             seqArray[nseqs1 + nseqs2] = vector<int>(alnPtr->getSequence(i + 1)->begin() + 1,
 194                                                      alnPtr->getSequence(i + 1)->end());
 195             seqArray[nseqs1 + nseqs2].resize(prfLength2 + extraEndElemNum);
 196
 197             for (j = len; j < prfLength2; j++)
 198             {
 199                 seqArray[nseqs1 + nseqs2][j] = userParameters->getGapPos1();
 200             }
 201
 202             seqArray[nseqs1 + nseqs2][j] = ENDALN;
 203             alnWeight[nseqs1 + nseqs2] = alnPtr->getSeqWeight(i);
 204             //cout << "weight " << nseqs1 + nseqs2 << alnWeight[nseqs1 + nseqs2] << "\n";
 205             nseqs2++;
 206         }
 207     }
 208
 209     // Change the Max alignment length in the Alignment Object!
 210     alnPtr->setMaxAlnLength(prfLength1 + prfLength2 + 2);
 211
 212     // calculate real length of profiles - removing gaps!
 213
 214     len1 = 0;
 215     for (i = 0; i < nseqs1; i++)
 216     {
 217         is = 0;
 218         for (j = 0; j < utilityObject->MIN(static_cast<int>(seqArray[i].size()), prfLength1); j++)
 219         {
 220             c = seqArray[i][j];
 221             if ((c != _gapPos1) && (c != _gapPos2))
 222             {
 223                 is++;
 224             }
 225         }
 226         len1 += is;
 227     }
 228     len1 = static_cast<int>(len1 / (float)nseqs1);
 229
 230     len2 = 0;
 231     for (i = nseqs1; i < nseqs2 + nseqs1; i++)
 232     {
 233         is = 0;
 234         for (j = 0; j < utilityObject->MIN(static_cast<int>(seqArray[i].size()), prfLength2);
 235              j++)
 236         {
 237             c = seqArray[i][j];
 238             if ((c != _gapPos1) && (c != _gapPos2))
 239             {
 240                 is++;
 241             }
 242         }
 243         len2 += is;
 244     }
 245     len2 = static_cast<int>(len2 / (float)nseqs2);
 246
 247     PrfScaleValues scaleVals;
 248     scaleVals.scale = 1.0;
 249     scaleVals.intScale = 100.0;
 250
 251     int matAvgScore = 0;
 252     minLen = utilityObject->MIN(len1, len2);
 253     maxRes = 0;
 254
 255     // Get the substitution matrix that will be stored in 'matrix'
 256     maxRes = subMatrix->getProfileAlignMatrix(matrix, pcid, minLen, scaleVals, matAvgScore);
 257
 258     if (maxRes == 0 || maxRes == -1)
 259     {
 260         return -1;
 261     }
 262     if (userParameters->getDNAFlag())
 263     {
 264         gapcoef1 = gapcoef2 = static_cast<int>(100.0 * userParameters->getGapOpen() * scaleVals.scale);
 265         lencoef1 = lencoef2 = static_cast<int>(100.0 * userParameters->getGapExtend() * scaleVals.scale);
 266     }
 267     else
 268     {
 269         if (len1 == 0 || len2 == 0)
 270         {
 271             logmin = 1.0;
 272             logdiff = 1.0;
 273         }
 274         else
 275         {
 276             logmin = 1.0 / log10((double)minLen);
 277             if (len2 < len1)
 278             {
 279                 logdiff = 1.0 + 0.5 * log10((double)((float)len2 / (float)len1));
 280             }
 281             else if (len1 < len2)
 282             {
 283                 logdiff = 1.0 + 0.5 * log10((double)((float)len1 / (float)len2));
 284             }
 285             else
 286             {
 287                 logdiff = 1.0;
 288             }
 289             if (logdiff < 0.9)
 290             {
 291                 logdiff = 0.9;
 292             }
 293         }
 294
 295
 296         if (negative)
 297         {
 298             gapcoef1 = gapcoef2 = static_cast<int>(100.0 *(float)(userParameters->getGapOpen()));
 299             lencoef1 = lencoef2 = static_cast<int>(100.0 * userParameters->getGapExtend());
 300         }
 301         else
 302         {
 303             if (matAvgScore <= 0)
 304             {
 305                 gapcoef1 = gapcoef2 = static_cast<int>(100.0 *(float)(userParameters->getGapOpen() + logmin));
 306             }
 307             else
 308             {
 309                 gapcoef1 = gapcoef2 = static_cast<int>(scaleVals.scale * matAvgScore *
 310                             (float)(userParameters->getGapOpen() / (logdiff * logmin)));
 311             }
 312             lencoef1 = lencoef2 = static_cast<int>(100.0 * userParameters->getGapExtend());
 313         }
 314     }
 315     // We need one profile with substitution matrix information and one without!
 316     // But this will change when we have the LE scoring function.
 317     profileWithSub = new ProfileWithSub(prfLength1, 0, nseqs1);
 318     profileStandard = new ProfileStandard(prfLength2, nseqs1, nseqs1 + nseqs2);
 319
 320     // Calculate the profile array with Substitution matrix info
 321     // calculate the Gap Coefficients.
 322
 323     gaps.resize(alnPtr->getMaxAlnLength() + 1);
 324
 325     bool profile1Pen;
 326     if (switchProfiles == false)
 327     {
 328         profile1Pen = userParameters->getStructPenalties1() && userParameters->getUseSS1();
 329         profileWithSub->calcGapCoeff(&seqArray, &gaps, profile1Pen,
 330                                      alnPtr->getGapPenaltyMask1(), gapcoef1, lencoef1);
 331     }
 332     else
 333     {
 334         profile1Pen = userParameters->getStructPenalties2() && userParameters->getUseSS2();
 335         profileWithSub->calcGapCoeff(&seqArray, &gaps, profile1Pen,
 336                                      alnPtr->getGapPenaltyMask2(), gapcoef1, lencoef1);
 337     }
 338     // calculate the profile matrix.
 339     profileWithSub->calcProfileWithSub(&seqArray, &gaps, matrix, &alnWeight);
 340     profile1 = profileWithSub->getProfilePtr();
 341
 342     if (userParameters->getDebug() > 4)
 343     {
 344         string aminoAcidCodes = userParameters->getAminoAcidCodes();
 345         for (j = 0; j <= userParameters->getMaxAA(); j++)
 346         {
 347             cout << aminoAcidCodes[j] << "    ";
 348         }
 349         cout << "\n";
 350         for (i = 0; i < prfLength1; i++)
 351         {
 352             for (j = 0; j <= userParameters->getMaxAA(); j++)
 353             {
 354                 cout << (*profile1)[i + 1][j] << " ";
 355             }
 356             cout << (*profile1)[i + 1][_gapPos1]<< " ";
 357             cout << (*profile1)[i + 1][_gapPos2]<< " ";
 358             cout << (*profile1)[i + 1][GAPCOL] << " " << (*profile1)[i + 1][LENCOL]
 359                  << "\n";
 360         }
 361     }
 362     // Calculate the standard profile array
 363     // calculate the Gap Coefficients.
 364     bool profile2Pen;
 365     if (switchProfiles == false)
 366     {
 367         profile2Pen = userParameters->getStructPenalties2() && userParameters->getUseSS2();
 368         profileStandard->calcGapCoeff(&seqArray, &gaps, profile2Pen,
 369                                       alnPtr->getGapPenaltyMask2(), gapcoef2, lencoef2);
 370     }
 371     else
 372     {
 373         profile2Pen = userParameters->getStructPenalties1() && userParameters->getUseSS1();
 374         profileStandard->calcGapCoeff(&seqArray, &gaps, profile2Pen,
 375                                       alnPtr->getGapPenaltyMask1(), gapcoef2, lencoef2);
 376     }
 377
 378     // calculate the profile matrix.
 379
 380     profileStandard->calcStandardProfile(&seqArray, &alnWeight);
 381     profile2 = profileStandard->getProfilePtr();
 382
 383     if (userParameters->getDebug() > 4)
 384     {
 385         string aminoAcidCodes = userParameters->getAminoAcidCodes();
 386         for (j = 0; j <= userParameters->getMaxAA(); j++)
 387         {
 388             cout << aminoAcidCodes[j] << "    ";
 389         }
 390         cout << "\n";
 391         for (i = 0; i < prfLength2; i++)
 392         {
 393             for (j = 0; j <= userParameters->getMaxAA(); j++)
 394             {
 395                 cout << (*profile2)[i + 1][j] << " ";
 396             }
 397             cout << (*profile2)[i + 1][_gapPos1]<< " ";
 398             cout << (*profile2)[i + 1][_gapPos2]<< " ";
 399             cout << (*profile2)[i + 1][GAPCOL] << " " << (*profile2)[i + 1][LENCOL]
 400                  << "\n";
 401         }
 402     }
 403     alnWeight.clear();
 404
 405     int _maxAlnLength = alnPtr->getMaxAlnLength();
 406
 407     alnPath1.resize(_maxAlnLength + 1);
 408     alnPath2.resize(_maxAlnLength + 1);
 409
 410     //align the profiles
 411
 412     // use Myers and Miller to align two sequences
 413
 414     lastPrint = 0;
 415     printPtr = 1;
 416
 417     sb1 = sb2 = 0;
 418     se1 = prfLength1;
 419     se2 = prfLength2;
 420
 421     HH.resize(_maxAlnLength + 1);
 422     DD.resize(_maxAlnLength + 1);
 423     RR.resize(_maxAlnLength + 1);
 424     SS.resize(_maxAlnLength + 1);
 425     gS.resize(_maxAlnLength + 1);
 426     displ.resize(_maxAlnLength + 1);
 427
 428     score = progDiff(sb1, sb2, se1 - sb1, se2 - sb2, (*profile1)[0][GAPCOL],
 429                     (*profile1)[prfLength1][GAPCOL]);
 430
 431     // May not need if we recreate the Profile every time!
 432     HH.clear();
 433     DD.clear();
 434     RR.clear();
 435     SS.clear();
 436     gS.clear();
 437
 438     alignmentLength = progTracepath();
 439
 440     displ.clear();
 441
 442     addGGaps(alnPtr, &seqArray);
 443
 444     profileStandard->resetProfile();
 445     profileWithSub->resetProfile();
 446
 447     prfLength1 = alignmentLength;
 448
 449     alnPath1.clear();
 450     alnPath2.clear();
 451
 452     if(userParameters->getDoRemoveFirstIteration() == TREE)
 453     {
 454         Iteration iterateObj;
 455         iterateObj.iterationOnTreeNode(numSeqsProf1, numSeqsProf2, prfLength1, prfLength2,
 456                                        &seqArray);
 457     }
 458
 459     //  Now we resize the SeqArray that holds the sequences in the alignment class
 460     //    and also update it with the new aligned sequences
 461
 462     SeqArray* newSequences = alnPtr->getSeqArrayForRealloc();
 463     seqNum = 0;
 464     for (j = 0; j < numSeq; j++)
 465     {
 466         if ((*group)[j + 1] == 1)
 467         {
 468             (*newSequences)[j + 1].clear();
 469             (*newSequences)[j + 1].resize(prfLength1 + 1);
 470             for (i = 0; i < prfLength1; i++)
 471             {
 472                 (*newSequences)[j + 1][i + 1] = seqArray[seqNum][i];
 473             }
 474             seqNum++;
 475         }
 476     }
 477     for (j = 0; j < numSeq; j++)
 478     {
 479         if ((*group)[j + 1] == 2)
 480         {
 481             (*newSequences)[j + 1].clear();
 482             (*newSequences)[j + 1].resize(prfLength1 + 1);
 483             for (i = 0; i < prfLength1; i++)
 484             {
 485                 (*newSequences)[j + 1][i + 1] = seqArray[seqNum][i];
 486             }
 487             seqNum++;
 488         }
 489     }
 490
 491     alnPtr->calculateMaxLengths(); // Mark change 20-6-07
 492
 493     gaps.clear();
 494     seqArray.clear();
 495
 496     delete profileWithSub;
 497     delete profileStandard;
 498
 499     int retScore = (score / 100);
 500
 501     return retScore;
 502
 503 }
 504
 505 /** ****************************************************************************************
 506  *                          Private functions                                              *
 507  *******************************************************************************************/
 508
 509 /**
 510  *
 511  * @param A
 512  * @param B
 513  * @param M
 514  * @param N
 515  * @param go1
 516  * @param go2
 517  * @return
 518  */
 519 int MyersMillerProfileAlign::progDiff(int A, int B, int M, int N, int go1, int go2)
 520 {
 521     int midi, midj, type;
 522     int midh;
 523
 524     static int t, tl, g, h;
 525
 526     {
 527         static int i, j;
 528         static int hh, f, e, s;
 529
 530         /* Boundary cases: M <= 1 or N == 0 */
 531         if (userParameters->getDebug() > 2)
 532         {
 533             cout << "A " << A << " B " << B << " M " << M << " N " << N
 534                  << " midi " << M / 2 << " go1 " << go1 << " go2 " << go2<< "\n";
 535         }
 536
 537         /* if sequence B is empty....                                            */
 538
 539         if (N <= 0)
 540         {
 541
 542             /* if sequence A is not empty....                                        */
 543
 544             if (M > 0)
 545             {
 546
 547                 /* delete residues A[1] to A[M]                                          */
 548
 549                 progDel(M);
 550             }
 551             return ( - gapPenalty1(A, B, M));
 552         }
 553
 554         /* if sequence A is empty....                                            */
 555
 556         if (M <= 1)
 557         {
 558             if (M <= 0)
 559             {
 560
 561                 /* insert residues B[1] to B[N]                                          */
 562
 563                 progAdd(N);
 564                 return ( - gapPenalty2(A, B, N));
 565             }
 566
 567             /* if sequence A has just one residue....                                */
 568
 569             if (go1 == 0)
 570             {
 571                 midh =  - gapPenalty1(A + 1, B + 1, N);
 572             }
 573             else
 574             {
 575                 midh =  - gapPenalty2(A + 1, B, 1) - gapPenalty1(A + 1, B + 1,
 576                     N);
 577             }
 578             midj = 0;
 579             for (j = 1; j <= N; j++)
 580             {
 581                 hh =  - gapPenalty1(A, B + 1, j - 1) + prfScore(A + 1, B + j)
 582                     - gapPenalty1(A + 1, B + j + 1, N - j);
 583                 if (hh > midh)
 584                 {
 585                     midh = hh;
 586                     midj = j;
 587                 }
 588             }
 589
 590             if (midj == 0)
 591             {
 592                 progAdd(N);
 593                 progDel(1);
 594             }
 595             else
 596             {
 597                 if (midj > 1)
 598                 {
 599                     progAdd(midj - 1);
 600                 }
 601                 progAlign();
 602                 if (midj < N)
 603                 {
 604                     progAdd(N - midj);
 605                 }
 606             }
 607             return midh;
 608         }
 609
 610
 611         /* Divide sequence A in half: midi */
 612
 613         midi = M / 2;
 614
 615         /* In a forward phase, calculate all HH[j] and HH[j] */
 616
 617         HH[0] = 0;
 618         t =  - openPenalty1(A, B + 1);
 619         tl =  - extPenalty1(A, B + 1);
 620         for (j = 1; j <= N; j++)
 621         {
 622             HH[j] = t = t + tl;
 623             DD[j] = t - openPenalty2(A + 1, B + j);
 624         }
 625
 626         if (go1 == 0)
 627         {
 628             t = 0;
 629         }
 630         else
 631         {
 632             t =  - openPenalty2(A + 1, B);
 633         }
 634         tl =  - extPenalty2(A + 1, B);
 635         for (i = 1; i <= midi; i++)
 636         {
 637             s = HH[0];
 638             HH[0] = hh = t = t + tl;
 639             f = t - openPenalty1(A + i, B + 1);
 640
 641             for (j = 1; j <= N; j++)
 642             {
 643                 g = openPenalty1(A + i, B + j);
 644                 h = extPenalty1(A + i, B + j);
 645                 if ((hh = hh - g - h) > (f = f - h))
 646                 {
 647                     f = hh;
 648                 }
 649                 g = openPenalty2(A + i, B + j);
 650                 h = extPenalty2(A + i, B + j);
 651                 if ((hh = HH[j] - g - h) > (e = DD[j] - h))
 652                 {
 653                     e = hh;
 654                 }
 655                 hh = s + prfScore(A + i, B + j);
 656                 if (f > hh)
 657                 {
 658                     hh = f;
 659                 }
 660                 if (e > hh)
 661                 {
 662                     hh = e;
 663                 }
 664
 665                 s = HH[j];
 666                 HH[j] = hh;
 667                 DD[j] = e;
 668
 669             }
 670         }
 671
 672         DD[0] = HH[0];
 673
 674         /* In a reverse phase, calculate all RR[j] and SS[j] */
 675
 676         RR[N] = 0;
 677         tl = 0;
 678         for (j = N - 1; j >= 0; j--)
 679         {
 680             g =  - openPenalty1(A + M, B + j + 1);
 681             tl -= extPenalty1(A + M, B + j + 1);
 682             RR[j] = g + tl;
 683             SS[j] = RR[j] - openPenalty2(A + M, B + j);
 684             gS[j] = openPenalty2(A + M, B + j);
 685         }
 686
 687         tl = 0;
 688         for (i = M - 1; i >= midi; i--)
 689         {
 690             s = RR[N];
 691             if (go2 == 0)
 692             {
 693                 g = 0;
 694             }
 695             else
 696             {
 697                 g =  - openPenalty2(A + i + 1, B + N);
 698             }
 699             tl -= extPenalty2(A + i + 1, B + N);
 700             RR[N] = hh = g + tl;
 701             t = openPenalty1(A + i, B + N);
 702             f = RR[N] - t;
 703
 704             for (j = N - 1; j >= 0; j--)
 705             {
 706                 g = openPenalty1(A + i, B + j + 1);
 707                 h = extPenalty1(A + i, B + j + 1);
 708                 if ((hh = hh - g - h) > (f = f - h - g + t))
 709                 {
 710                     f = hh;
 711                 }
 712                 t = g;
 713                 g = openPenalty2(A + i + 1, B + j);
 714                 h = extPenalty2(A + i + 1, B + j);
 715                 hh = RR[j] - g - h;
 716                 if (i == (M - 1))
 717                 {
 718                     e = SS[j] - h;
 719                 }
 720                 else
 721                 {
 722                     e = SS[j] - h - g + openPenalty2(A + i + 2, B + j);
 723                     gS[j] = g;
 724                 }
 725                 if (hh > e)
 726                 {
 727                     e = hh;
 728                 }
 729                 hh = s + prfScore(A + i + 1, B + j + 1);
 730                 if (f > hh)
 731                 {
 732                     hh = f;
 733                 }
 734                 if (e > hh)
 735                 {
 736                     hh = e;
 737                 }
 738
 739                 s = RR[j];
 740                 RR[j] = hh;
 741                 SS[j] = e;
 742
 743             }
 744         }
 745         SS[N] = RR[N];
 746         gS[N] = openPenalty2(A + midi + 1, B + N);
 747
 748         /* find midj, such that HH[j]+RR[j] or DD[j]+SS[j]+gap is the maximum */
 749
 750         midh = HH[0] + RR[0];
 751         midj = 0;
 752         type = 1;
 753         for (j = 0; j <= N; j++)
 754         {
 755             hh = HH[j] + RR[j];
 756             if (hh >= midh)
 757             if (hh > midh || (HH[j] != DD[j] && RR[j] == SS[j]))
 758             {
 759                 midh = hh;
 760                 midj = j;
 761             }
 762         }
 763
 764         for (j = N; j >= 0; j--)
 765         {
 766             hh = DD[j] + SS[j] + gS[j];
 767             if (hh > midh)
 768             {
 769                 midh = hh;
 770                 midj = j;
 771                 type = 2;
 772             }
 773         }
 774     }
 775
 776     /* Conquer recursively around midpoint                                   */
 777
 778
 779     if (type == 1)
 780     {
 781          /* Type 1 gaps  */
 782         if (userParameters->getDebug() > 2)
 783         {
 784             cout << "Type 1,1: midj " << midj << "\n";
 785         }
 786         progDiff(A, B, midi, midj, go1, 1);
 787         if (userParameters->getDebug() > 2)
 788         {
 789             cout << "Type 1,2: midj " << midj << "\n";
 790         }
 791         progDiff(A + midi, B + midj, M - midi, N - midj, 1, go2);
 792     }
 793     else
 794     {
 795         if (userParameters->getDebug() > 2)
 796         {
 797             cout << "Type 2,1: midj " << midj << "\n";
 798         }
 799         progDiff(A, B, midi - 1, midj, go1, 0);
 800         progDel(2);
 801         if (userParameters->getDebug() > 2)
 802         {
 803             cout << "Type 2,2: midj " << midj << "\n";
 804         }
 805         progDiff(A + midi + 1, B + midj, M - midi - 1, N - midj, 0, go2);
 806     }
 807
 808     return midh; /* Return the score of the best alignment */
 809 }
 810
 811 /**
 812  *
 813  * @param alnPtr
 814  * @param seqArray
 815  */
 816 void MyersMillerProfileAlign::addGGaps(Alignment* alnPtr, SeqArray* seqArray)
 817 {
 818     int j;
 819     int i, ix;
 820     int len;
 821     vector<int> ta;
 822
 823     ta.resize(alignmentLength + 1);
 824
 825     for (j = 0; j < nseqs1; j++)
 826     {
 827         ix = 0;
 828         for (i = 0; i < alignmentLength; i++)
 829         {
 830             if (alnPath1[i] == 2)
 831             {
 832                 if (ix < ((int)(*seqArray)[j].size() - extraEndElemNum))
 833                 {
 834                     ta[i] = (*seqArray)[j][ix];
 835                 }
 836                 else
 837                 {
 838                     ta[i] = ENDALN;
 839                 }
 840                 ix++;
 841             }
 842             else if (alnPath1[i] == 1)
 843             {
 844                 /*
 845                 insertion in first alignment...
 846                  */
 847                 ta[i] = _gapPos1;
 848             }
 849             else
 850             {
 851                 cerr << "Error in aln_path\n";
 852             }
 853         }
 854         ta[i] = ENDALN;
 855
 856         len = alignmentLength;
 857
 858         (*seqArray)[j].resize(len + 2);
 859
 860         for (i = 0; i < len; i++)
 861         {
 862             (*seqArray)[j][i] = ta[i];
 863         }
 864         (*seqArray)[j][len] = ENDALN;
 865
 866     }
 867
 868     for (j = nseqs1; j < nseqs1 + nseqs2; j++)
 869     {
 870         ix = 0;
 871         for (i = 0; i < alignmentLength; i++)
 872         {
 873             if (alnPath2[i] == 2)
 874             {
 875                 if (ix < ((int)(*seqArray)[j].size() - extraEndElemNum))
 876                 {
 877                     ta[i] = (*seqArray)[j][ix];
 878                 }
 879                 else
 880                 {
 881                     ta[i] = ENDALN;
 882                 }
 883                 ix++;
 884             }
 885             else if (alnPath2[i] == 1)
 886             {
 887                 /*
 888                 insertion in second alignment...
 889                  */
 890                 ta[i] = _gapPos1;
 891             }
 892             else
 893             {
 894                 cerr << "Error in alnPath\n";
 895             }
 896         }
 897         ta[i] = ENDALN;
 898
 899         len = alignmentLength;
 900
 901         (*seqArray)[j].resize(len + 2);
 902
 903         for (i = 0; i < len; i++)
 904         {
 905             (*seqArray)[j][i] = ta[i];
 906         }
 907         (*seqArray)[j][len] = ENDALN;
 908     }
 909
 910
 911     if (userParameters->getStructPenalties1() != NONE)
 912     {
 913         addGGapsMask(alnPtr->getGapPenaltyMask1(), alignmentLength,
 914                      &alnPath1, &alnPath2);
 915     }
 916     if (userParameters->getStructPenalties1() == SECST)
 917     {
 918         addGGapsMask(alnPtr->getSecStructMask1(), alignmentLength,
 919                      &alnPath1, &alnPath2);
 920     }
 921
 922     if (userParameters->getStructPenalties2() != NONE)
 923     {
 924         addGGapsMask(alnPtr->getGapPenaltyMask2(), alignmentLength,
 925                      &alnPath2, &alnPath1);
 926     }
 927     if (userParameters->getStructPenalties2() == SECST)
 928     {
 929         addGGapsMask(alnPtr->getSecStructMask2(), alignmentLength,
 930                      &alnPath2, &alnPath1);
 931     }
 932 }
 933
 934
 935 /**
 936  *
 937  * @param mask
 938  * @param len
 939  * @param path1
 940  * @param path2
 941  */
 942 void MyersMillerProfileAlign::addGGapsMask(vector<char>* mask, int len, vector<int>* path1,
 943                                            vector<int>* path2)
 944 {
 945     int i, ix;
 946     char *ta;
 947
 948     ta = new char[len + 1];
 949
 950     ix = 0;
 951     if (switchProfiles == false)
 952     {
 953         for (i = 0; i < len; i++)
 954         {
 955             if ((*path1)[i] == 2)
 956             {
 957                 ta[i] = (*mask)[ix];
 958                 ix++;
 959             }
 960             else if ((*path1)[i] == 1)
 961             {
 962                 ta[i] = _gapPos1;
 963             }
 964         }
 965     }
 966     else
 967     {
 968         for (i = 0; i < len; i++)
 969         {
 970             if ((*path2)[i] == 2)
 971             {
 972                 ta[i] = (*mask)[ix];
 973                 ix++;
 974             }
 975             else if ((*path2)[i] == 1)
 976             {
 977                 ta[i] = _gapPos1;
 978             }
 979         }
 980     }
 981     mask->resize(len + 2);
 982
 983     for (i = 0; i < len; i++)
 984     {
 985         (*mask)[i] = ta[i];
 986     }
 987
 988
 989     delete [] ta;
 990     ta  = NULL;
 991 }
 992
 993
 994 /**
 995  *
 996  * @param n
 997  * @param m
 998  * @return
 999  */
1000 inline int MyersMillerProfileAlign::prfScore(int n, int m)
1001 {
1002     int ix;
1003     int score;
1004
1005     score = 0;
1006     int _maxAA = userParameters->getMaxAA(); // NOTE Change here!
1007     for (ix = 0; ix <= _maxAA; ix++)
1008     {
1009         score += ((*profile1)[n][ix] * (*profile2)[m][ix]);
1010     }
1011     score += ((*profile1)[n][_gapPos1] * (*profile2)[m][_gapPos1]);
1012     score += ((*profile1)[n][_gapPos2] * (*profile2)[m][_gapPos2]);
1013     return (score / 10);
1014 }
1015
1016
1017 /**
1018  *
1019  * @return
1020  */
1021 int MyersMillerProfileAlign::progTracepath()
1022 {
1023     int i, j, k, pos, toDo;
1024     int alignLen;
1025     pos = 0;
1026
1027     toDo = printPtr - 1;
1028
1029     for (i = 1; i <= toDo; ++i)
1030     {
1031         if (userParameters->getDebug() > 1)
1032         {
1033             cout << displ[i] << " ";
1034         }
1035         if (displ[i] == 0)
1036         {
1037             alnPath1[pos] = 2;
1038             alnPath2[pos] = 2;
1039             ++pos;
1040         }
1041         else
1042         {
1043             if ((k = displ[i]) > 0)
1044             {
1045                 for (j = 0; j <= k - 1; ++j)
1046                 {
1047                     alnPath2[pos + j] = 2;
1048                     alnPath1[pos + j] = 1;
1049                 }
1050                 pos += k;
1051             }
1052             else
1053             {
1054                 k = (displ[i] < 0) ? displ[i] * - 1: displ[i];
1055                 for (j = 0; j <= k - 1; ++j)
1056                 {
1057                     alnPath1[pos + j] = 2;
1058                     alnPath2[pos + j] = 1;
1059                 }
1060                 pos += k;
1061             }
1062         }
1063     }
1064     if (userParameters->getDebug() > 1)
1065     {
1066         cout << "\n";
1067     }
1068
1069     alignLen = pos;
1070     return alignLen;
1071 }
1072
1073
1074 /**
1075  *
1076  * @param k
1077  */
1078 void MyersMillerProfileAlign::progDel(int k)
1079 {
1080     if (lastPrint < 0)
1081     {
1082         lastPrint = displ[printPtr - 1] -= k;
1083     }
1084     else
1085     {
1086         lastPrint = displ[printPtr++] =  - (k);
1087     }
1088 }
1089
1090
1091 /**
1092  *
1093  * @param k
1094  */
1095 void MyersMillerProfileAlign::progAdd(int k)
1096 {
1097     if (lastPrint < 0)
1098     {
1099         displ[printPtr - 1] = k;
1100         displ[printPtr++] = lastPrint;
1101     }
1102     else
1103     {
1104         lastPrint = displ[printPtr++] = k;
1105     }
1106 }
1107
1108
1109 /**
1110  *
1111  */
1112 void MyersMillerProfileAlign::progAlign()
1113 {
1114     displ[printPtr++] = lastPrint = 0;
1115 }
1116
1117
1118 /**
1119  *
1120  * @param i
1121  * @param j
1122  * @return
1123  */
1124 inline int MyersMillerProfileAlign::openPenalty1(int i, int j) // NOTE Change here!
1125 {
1126     int g;
1127     if (!userParameters->getEndGapPenalties() && (i == 0 || i == prfLength1))
1128     {
1129         return (0);
1130     }
1131
1132     g = (*profile2)[j][GAPCOL] + (*profile1)[i][GAPCOL];
1133     return (g);
1134 }
1135
1136
1137 /**
1138  *
1139  * @param i
1140  * @param j
1141  * @return
1142  */
1143 inline int MyersMillerProfileAlign::extPenalty1(int i, int j) // NOTE Change here!
1144 {
1145     int h;
1146
1147     if (!userParameters->getEndGapPenalties() && (i == 0 || i == prfLength1))
1148     {
1149         return (0);
1150     }
1151
1152     h = (*profile2)[j][LENCOL];
1153     return (h);
1154 }
1155
1156
1157 /**
1158  *
1159  * @param i
1160  * @param j
1161  * @param k
1162  * @return
1163  */
1164 int MyersMillerProfileAlign::gapPenalty1(int i, int j, int k)
1165 {
1166     int ix;
1167     int gp;
1168     int g, h = 0;
1169
1170     if (k <= 0)
1171     {
1172         return (0);
1173     }
1174     if (!userParameters->getEndGapPenalties() && (i == 0 || i == prfLength1))
1175     {
1176         return (0);
1177     }
1178
1179     g = (*profile2)[j][GAPCOL] + (*profile1)[i][GAPCOL];
1180     for (ix = 0; ix < k && ix + j < prfLength2; ix++)
1181     {
1182         h += (*profile2)[ix + j][LENCOL];
1183     }
1184
1185     gp = g + h;
1186     return (gp);
1187 }
1188
1189
1190 /**
1191  *
1192  * @param i
1193  * @param j
1194  * @return
1195  */
1196 inline int MyersMillerProfileAlign::openPenalty2(int i, int j) // NOTE Change here!
1197 {
1198     int g;
1199
1200     if (!userParameters->getEndGapPenalties() && (j == 0 || j == prfLength2))
1201     {
1202         return (0);
1203     }
1204
1205     g = (*profile1)[i][GAPCOL] + (*profile2)[j][GAPCOL];
1206     return (g);
1207 }
1208
1209
1210 /**
1211  *
1212  * @param i
1213  * @param j
1214  * @return
1215  */
1216 inline int MyersMillerProfileAlign::extPenalty2(int i, int j) // NOTE Change here!
1217 {
1218     int h;
1219
1220     if (!userParameters->getEndGapPenalties() && (j == 0 || j == prfLength2))
1221     {
1222         return (0);
1223     }
1224
1225     h = (*profile1)[i][LENCOL];
1226     return (h);
1227 }
1228
1229
1230 /**
1231  *
1232  * @param i
1233  * @param j
1234  * @param k
1235  * @return
1236  */
1237 int MyersMillerProfileAlign::gapPenalty2(int i, int j, int k)
1238 {
1239     int ix;
1240     int gp;
1241     int g, h = 0;
1242
1243     if (k <= 0)
1244     {
1245         return (0);
1246     }
1247     if (!userParameters->getEndGapPenalties() && (j == 0 || j == prfLength2))
1248     {
1249         return (0);
1250     }
1251
1252     g = (*profile1)[i][GAPCOL] + (*profile2)[j][GAPCOL];
1253     for (ix = 0; ix < k && ix + i < prfLength1; ix++)
1254     {
1255         h += (*profile1)[ix + i][LENCOL];
1256     }
1257
1258     gp = g + h;
1259     return (gp);
1260 }
1261
1262 }