+++ /dev/null
-/**
- * Author: Mark Larkin
- *
- * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
- * Changes: Mark 20-6-07, I added a call to calculateMaxlengths
- */
-#ifdef HAVE_CONFIG_H
- #include "config.h"
-#endif
-#include "MyersMillerProfileAlign.h"
-#include "Iteration.h"
-#include <math.h>
-
-namespace clustalw
-{
-
-/**
- *
- * @return
- */
-MyersMillerProfileAlign::MyersMillerProfileAlign()
- : _gapPos1(userParameters->getGapPos1()),
- _gapPos2(userParameters->getGapPos2())
-{
-
-}
-
-/**
- *
- * @param alnPtr
- * @param distMat
- * @param group
- * @param aligned
- * @return
- */
-int MyersMillerProfileAlign::profileAlign(Alignment* alnPtr, DistMatrix* distMat,
- vector<int>* group, int* aligned)
-{
- bool negative = false;
- int i = 0, j = 0, count = 0;
- int numSeq = 0;
- int seqNum = 0;
- int len = 0, len1 = 0, len2 = 0, is = 0, minLen = 0;
- int se1 = 0, se2 = 0, sb1 = 0, sb2 = 0;
- int maxRes = 0;
- int c = 0;
- int score = 0;
- double logmin = 0.0, logdiff = 0.0;
- double pcid = 0.0;
- int matrix[NUMRES][NUMRES];
- int numSeqsProf1 = 0, numSeqsProf2 = 0;
- // Initialise the matrix. To get rid of a valgrind error.
- for(int i = 0; i < NUMRES; i++)
- {
- for(int j = 0; j < NUMRES; j++)
- {
- matrix[i][j] = 0;
- }
- }
-
- numSeq = alnPtr->getNumSeqs();
- seqArray.resize(numSeq);
- alnWeight.resize(numSeq);
-
- for (i = 0; i < numSeq; i++)
- {
- if (aligned[i + 1] == 0)
- {
- (*group)[i + 1] = 0;
- }
- }
-
- nseqs1 = nseqs2 = 0;
- for (i = 0; i < numSeq; i++)
- {
- if ((*group)[i + 1] == 1)
- {
- nseqs1++;
- }
- else if ((*group)[i + 1] == 2)
- {
- nseqs2++;
- }
- }
-
- if ((nseqs1 == 0) || (nseqs2 == 0))
- {
- return (0);
- }
- numSeqsProf1 = nseqs1;
- numSeqsProf2 = nseqs2;
-
- if (nseqs2 > nseqs1)
- {
- switchProfiles = true;
- for (i = 0; i < numSeq; i++)
- {
- if ((*group)[i + 1] == 1)
- {
- (*group)[i + 1] = 2;
- }
- else if ((*group)[i + 1] == 2)
- {
- (*group)[i + 1] = 1;
- }
- }
- }
- else
- {
- switchProfiles = false;
- }
-
- // calculate the mean of the sequence pc identities between the two groups
-
- count = 0;
- pcid = 0.0;
- negative = userParameters->getUseNegMatrix();
-
- for (i = 0; i < numSeq; i++)
- {
- if ((*group)[i + 1] == 1)
- {
- for (j = 0; j < numSeq; j++)
- {
- if ((*group)[j + 1] == 2)
- {
- count++;
- pcid += (*distMat)(i + 1, j + 1);
- }
- }
- }
- }
-
- pcid = pcid / (float)count;
-
- // Make the first profile.
-
- prfLength1 = 0;
- for (i = 0; i < numSeq; i++)
- {
- if ((*group)[i + 1] == 1)
- {
- if (alnPtr->getSeqLength(i + 1) > prfLength1)
- {
- prfLength1 = alnPtr->getSeqLength(i + 1);
- }
- }
- }
- nseqs1 = 0;
-
- for (i = 0; i < numSeq; i++)
- {
- if ((*group)[i + 1] == 1)
- {
- len = alnPtr->getSeqLength(i + 1);
-
- // initialise with the other vector!
- seqArray[nseqs1] = vector<int>(alnPtr->getSequence(i + 1)->begin() + 1,
- alnPtr->getSequence(i + 1)->end());
- seqArray[nseqs1].resize(prfLength1 + extraEndElemNum);
-
- for (j = len; j < prfLength1; j++)
- {
- seqArray[nseqs1][j] = userParameters->getGapPos1();
- }
- alnWeight[nseqs1] = alnPtr->getSeqWeight(i);
- nseqs1++;
- }
- }
-
-
- // Make the second profile.
-
- prfLength2 = 0;
- for (i = 0; i < numSeq; i++)
- {
- if ((*group)[i + 1] == 2)
- {
- if (alnPtr->getSeqLength(i + 1) > prfLength2)
- {
- prfLength2 = alnPtr->getSeqLength(i + 1);
- }
- }
- }
- nseqs2 = 0;
-
- for (i = 0; i < numSeq; i++)
- {
- if ((*group)[i + 1] == 2)
- {
- len = alnPtr->getSeqLength(i + 1);
-
- seqArray[nseqs1 + nseqs2] = vector<int>(alnPtr->getSequence(i + 1)->begin() + 1,
- alnPtr->getSequence(i + 1)->end());
- seqArray[nseqs1 + nseqs2].resize(prfLength2 + extraEndElemNum);
-
- for (j = len; j < prfLength2; j++)
- {
- seqArray[nseqs1 + nseqs2][j] = userParameters->getGapPos1();
- }
-
- seqArray[nseqs1 + nseqs2][j] = ENDALN;
- alnWeight[nseqs1 + nseqs2] = alnPtr->getSeqWeight(i);
- //cout << "weight " << nseqs1 + nseqs2 << alnWeight[nseqs1 + nseqs2] << "\n";
- nseqs2++;
- }
- }
-
- // Change the Max alignment length in the Alignment Object!
- alnPtr->setMaxAlnLength(prfLength1 + prfLength2 + 2);
-
- // calculate real length of profiles - removing gaps!
-
- len1 = 0;
- for (i = 0; i < nseqs1; i++)
- {
- is = 0;
- for (j = 0; j < utilityObject->MIN(static_cast<int>(seqArray[i].size()), prfLength1); j++)
- {
- c = seqArray[i][j];
- if ((c != _gapPos1) && (c != _gapPos2))
- {
- is++;
- }
- }
- len1 += is;
- }
- len1 = static_cast<int>(len1 / (float)nseqs1);
-
- len2 = 0;
- for (i = nseqs1; i < nseqs2 + nseqs1; i++)
- {
- is = 0;
- for (j = 0; j < utilityObject->MIN(static_cast<int>(seqArray[i].size()), prfLength2);
- j++)
- {
- c = seqArray[i][j];
- if ((c != _gapPos1) && (c != _gapPos2))
- {
- is++;
- }
- }
- len2 += is;
- }
- len2 = static_cast<int>(len2 / (float)nseqs2);
-
- PrfScaleValues scaleVals;
- scaleVals.scale = 1.0;
- scaleVals.intScale = 100.0;
-
- int matAvgScore = 0;
- minLen = utilityObject->MIN(len1, len2);
- maxRes = 0;
-
- // Get the substitution matrix that will be stored in 'matrix'
- maxRes = subMatrix->getProfileAlignMatrix(matrix, pcid, minLen, scaleVals, matAvgScore);
-
- if (maxRes == 0 || maxRes == -1)
- {
- return -1;
- }
- if (userParameters->getDNAFlag())
- {
- gapcoef1 = gapcoef2 = static_cast<int>(100.0 * userParameters->getGapOpen() * scaleVals.scale);
- lencoef1 = lencoef2 = static_cast<int>(100.0 * userParameters->getGapExtend() * scaleVals.scale);
- }
- else
- {
- if (len1 == 0 || len2 == 0)
- {
- logmin = 1.0;
- logdiff = 1.0;
- }
- else
- {
- logmin = 1.0 / log10((double)minLen);
- if (len2 < len1)
- {
- logdiff = 1.0 + 0.5 * log10((double)((float)len2 / (float)len1));
- }
- else if (len1 < len2)
- {
- logdiff = 1.0 + 0.5 * log10((double)((float)len1 / (float)len2));
- }
- else
- {
- logdiff = 1.0;
- }
- if (logdiff < 0.9)
- {
- logdiff = 0.9;
- }
- }
-
-
- if (negative)
- {
- gapcoef1 = gapcoef2 = static_cast<int>(100.0 *(float)(userParameters->getGapOpen()));
- lencoef1 = lencoef2 = static_cast<int>(100.0 * userParameters->getGapExtend());
- }
- else
- {
- if (matAvgScore <= 0)
- {
- gapcoef1 = gapcoef2 = static_cast<int>(100.0 *(float)(userParameters->getGapOpen() + logmin));
- }
- else
- {
- gapcoef1 = gapcoef2 = static_cast<int>(scaleVals.scale * matAvgScore *
- (float)(userParameters->getGapOpen() / (logdiff * logmin)));
- }
- lencoef1 = lencoef2 = static_cast<int>(100.0 * userParameters->getGapExtend());
- }
- }
- // We need one profile with substitution matrix information and one without!
- // But this will change when we have the LE scoring function.
- profileWithSub = new ProfileWithSub(prfLength1, 0, nseqs1);
- profileStandard = new ProfileStandard(prfLength2, nseqs1, nseqs1 + nseqs2);
-
- // Calculate the profile array with Substitution matrix info
- // calculate the Gap Coefficients.
-
- gaps.resize(alnPtr->getMaxAlnLength() + 1);
-
- bool profile1Pen;
- if (switchProfiles == false)
- {
- profile1Pen = userParameters->getStructPenalties1() && userParameters->getUseSS1();
- profileWithSub->calcGapCoeff(&seqArray, &gaps, profile1Pen,
- alnPtr->getGapPenaltyMask1(), gapcoef1, lencoef1);
- }
- else
- {
- profile1Pen = userParameters->getStructPenalties2() && userParameters->getUseSS2();
- profileWithSub->calcGapCoeff(&seqArray, &gaps, profile1Pen,
- alnPtr->getGapPenaltyMask2(), gapcoef1, lencoef1);
- }
- // calculate the profile matrix.
- profileWithSub->calcProfileWithSub(&seqArray, &gaps, matrix, &alnWeight);
- profile1 = profileWithSub->getProfilePtr();
-
- if (userParameters->getDebug() > 4)
- {
- string aminoAcidCodes = userParameters->getAminoAcidCodes();
- for (j = 0; j <= userParameters->getMaxAA(); j++)
- {
- cout << aminoAcidCodes[j] << " ";
- }
- cout << "\n";
- for (i = 0; i < prfLength1; i++)
- {
- for (j = 0; j <= userParameters->getMaxAA(); j++)
- {
- cout << (*profile1)[i + 1][j] << " ";
- }
- cout << (*profile1)[i + 1][_gapPos1]<< " ";
- cout << (*profile1)[i + 1][_gapPos2]<< " ";
- cout << (*profile1)[i + 1][GAPCOL] << " " << (*profile1)[i + 1][LENCOL]
- << "\n";
- }
- }
- // Calculate the standard profile array
- // calculate the Gap Coefficients.
- bool profile2Pen;
- if (switchProfiles == false)
- {
- profile2Pen = userParameters->getStructPenalties2() && userParameters->getUseSS2();
- profileStandard->calcGapCoeff(&seqArray, &gaps, profile2Pen,
- alnPtr->getGapPenaltyMask2(), gapcoef2, lencoef2);
- }
- else
- {
- profile2Pen = userParameters->getStructPenalties1() && userParameters->getUseSS1();
- profileStandard->calcGapCoeff(&seqArray, &gaps, profile2Pen,
- alnPtr->getGapPenaltyMask1(), gapcoef2, lencoef2);
- }
-
- // calculate the profile matrix.
-
- profileStandard->calcStandardProfile(&seqArray, &alnWeight);
- profile2 = profileStandard->getProfilePtr();
-
- if (userParameters->getDebug() > 4)
- {
- string aminoAcidCodes = userParameters->getAminoAcidCodes();
- for (j = 0; j <= userParameters->getMaxAA(); j++)
- {
- cout << aminoAcidCodes[j] << " ";
- }
- cout << "\n";
- for (i = 0; i < prfLength2; i++)
- {
- for (j = 0; j <= userParameters->getMaxAA(); j++)
- {
- cout << (*profile2)[i + 1][j] << " ";
- }
- cout << (*profile2)[i + 1][_gapPos1]<< " ";
- cout << (*profile2)[i + 1][_gapPos2]<< " ";
- cout << (*profile2)[i + 1][GAPCOL] << " " << (*profile2)[i + 1][LENCOL]
- << "\n";
- }
- }
- alnWeight.clear();
-
- int _maxAlnLength = alnPtr->getMaxAlnLength();
-
- alnPath1.resize(_maxAlnLength + 1);
- alnPath2.resize(_maxAlnLength + 1);
-
- //align the profiles
-
- // use Myers and Miller to align two sequences
-
- lastPrint = 0;
- printPtr = 1;
-
- sb1 = sb2 = 0;
- se1 = prfLength1;
- se2 = prfLength2;
-
- HH.resize(_maxAlnLength + 1);
- DD.resize(_maxAlnLength + 1);
- RR.resize(_maxAlnLength + 1);
- SS.resize(_maxAlnLength + 1);
- gS.resize(_maxAlnLength + 1);
- displ.resize(_maxAlnLength + 1);
-
- score = progDiff(sb1, sb2, se1 - sb1, se2 - sb2, (*profile1)[0][GAPCOL],
- (*profile1)[prfLength1][GAPCOL]);
-
- // May not need if we recreate the Profile every time!
- HH.clear();
- DD.clear();
- RR.clear();
- SS.clear();
- gS.clear();
-
- alignmentLength = progTracepath();
-
- displ.clear();
-
- addGGaps(alnPtr, &seqArray);
-
- profileStandard->resetProfile();
- profileWithSub->resetProfile();
-
- prfLength1 = alignmentLength;
-
- alnPath1.clear();
- alnPath2.clear();
-
- if(userParameters->getDoRemoveFirstIteration() == TREE)
- {
- Iteration iterateObj;
- iterateObj.iterationOnTreeNode(numSeqsProf1, numSeqsProf2, prfLength1, prfLength2,
- &seqArray);
- }
-
- // Now we resize the SeqArray that holds the sequences in the alignment class
- // and also update it with the new aligned sequences
-
- SeqArray* newSequences = alnPtr->getSeqArrayForRealloc();
- seqNum = 0;
- for (j = 0; j < numSeq; j++)
- {
- if ((*group)[j + 1] == 1)
- {
- (*newSequences)[j + 1].clear();
- (*newSequences)[j + 1].resize(prfLength1 + 1);
- for (i = 0; i < prfLength1; i++)
- {
- (*newSequences)[j + 1][i + 1] = seqArray[seqNum][i];
- }
- seqNum++;
- }
- }
- for (j = 0; j < numSeq; j++)
- {
- if ((*group)[j + 1] == 2)
- {
- (*newSequences)[j + 1].clear();
- (*newSequences)[j + 1].resize(prfLength1 + 1);
- for (i = 0; i < prfLength1; i++)
- {
- (*newSequences)[j + 1][i + 1] = seqArray[seqNum][i];
- }
- seqNum++;
- }
- }
-
- alnPtr->calculateMaxLengths(); // Mark change 20-6-07
-
- gaps.clear();
- seqArray.clear();
-
- delete profileWithSub;
- delete profileStandard;
-
- int retScore = (score / 100);
-
- return retScore;
-
-}
-
-/** ****************************************************************************************
- * Private functions *
- *******************************************************************************************/
-
-/**
- *
- * @param A
- * @param B
- * @param M
- * @param N
- * @param go1
- * @param go2
- * @return
- */
-int MyersMillerProfileAlign::progDiff(int A, int B, int M, int N, int go1, int go2)
-{
- int midi, midj, type;
- int midh;
-
- static int t, tl, g, h;
-
- {
- static int i, j;
- static int hh, f, e, s;
-
- /* Boundary cases: M <= 1 or N == 0 */
- if (userParameters->getDebug() > 2)
- {
- cout << "A " << A << " B " << B << " M " << M << " N " << N
- << " midi " << M / 2 << " go1 " << go1 << " go2 " << go2<< "\n";
- }
-
- /* if sequence B is empty.... */
-
- if (N <= 0)
- {
-
- /* if sequence A is not empty.... */
-
- if (M > 0)
- {
-
- /* delete residues A[1] to A[M] */
-
- progDel(M);
- }
- return ( - gapPenalty1(A, B, M));
- }
-
- /* if sequence A is empty.... */
-
- if (M <= 1)
- {
- if (M <= 0)
- {
-
- /* insert residues B[1] to B[N] */
-
- progAdd(N);
- return ( - gapPenalty2(A, B, N));
- }
-
- /* if sequence A has just one residue.... */
-
- if (go1 == 0)
- {
- midh = - gapPenalty1(A + 1, B + 1, N);
- }
- else
- {
- midh = - gapPenalty2(A + 1, B, 1) - gapPenalty1(A + 1, B + 1,
- N);
- }
- midj = 0;
- for (j = 1; j <= N; j++)
- {
- hh = - gapPenalty1(A, B + 1, j - 1) + prfScore(A + 1, B + j)
- - gapPenalty1(A + 1, B + j + 1, N - j);
- if (hh > midh)
- {
- midh = hh;
- midj = j;
- }
- }
-
- if (midj == 0)
- {
- progAdd(N);
- progDel(1);
- }
- else
- {
- if (midj > 1)
- {
- progAdd(midj - 1);
- }
- progAlign();
- if (midj < N)
- {
- progAdd(N - midj);
- }
- }
- return midh;
- }
-
-
- /* Divide sequence A in half: midi */
-
- midi = M / 2;
-
- /* In a forward phase, calculate all HH[j] and HH[j] */
-
- HH[0] = 0;
- t = - openPenalty1(A, B + 1);
- tl = - extPenalty1(A, B + 1);
- for (j = 1; j <= N; j++)
- {
- HH[j] = t = t + tl;
- DD[j] = t - openPenalty2(A + 1, B + j);
- }
-
- if (go1 == 0)
- {
- t = 0;
- }
- else
- {
- t = - openPenalty2(A + 1, B);
- }
- tl = - extPenalty2(A + 1, B);
- for (i = 1; i <= midi; i++)
- {
- s = HH[0];
- HH[0] = hh = t = t + tl;
- f = t - openPenalty1(A + i, B + 1);
-
- for (j = 1; j <= N; j++)
- {
- g = openPenalty1(A + i, B + j);
- h = extPenalty1(A + i, B + j);
- if ((hh = hh - g - h) > (f = f - h))
- {
- f = hh;
- }
- g = openPenalty2(A + i, B + j);
- h = extPenalty2(A + i, B + j);
- if ((hh = HH[j] - g - h) > (e = DD[j] - h))
- {
- e = hh;
- }
- hh = s + prfScore(A + i, B + j);
- if (f > hh)
- {
- hh = f;
- }
- if (e > hh)
- {
- hh = e;
- }
-
- s = HH[j];
- HH[j] = hh;
- DD[j] = e;
-
- }
- }
-
- DD[0] = HH[0];
-
- /* In a reverse phase, calculate all RR[j] and SS[j] */
-
- RR[N] = 0;
- tl = 0;
- for (j = N - 1; j >= 0; j--)
- {
- g = - openPenalty1(A + M, B + j + 1);
- tl -= extPenalty1(A + M, B + j + 1);
- RR[j] = g + tl;
- SS[j] = RR[j] - openPenalty2(A + M, B + j);
- gS[j] = openPenalty2(A + M, B + j);
- }
-
- tl = 0;
- for (i = M - 1; i >= midi; i--)
- {
- s = RR[N];
- if (go2 == 0)
- {
- g = 0;
- }
- else
- {
- g = - openPenalty2(A + i + 1, B + N);
- }
- tl -= extPenalty2(A + i + 1, B + N);
- RR[N] = hh = g + tl;
- t = openPenalty1(A + i, B + N);
- f = RR[N] - t;
-
- for (j = N - 1; j >= 0; j--)
- {
- g = openPenalty1(A + i, B + j + 1);
- h = extPenalty1(A + i, B + j + 1);
- if ((hh = hh - g - h) > (f = f - h - g + t))
- {
- f = hh;
- }
- t = g;
- g = openPenalty2(A + i + 1, B + j);
- h = extPenalty2(A + i + 1, B + j);
- hh = RR[j] - g - h;
- if (i == (M - 1))
- {
- e = SS[j] - h;
- }
- else
- {
- e = SS[j] - h - g + openPenalty2(A + i + 2, B + j);
- gS[j] = g;
- }
- if (hh > e)
- {
- e = hh;
- }
- hh = s + prfScore(A + i + 1, B + j + 1);
- if (f > hh)
- {
- hh = f;
- }
- if (e > hh)
- {
- hh = e;
- }
-
- s = RR[j];
- RR[j] = hh;
- SS[j] = e;
-
- }
- }
- SS[N] = RR[N];
- gS[N] = openPenalty2(A + midi + 1, B + N);
-
- /* find midj, such that HH[j]+RR[j] or DD[j]+SS[j]+gap is the maximum */
-
- midh = HH[0] + RR[0];
- midj = 0;
- type = 1;
- for (j = 0; j <= N; j++)
- {
- hh = HH[j] + RR[j];
- if (hh >= midh)
- if (hh > midh || (HH[j] != DD[j] && RR[j] == SS[j]))
- {
- midh = hh;
- midj = j;
- }
- }
-
- for (j = N; j >= 0; j--)
- {
- hh = DD[j] + SS[j] + gS[j];
- if (hh > midh)
- {
- midh = hh;
- midj = j;
- type = 2;
- }
- }
- }
-
- /* Conquer recursively around midpoint */
-
-
- if (type == 1)
- {
- /* Type 1 gaps */
- if (userParameters->getDebug() > 2)
- {
- cout << "Type 1,1: midj " << midj << "\n";
- }
- progDiff(A, B, midi, midj, go1, 1);
- if (userParameters->getDebug() > 2)
- {
- cout << "Type 1,2: midj " << midj << "\n";
- }
- progDiff(A + midi, B + midj, M - midi, N - midj, 1, go2);
- }
- else
- {
- if (userParameters->getDebug() > 2)
- {
- cout << "Type 2,1: midj " << midj << "\n";
- }
- progDiff(A, B, midi - 1, midj, go1, 0);
- progDel(2);
- if (userParameters->getDebug() > 2)
- {
- cout << "Type 2,2: midj " << midj << "\n";
- }
- progDiff(A + midi + 1, B + midj, M - midi - 1, N - midj, 0, go2);
- }
-
- return midh; /* Return the score of the best alignment */
-}
-
-/**
- *
- * @param alnPtr
- * @param seqArray
- */
-void MyersMillerProfileAlign::addGGaps(Alignment* alnPtr, SeqArray* seqArray)
-{
- int j;
- int i, ix;
- int len;
- vector<int> ta;
-
- ta.resize(alignmentLength + 1);
-
- for (j = 0; j < nseqs1; j++)
- {
- ix = 0;
- for (i = 0; i < alignmentLength; i++)
- {
- if (alnPath1[i] == 2)
- {
- if (ix < ((int)(*seqArray)[j].size() - extraEndElemNum))
- {
- ta[i] = (*seqArray)[j][ix];
- }
- else
- {
- ta[i] = ENDALN;
- }
- ix++;
- }
- else if (alnPath1[i] == 1)
- {
- /*
- insertion in first alignment...
- */
- ta[i] = _gapPos1;
- }
- else
- {
- cerr << "Error in aln_path\n";
- }
- }
- ta[i] = ENDALN;
-
- len = alignmentLength;
-
- (*seqArray)[j].resize(len + 2);
-
- for (i = 0; i < len; i++)
- {
- (*seqArray)[j][i] = ta[i];
- }
- (*seqArray)[j][len] = ENDALN;
-
- }
-
- for (j = nseqs1; j < nseqs1 + nseqs2; j++)
- {
- ix = 0;
- for (i = 0; i < alignmentLength; i++)
- {
- if (alnPath2[i] == 2)
- {
- if (ix < ((int)(*seqArray)[j].size() - extraEndElemNum))
- {
- ta[i] = (*seqArray)[j][ix];
- }
- else
- {
- ta[i] = ENDALN;
- }
- ix++;
- }
- else if (alnPath2[i] == 1)
- {
- /*
- insertion in second alignment...
- */
- ta[i] = _gapPos1;
- }
- else
- {
- cerr << "Error in alnPath\n";
- }
- }
- ta[i] = ENDALN;
-
- len = alignmentLength;
-
- (*seqArray)[j].resize(len + 2);
-
- for (i = 0; i < len; i++)
- {
- (*seqArray)[j][i] = ta[i];
- }
- (*seqArray)[j][len] = ENDALN;
- }
-
-
- if (userParameters->getStructPenalties1() != NONE)
- {
- addGGapsMask(alnPtr->getGapPenaltyMask1(), alignmentLength,
- &alnPath1, &alnPath2);
- }
- if (userParameters->getStructPenalties1() == SECST)
- {
- addGGapsMask(alnPtr->getSecStructMask1(), alignmentLength,
- &alnPath1, &alnPath2);
- }
-
- if (userParameters->getStructPenalties2() != NONE)
- {
- addGGapsMask(alnPtr->getGapPenaltyMask2(), alignmentLength,
- &alnPath2, &alnPath1);
- }
- if (userParameters->getStructPenalties2() == SECST)
- {
- addGGapsMask(alnPtr->getSecStructMask2(), alignmentLength,
- &alnPath2, &alnPath1);
- }
-}
-
-
-/**
- *
- * @param mask
- * @param len
- * @param path1
- * @param path2
- */
-void MyersMillerProfileAlign::addGGapsMask(vector<char>* mask, int len, vector<int>* path1,
- vector<int>* path2)
-{
- int i, ix;
- char *ta;
-
- ta = new char[len + 1];
-
- ix = 0;
- if (switchProfiles == false)
- {
- for (i = 0; i < len; i++)
- {
- if ((*path1)[i] == 2)
- {
- ta[i] = (*mask)[ix];
- ix++;
- }
- else if ((*path1)[i] == 1)
- {
- ta[i] = _gapPos1;
- }
- }
- }
- else
- {
- for (i = 0; i < len; i++)
- {
- if ((*path2)[i] == 2)
- {
- ta[i] = (*mask)[ix];
- ix++;
- }
- else if ((*path2)[i] == 1)
- {
- ta[i] = _gapPos1;
- }
- }
- }
- mask->resize(len + 2);
-
- for (i = 0; i < len; i++)
- {
- (*mask)[i] = ta[i];
- }
-
-
- delete [] ta;
- ta = NULL;
-}
-
-
-/**
- *
- * @param n
- * @param m
- * @return
- */
-inline int MyersMillerProfileAlign::prfScore(int n, int m)
-{
- int ix;
- int score;
-
- score = 0;
- int _maxAA = userParameters->getMaxAA(); // NOTE Change here!
- for (ix = 0; ix <= _maxAA; ix++)
- {
- score += ((*profile1)[n][ix] * (*profile2)[m][ix]);
- }
- score += ((*profile1)[n][_gapPos1] * (*profile2)[m][_gapPos1]);
- score += ((*profile1)[n][_gapPos2] * (*profile2)[m][_gapPos2]);
- return (score / 10);
-}
-
-
-/**
- *
- * @return
- */
-int MyersMillerProfileAlign::progTracepath()
-{
- int i, j, k, pos, toDo;
- int alignLen;
- pos = 0;
-
- toDo = printPtr - 1;
-
- for (i = 1; i <= toDo; ++i)
- {
- if (userParameters->getDebug() > 1)
- {
- cout << displ[i] << " ";
- }
- if (displ[i] == 0)
- {
- alnPath1[pos] = 2;
- alnPath2[pos] = 2;
- ++pos;
- }
- else
- {
- if ((k = displ[i]) > 0)
- {
- for (j = 0; j <= k - 1; ++j)
- {
- alnPath2[pos + j] = 2;
- alnPath1[pos + j] = 1;
- }
- pos += k;
- }
- else
- {
- k = (displ[i] < 0) ? displ[i] * - 1: displ[i];
- for (j = 0; j <= k - 1; ++j)
- {
- alnPath1[pos + j] = 2;
- alnPath2[pos + j] = 1;
- }
- pos += k;
- }
- }
- }
- if (userParameters->getDebug() > 1)
- {
- cout << "\n";
- }
-
- alignLen = pos;
- return alignLen;
-}
-
-
-/**
- *
- * @param k
- */
-void MyersMillerProfileAlign::progDel(int k)
-{
- if (lastPrint < 0)
- {
- lastPrint = displ[printPtr - 1] -= k;
- }
- else
- {
- lastPrint = displ[printPtr++] = - (k);
- }
-}
-
-
-/**
- *
- * @param k
- */
-void MyersMillerProfileAlign::progAdd(int k)
-{
- if (lastPrint < 0)
- {
- displ[printPtr - 1] = k;
- displ[printPtr++] = lastPrint;
- }
- else
- {
- lastPrint = displ[printPtr++] = k;
- }
-}
-
-
-/**
- *
- */
-void MyersMillerProfileAlign::progAlign()
-{
- displ[printPtr++] = lastPrint = 0;
-}
-
-
-/**
- *
- * @param i
- * @param j
- * @return
- */
-inline int MyersMillerProfileAlign::openPenalty1(int i, int j) // NOTE Change here!
-{
- int g;
- if (!userParameters->getEndGapPenalties() && (i == 0 || i == prfLength1))
- {
- return (0);
- }
-
- g = (*profile2)[j][GAPCOL] + (*profile1)[i][GAPCOL];
- return (g);
-}
-
-
-/**
- *
- * @param i
- * @param j
- * @return
- */
-inline int MyersMillerProfileAlign::extPenalty1(int i, int j) // NOTE Change here!
-{
- int h;
-
- if (!userParameters->getEndGapPenalties() && (i == 0 || i == prfLength1))
- {
- return (0);
- }
-
- h = (*profile2)[j][LENCOL];
- return (h);
-}
-
-
-/**
- *
- * @param i
- * @param j
- * @param k
- * @return
- */
-int MyersMillerProfileAlign::gapPenalty1(int i, int j, int k)
-{
- int ix;
- int gp;
- int g, h = 0;
-
- if (k <= 0)
- {
- return (0);
- }
- if (!userParameters->getEndGapPenalties() && (i == 0 || i == prfLength1))
- {
- return (0);
- }
-
- g = (*profile2)[j][GAPCOL] + (*profile1)[i][GAPCOL];
- for (ix = 0; ix < k && ix + j < prfLength2; ix++)
- {
- h += (*profile2)[ix + j][LENCOL];
- }
-
- gp = g + h;
- return (gp);
-}
-
-
-/**
- *
- * @param i
- * @param j
- * @return
- */
-inline int MyersMillerProfileAlign::openPenalty2(int i, int j) // NOTE Change here!
-{
- int g;
-
- if (!userParameters->getEndGapPenalties() && (j == 0 || j == prfLength2))
- {
- return (0);
- }
-
- g = (*profile1)[i][GAPCOL] + (*profile2)[j][GAPCOL];
- return (g);
-}
-
-
-/**
- *
- * @param i
- * @param j
- * @return
- */
-inline int MyersMillerProfileAlign::extPenalty2(int i, int j) // NOTE Change here!
-{
- int h;
-
- if (!userParameters->getEndGapPenalties() && (j == 0 || j == prfLength2))
- {
- return (0);
- }
-
- h = (*profile1)[i][LENCOL];
- return (h);
-}
-
-
-/**
- *
- * @param i
- * @param j
- * @param k
- * @return
- */
-int MyersMillerProfileAlign::gapPenalty2(int i, int j, int k)
-{
- int ix;
- int gp;
- int g, h = 0;
-
- if (k <= 0)
- {
- return (0);
- }
- if (!userParameters->getEndGapPenalties() && (j == 0 || j == prfLength2))
- {
- return (0);
- }
-
- g = (*profile1)[i][GAPCOL] + (*profile2)[j][GAPCOL];
- for (ix = 0; ix < k && ix + i < prfLength1; ix++)
- {
- h += (*profile1)[ix + i][LENCOL];
- }
-
- gp = g + h;
- return (gp);
-}
-
-}
git://source.jalview.org / jabaws.git / blobdiff