--- /dev/null
+/**
+ * Author: Mark Larkin
+ *
+ * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
+ */
+#ifdef HAVE_CONFIG_H
+ #include "config.h"
+#endif
+#include "FullPairwiseAlign.h"
+#include <math.h>
+
+namespace clustalw
+{
+
+FullPairwiseAlign::FullPairwiseAlign()
+: _maxAlnLength(0),
+ intScale(0),
+ mmScore(0),
+ printPtr(0),
+ lastPrint(0),
+ _gapOpen(0),
+ _gapExtend(0),
+ seq1(0),
+ seq2(0),
+ maxScore(0),
+ sb1(0),
+ sb2(0),
+ se1(0),
+ se2(0)
+{
+
+}
+
+void FullPairwiseAlign::pairwiseAlign(Alignment *alignPtr, DistMatrix *distMat, int iStart, int iEnd, int jStart, int jEnd)
+{
+ int si, sj, i;
+ int n, m, len1, len2;
+ int maxRes;
+ int _matAvgScore;
+ int res;
+ double _score;
+ float gapOpenScale, gapExtendScale;
+
+ try
+ {
+
+ if(distMat->getSize() != alignPtr->getNumSeqs() + 1)
+ {
+ cerr << "The distance matrix is not the right size!\n"
+ << "Need to terminate program.\n";
+ exit(1);
+ }
+ if((iStart < 0) || (iEnd < iStart) || (jStart < 0) || (jEnd < jStart))
+ {
+ cerr << "The range for pairwise Alignment is incorrect.\n"
+ << "Need to terminate program.\n";
+ exit(1);
+ }
+
+ _maxAlnLength = alignPtr->getMaxAlnLength();
+
+ int _numSeqs = alignPtr->getNumSeqs();
+ if(_numSeqs == 0)
+ {
+ return;
+ }
+
+ int num = (2 * _maxAlnLength) + 1;
+ bool _DNAFlag = userParameters->getDNAFlag();
+ float _pwGapOpen, _pwGapExtend;
+ _pwGapOpen = userParameters->getPWGapOpen();
+ _pwGapExtend = userParameters->getPWGapExtend();
+
+ displ.resize(num);
+ HH.resize(_maxAlnLength);
+ DD.resize(_maxAlnLength);
+ RR.resize(_maxAlnLength);
+ SS.resize(_maxAlnLength);
+ // Note these 2 lines replace the stuff above because it is all done in the SubMatrix
+ PairScaleValues scaleValues;
+ maxRes = subMatrix->getPairwiseMatrix(matrix, scaleValues, _matAvgScore);
+ if (maxRes == 0)
+ {
+ cerr << "Could not get the substitution matrix\n";
+ return;
+ }
+
+ intScale = scaleValues.intScale;
+ gapOpenScale = scaleValues.gapOpenScale;
+ gapExtendScale = scaleValues.gapExtendScale;
+
+ int _gapPos1, _gapPos2;
+ _gapPos1 = userParameters->getGapPos1();
+ _gapPos2 = userParameters->getGapPos2();
+ const SeqArray* _ptrToSeqArray = alignPtr->getSeqArray(); //This is faster!
+
+ for (si = utilityObject->MAX(0, iStart); si < _numSeqs && si < iEnd; si++)
+ {
+ n = alignPtr->getSeqLength(si + 1);
+ len1 = 0;
+ for (i = 1; i <= n; i++)
+ {
+ res = (*_ptrToSeqArray)[si + 1][i];
+ if ((res != _gapPos1) && (res != _gapPos2))
+ {
+ len1++;
+ }
+ }
+
+ for (sj = utilityObject->MAX(si+1, jStart+1); sj < _numSeqs && sj < jEnd; sj++)
+ {
+ m = alignPtr->getSeqLength(sj + 1);
+ if (n == 0 || m == 0)
+ {
+ distMat->SetAt(si + 1, sj + 1, 1.0);
+ distMat->SetAt(sj + 1, si + 1, 1.0);
+ continue;
+ }
+ len2 = 0;
+ for (i = 1; i <= m; i++)
+ {
+ res = (*_ptrToSeqArray)[sj + 1][i];
+ if ((res != _gapPos1) && (res != _gapPos2))
+ {
+ len2++;
+ }
+ }
+
+ if (_DNAFlag)
+ {
+ _gapOpen = static_cast<int>(2 * _pwGapOpen * intScale *
+ gapOpenScale);
+ _gapExtend = static_cast<int>(_pwGapExtend * intScale * gapExtendScale);
+ }
+ else
+ {
+ if (_matAvgScore <= 0)
+ {
+ _gapOpen = 2 * static_cast<int>((_pwGapOpen +
+ log(static_cast<double>(utilityObject->MIN(n, m)))) * intScale);
+ }
+ else
+ {
+ _gapOpen = static_cast<int>(2 * _matAvgScore * (_pwGapOpen +
+ log(static_cast<double>(utilityObject->MIN(n, m)))) * gapOpenScale);
+ }
+ _gapExtend = static_cast<int>(_pwGapExtend * intScale);
+ }
+ // align the sequences
+
+ seq1 = si + 1;
+ seq2 = sj + 1;
+
+ _ptrToSeq1 = alignPtr->getSequence(seq1);
+ _ptrToSeq2 = alignPtr->getSequence(seq2);
+
+ forwardPass(_ptrToSeq1, _ptrToSeq2, n, m);
+ reversePass(_ptrToSeq1, _ptrToSeq2);
+
+ lastPrint = 0;
+ printPtr = 1;
+
+ // use Myers and Miller to align two sequences
+
+ maxScore = diff(sb1 - 1, sb2 - 1, se1 - sb1 + 1, se2 - sb2 + 1,
+ (int)0, (int)0);
+
+ // calculate percentage residue identity
+
+ mmScore = tracePath(sb1, sb2);
+
+ if (len1 == 0 || len2 == 0)
+ {
+ mmScore = 0;
+ }
+ else
+ {
+ mmScore /= (float)utilityObject->MIN(len1, len2);
+ }
+
+ _score = ((float)100.0 - mmScore) / (float)100.0;
+ distMat->SetAt(si + 1, sj + 1, _score);
+ distMat->SetAt(sj + 1, si + 1, _score);
+
+ if(userParameters->getDisplayInfo())
+ {
+ utilityObject->info("Sequences (%d:%d) Aligned. Score: %d",
+ si+1, sj+1, (int)mmScore);
+ }
+ }
+ }
+
+ displ.clear();
+ HH.clear();
+ DD.clear();
+ RR.clear();
+ SS.clear();
+ }
+ catch(const exception& e)
+ {
+ cerr << "An exception has occured in the FullPairwiseAlign class.\n"
+ << e.what() << "\n";
+ exit(1);
+ }
+}
+
+void FullPairwiseAlign::add(int v)
+{
+ if (lastPrint < 0)
+ {
+ displ[printPtr - 1] = v;
+ displ[printPtr++] = lastPrint;
+ }
+ else
+ {
+ lastPrint = displ[printPtr++] = v;
+ }
+}
+
+inline int FullPairwiseAlign::calcScore(int iat, int jat, int v1, int v2)
+{
+ return matrix[(*_ptrToSeq1)[v1 + iat]][(*_ptrToSeq2)[v2 + jat]];
+}
+
+float FullPairwiseAlign::tracePath(int tsb1, int tsb2)
+{
+ int res1, res2;
+ int i1, i2;
+ int i, k, pos, toDo;
+ int count;
+ float score;
+
+ toDo = printPtr - 1;
+ i1 = tsb1;
+ i2 = tsb2;
+
+ pos = 0;
+ count = 0;
+ for (i = 1; i <= toDo; ++i)
+ {
+ if (displ[i] == 0)
+ {
+ res1 = (*_ptrToSeq1)[i1];
+ res2 = (*_ptrToSeq2)[i2];
+
+ if ((res1 != userParameters->getGapPos1()) &&
+ (res2 != userParameters->getGapPos2()) && (res1 == res2))
+ {
+ count++;
+ }
+ ++i1;
+ ++i2;
+ ++pos;
+ }
+ else
+ {
+ if ((k = displ[i]) > 0)
+ {
+ i2 += k;
+ pos += k;
+ }
+ else
+ {
+ i1 -= k;
+ pos -= k;
+ }
+ }
+ }
+
+ score = 100.0 *(float)count;
+ return (score);
+}
+
+void FullPairwiseAlign::forwardPass(const vector<int>* seq1, const vector<int>* seq2, int n, int m)
+{
+ int i, j;
+ int f, hh, p, t;
+
+ maxScore = 0;
+ se1 = se2 = 0;
+ for (i = 0; i <= m; i++)
+ {
+ HH[i] = 0;
+ DD[i] = -_gapOpen;
+ }
+
+ for (i = 1; i <= n; i++)
+ {
+ hh = p = 0;
+ f = -_gapOpen;
+
+ for (j = 1; j <= m; j++)
+ {
+
+ f -= _gapExtend;
+ t = hh - _gapOpen - _gapExtend;
+ if (f < t)
+ {
+ f = t;
+ }
+
+ DD[j] -= _gapExtend;
+ t = HH[j] - _gapOpen - _gapExtend;
+ if (DD[j] < t)
+ {
+ DD[j] = t;
+ }
+
+ hh = p + matrix[(*seq1)[i]][(*seq2)[j]];
+ if (hh < f)
+ {
+ hh = f;
+ }
+ if (hh < DD[j])
+ {
+ hh = DD[j];
+ }
+ if (hh < 0)
+ {
+ hh = 0;
+ }
+
+ p = HH[j];
+ HH[j] = hh;
+
+ if (hh > maxScore)
+ {
+ maxScore = hh;
+ se1 = i;
+ se2 = j;
+ }
+ }
+ }
+
+}
+
+void FullPairwiseAlign::reversePass(const vector<int>* seq1, const vector<int>* seq2)
+{
+ int i, j;
+ int f, hh, p, t;
+ int cost;
+
+ cost = 0;
+ sb1 = sb2 = 1;
+ for (i = se2; i > 0; i--)
+ {
+ HH[i] = - 1;
+ DD[i] = - 1;
+ }
+
+ for (i = se1; i > 0; i--)
+ {
+ hh = f = - 1;
+ if (i == se1)
+ {
+ p = 0;
+ }
+ else
+ {
+ p = - 1;
+ }
+
+ for (j = se2; j > 0; j--)
+ {
+
+ f -= _gapExtend;
+ t = hh - _gapOpen - _gapExtend;
+ if (f < t)
+ {
+ f = t;
+ }
+
+ DD[j] -= _gapExtend;
+ t = HH[j] - _gapOpen - _gapExtend;
+ if (DD[j] < t)
+ {
+ DD[j] = t;
+ }
+
+ hh = p + matrix[(*seq1)[i]][(*seq2)[j]];
+ if (hh < f)
+ {
+ hh = f;
+ }
+ if (hh < DD[j])
+ {
+ hh = DD[j];
+ }
+
+ p = HH[j];
+ HH[j] = hh;
+
+ if (hh > cost)
+ {
+ cost = hh;
+ sb1 = i;
+ sb2 = j;
+ if (cost >= maxScore)
+ {
+ break;
+ }
+ }
+ }
+ if (cost >= maxScore)
+ {
+ break;
+ }
+ }
+
+}
+
+int FullPairwiseAlign::diff(int A, int B, int M, int N, int tb, int te)
+{
+ int type;
+ int midi, midj, i, j;
+ int midh;
+ static int f, hh, e, s, t;
+
+ if (N <= 0)
+ {
+ if (M > 0)
+ {
+ del(M);
+ }
+
+ return ( -(int)tbgap(M, tb));
+ }
+
+ if (M <= 1)
+ {
+ if (M <= 0)
+ {
+ add(N);
+ return ( -(int)tbgap(N, tb));
+ }
+
+ midh = - (tb + _gapExtend) - tegap(N, te);
+ hh = - (te + _gapExtend) - tbgap(N, tb);
+ if (hh > midh)
+ {
+ midh = hh;
+ }
+ midj = 0;
+ for (j = 1; j <= N; j++)
+ {
+ hh = calcScore(1, j, A, B) - tegap(N - j, te) - tbgap(j - 1, tb);
+ if (hh > midh)
+ {
+ midh = hh;
+ midj = j;
+ }
+ }
+
+ if (midj == 0)
+ {
+ del(1);
+ add(N);
+ }
+ else
+ {
+ if (midj > 1)
+ {
+ add(midj - 1);
+ }
+ displ[printPtr++] = lastPrint = 0;
+ if (midj < N)
+ {
+ add(N - midj);
+ }
+ }
+ return midh;
+ }
+
+ // Divide: Find optimum midpoint (midi,midj) of cost midh
+
+ midi = M / 2;
+ HH[0] = 0;
+ t = - tb;
+ for (j = 1; j <= N; j++)
+ {
+ HH[j] = t = t - _gapExtend;
+ DD[j] = t - _gapOpen;
+ }
+
+ t = - tb;
+ for (i = 1; i <= midi; i++)
+ {
+ s = HH[0];
+ HH[0] = hh = t = t - _gapExtend;
+ f = t - _gapOpen;
+ for (j = 1; j <= N; j++)
+ {
+ if ((hh = hh - _gapOpen - _gapExtend) > (f = f - _gapExtend))
+ {
+ f = hh;
+ }
+ if ((hh = HH[j] - _gapOpen - _gapExtend) > (e = DD[j] - _gapExtend))
+ {
+ e = hh;
+ }
+ hh = s + calcScore(i, j, A, B);
+ if (f > hh)
+ {
+ hh = f;
+ }
+ if (e > hh)
+ {
+ hh = e;
+ }
+
+ s = HH[j];
+ HH[j] = hh;
+ DD[j] = e;
+ }
+ }
+
+ DD[0] = HH[0];
+
+ RR[N] = 0;
+ t = - te;
+ for (j = N - 1; j >= 0; j--)
+ {
+ RR[j] = t = t - _gapExtend;
+ SS[j] = t - _gapOpen;
+ }
+
+ t = - te;
+ for (i = M - 1; i >= midi; i--)
+ {
+ s = RR[N];
+ RR[N] = hh = t = t - _gapExtend;
+ f = t - _gapOpen;
+
+ for (j = N - 1; j >= 0; j--)
+ {
+
+ if ((hh = hh - _gapOpen - _gapExtend) > (f = f - _gapExtend))
+ {
+ f = hh;
+ }
+ if ((hh = RR[j] - _gapOpen - _gapExtend) > (e = SS[j] - _gapExtend))
+ {
+ e = hh;
+ }
+ hh = s + calcScore(i + 1, j + 1, A, B);
+ if (f > hh)
+ {
+ hh = f;
+ }
+ if (e > hh)
+ {
+ hh = e;
+ }
+
+ s = RR[j];
+ RR[j] = hh;
+ SS[j] = e;
+
+ }
+ }
+
+ SS[N] = RR[N];
+
+ 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] + _gapOpen;
+ if (hh > midh)
+ {
+ midh = hh;
+ midj = j;
+ type = 2;
+ }
+ }
+
+ // Conquer recursively around midpoint
+
+
+ if (type == 1)
+ {
+ // Type 1 gaps
+ diff(A, B, midi, midj, tb, _gapOpen);
+ diff(A + midi, B + midj, M - midi, N - midj, _gapOpen, te);
+ }
+ else
+ {
+ diff(A, B, midi - 1, midj, tb, 0);
+ del(2);
+ diff(A + midi + 1, B + midj, M - midi - 1, N - midj, 0, te);
+ }
+
+ return midh; // Return the score of the best alignment
+}
+
+void FullPairwiseAlign::del(int k)
+{
+ if (lastPrint < 0)
+ {
+ lastPrint = displ[printPtr - 1] -= k;
+ }
+ else
+ {
+ lastPrint = displ[printPtr++] = - (k);
+ }
+}
+
+int FullPairwiseAlign::gap(int k)
+{
+ if(k <= 0)
+ {
+ return 0;
+ }
+ else
+ {
+ return _gapOpen + _gapExtend * k;
+ }
+}
+
+int FullPairwiseAlign::tbgap(int k, int tb)
+{
+ if(k <= 0)
+ {
+ return 0;
+ }
+ else
+ {
+ return tb + _gapExtend * k;
+ }
+}
+
+int FullPairwiseAlign::tegap(int k, int te)
+{
+ if(k <= 0)
+ {
+ return 0;
+ }
+ else
+ {
+ return te + _gapExtend * k;
+ }
+}
+
+}
git://source.jalview.org / jabaws.git / blobdiff