+++ /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