+++ /dev/null
-/**
- * Author: Mark Larkin
- *
- * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
- */
-#ifdef HAVE_CONFIG_H
- #include "config.h"
-#endif
-#include <math.h>
-#include "FastPairwiseAlign.h"
-
-namespace clustalw
-{
-
-FastPairwiseAlign::FastPairwiseAlign()
-{
- _maxAlnLength = 0;
-}
-
-void FastPairwiseAlign::pairwiseAlign(Alignment *alignPtr, DistMatrix *distMat, int iStart,
- int iEnd, int jStart, int jEnd)
-{
- 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))
- {
- cout << "The range for pairwise Alignment is incorrect.\n"
- << "Need to terminate program.\n";
- exit(1);
- }
-
- int i, j, dsr;
- double calcScore;
- bool _DNAFlag = userParameters->getDNAFlag();
- _maxAlnLength = alignPtr->getMaxAlnLength();
- int num = (2 * _maxAlnLength) + 1;
- accum.ResizeRect(5, num);
-
- displ.resize(num);
- slopes.resize(num);
- diagIndex.resize(num);
-
- zza.resize(_maxAlnLength + 1);
- zzb.resize(_maxAlnLength + 1);
- zzc.resize(_maxAlnLength + 1);
- zzd.resize(_maxAlnLength + 1);
-
- if (_DNAFlag)
- {
- userParameters->setDNAParams();
- }
- else
- {
- userParameters->setProtParams();
- }
-
- cout << "\n\n";
-
- for (i = iStart + 1; i <= iEnd; ++i)
- {
- const vector<int>* _seqIPtr = alignPtr->getSequence(i);
- int _seqILength = alignPtr->getSeqLength(i);
- if (_DNAFlag)
- {
- makeNPtrs(zza, zzc, _seqIPtr, _seqILength);
- }
- else
- {
- makePPtrs(zza, zzc, _seqIPtr, _seqILength);
- }
- double _score;
- for (j = jStart + 2 > i+1 ? jStart + 2 : i+1; j <= jEnd; ++j)
- {
- const vector<int>* _seqJPtr = alignPtr->getSequence(j);
- int _seqJLength = alignPtr->getSeqLength(j);
- if (_DNAFlag)
- {
- makeNPtrs(zzb, zzd, _seqJPtr, _seqJLength);
- }
- else
- {
- makePPtrs(zzb, zzd, _seqJPtr, _seqJLength);
- }
- pairAlign(_seqIPtr, _seqILength, _seqJLength);
- if (!maxSoFar)
- {
- calcScore = 0.0;
- }
- else
- {
- calcScore = (double)accum[0][maxSoFar];
- if (userParameters->getPercent())
- {
- dsr = (_seqILength < _seqJLength) ? _seqILength
- : _seqJLength;
- calcScore = (calcScore / (double)dsr) *100.0;
- }
- }
- _score = (100.0 - calcScore) / 100.0;
- distMat->SetAt(i, j, _score);
- //distMat->SetAt(j, i, _score); /* distMat symmetric, FS, 2009-04-06 */
-
-
- if(userParameters->getDisplayInfo())
- {
- if (calcScore > 0.1)
- {
- utilityObject->info("Sequences (%d:%d) Aligned. Score: %lg",
- i, j, calcScore);
- }
- else
- {
- utilityObject->info("Sequences (%d:%d) Not Aligned", i, j);
- }
- }
- }
- }
- accum.clearArray();
- displ.clear();
- slopes.clear();
- diagIndex.clear();
-
- zza.clear();
- zzb.clear();
- zzc.clear();
- zzd.clear();
- }
- catch(const exception& e)
- {
- cerr << "An exception has occured in the FastPairwiseAlign class.\n"
- << e.what() << "\n";
- exit(1);
- }
-}
-
-
-/*
- * Note: There is a problem with the treatment of DNA/RNA.
- * During file reading all residues are encoded as AminoAcids,
- * even before it has been established if they are AA or not.
- * This is bad and will have to be changed (later).
- * 'A' is assigned code 0, C is assigned 2, G = 6, T=18, U=19.
- * However, the fast alignment routines require that
- * A=0, C=1, G=2, T=U=3. In the best case the results of the
- * fast alignment (of DNA) will simply be meaningless, the
- * worst case is a core dump.
- * As a quick fix I implemented the following mask, that
- * (for DNA/RNA) translates 0->0, 2->1, 6->2, 18->3, 19->3.
- * This is awfull (it is not OO compliant) but it was quick,
- * uses (much) less memory than a second DNA array, and is
- * (much) faster than calling a translation function.
- * Ideally this will be removed, but this requires changes
- * to (i) the sequence encoding during file-reading AND
- * (ii) all the DNA substitution matrices. (Fabian, 2009-02-25)
- *
- * A B C D E F G H I K L M 0 N P Q R S
- * T U W X Y Z (goodmeasuregoodmeasuregood) */
-int ziAA2DNA[] = {0,-1, 1,-1,-1,-1, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
- 3, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
-const int ciMaxResID = 3;
-
-
-
-
-void FastPairwiseAlign::pairAlign(const vector<int>* seq, int l1, int l2)
-{
- int pot[8], i, j, l, m, limit, pos, tl1, vn1, vn2, flen, osptr, fs;
- int tv1, tv2, encrypt, subt1, subt2, rmndr;
- bool flag;
- int residue;
- bool _DNAFlag = userParameters->getDNAFlag();
- int _ktup = userParameters->getKtup();
- int _maxAA = userParameters->getMaxAA();
- int _windowGap = userParameters->getWindowGap();
- int _window = userParameters->getWindow();
- int _signif = userParameters->getSignif();
-
- if (_DNAFlag)
- {
- for (i = 1; i <= _ktup; ++i)
- {
- pot[i] = (int)pow((double)4, (double)(i - 1));
- }
- limit = (int)pow((double)4, (double)_ktup);
- }
- else
- {
- for (i = 1; i <= _ktup; i++)
- {
- pot[i] = (int)pow((double)(_maxAA + 1), (double)(i - 1));
- }
- limit = (int)pow((double)(_maxAA + 1), (double)_ktup);
- }
- tl1 = (l1 + l2) - 1;
-
- for (i = 1; i <= tl1; ++i)
- {
- slopes[i] = displ[i] = 0;
- diagIndex[i] = i;
- }
-
-
- // increment diagonal score for each k_tuple match
-
- for (i = 1; i <= limit; ++i)
- {
- vn1 = zzc[i];
- while (true)
- {
- if (!vn1)
- {
- break;
- }
- vn2 = zzd[i];
- while (vn2 != 0)
- {
- osptr = vn1 - vn2 + l2;
- ++displ[osptr];
- vn2 = zzb[vn2];
- }
- vn1 = zza[vn1];
- }
- }
-
- // choose the top SIGNIF diagonals
-
- desQuickSort(displ, diagIndex, tl1);
-
- j = tl1 - _signif + 1;
- if (j < 1)
- {
- j = 1;
- }
-
- // flag all diagonals within WINDOW of a top diagonal
-
- for (i = tl1; i >= j; i--)
- if (displ[i] > 0)
- {
- pos = diagIndex[i];
- l = (1 > pos - _window) ? 1 : pos - _window;
- m = (tl1 < pos + _window) ? tl1 : pos + _window;
- for (; l <= m; l++)
- {
- slopes[l] = 1;
- }
- }
-
- for (i = 1; i <= tl1; i++)
- {
- displ[i] = 0;
- }
-
-
- currFrag = maxSoFar = 0;
-
- for (i = 1; i <= (l1 - _ktup + 1); ++i)
- {
- encrypt = flag = 0;
- if (_DNAFlag){
- for (j = 1; j <= _ktup; ++j)
- {
- residue = ziAA2DNA[(*seq)[i + j - 1]];
- if ((residue < 0) || (residue > ciMaxResID))
- {
- flag = true;
- break;
- }
- encrypt += ((residue) * pot[j]);
- }
- }
- else {
- for (j = 1; j <= _ktup; ++j)
- {
- residue = (*seq)[i + j - 1];
- if ((residue < 0) || (residue > _maxAA))
- {
- flag = true;
- break;
- }
- encrypt += ((residue) * pot[j]);
- }
- }
-
- if (flag)
- {
- continue;
- }
- ++encrypt;
-
- vn2 = zzd[encrypt];
-
- flag = false;
- while (true)
- {
- if (!vn2)
- {
- flag = true;
- break;
- }
- osptr = i - vn2 + l2;
- if (slopes[osptr] != 1)
- {
- vn2 = zzb[vn2];
- continue;
- }
- flen = 0;
- fs = _ktup;
- next = maxSoFar;
-
-
- // A-loop
-
-
- while (true)
- {
- if (!next)
- {
- ++currFrag;
- if (currFrag >= 2 * _maxAlnLength)
- {
- utilityObject->info("(Partial alignment)");
- vatend = 1;
- return ;
- }
- displ[osptr] = currFrag;
- putFrag(fs, i, vn2, flen);
- }
- else
- {
- tv1 = accum[1][next];
- tv2 = accum[2][next];
- if (fragRelPos(i, vn2, tv1, tv2))
- {
- if (i - vn2 == accum[1][next] - accum[2][next])
- {
- if (i > accum[1][next] + (_ktup - 1))
- {
- fs = accum[0][next] + _ktup;
- }
- else
- {
- rmndr = i - accum[1][next];
- fs = accum[0][next] + rmndr;
- }
- flen = next;
- next = 0;
- continue;
- }
- else
- {
- if (displ[osptr] == 0)
- {
- subt1 = _ktup;
- }
- else
- {
- if (i > accum[1][displ[osptr]] + (_ktup - 1))
- {
- subt1 = accum[0][displ[osptr]] + _ktup;
- }
- else
- {
- rmndr = i - accum[1][displ[osptr]];
- subt1 = accum[0][displ[osptr]] + rmndr;
- }
- }
- subt2 = accum[0][next] - _windowGap + _ktup;
- if (subt2 > subt1)
- {
- flen = next;
- fs = subt2;
- }
- else
- {
- flen = displ[osptr];
- fs = subt1;
- }
- next = 0;
- continue;
- }
- }
- else
- {
- next = accum[4][next];
- continue;
- }
- }
- break;
- }
- // End of Aloop
- vn2 = zzb[vn2];
- }
- }
- vatend = 0;
-}
-
-void FastPairwiseAlign::makePPtrs(vector<int>& tptr, vector<int>& pl, const vector<int>* seq, int length)
-{
- int a[10];
- int i, j, limit, code;
- bool flag;
- int residue;
- int _ktup = userParameters->getKtup();
- int _maxAA = userParameters->getMaxAA();
-
- for (i = 1; i <= _ktup; i++)
- {
- a[i] = (int)pow((double)(_maxAA + 1), (double)(i - 1));
- }
-
- limit = (int)pow((double)(_maxAA + 1), (double)_ktup);
- if(limit >= (int)pl.size())
- {
- pl.resize(limit + 1);
- }
- if(length >= (int)tptr.size())
- {
- tptr.resize(length + 1);
- }
-
- for (i = 1; i <= limit; ++i)
- {
- pl[i] = 0;
- }
- for (i = 1; i <= length; ++i) // NOTE changed this
- {
- tptr[i] = 0;
- }
-
- for (i = 1; i <= (length - _ktup + 1); ++i)
- {
- code = 0;
- flag = false;
- for (j = 1; j <= _ktup; ++j)
- {
- residue = (*seq)[i + j - 1];
- if ((residue < 0) || (residue > _maxAA))
- {
- flag = true;
- break;
- }
- code += ((residue) * a[j]);
- }
- if (flag)
- {
- continue;
- }
- ++code;
- if (pl[code] != 0)
- {
- tptr[i] = pl[code];
- }
- pl[code] = i; //Ktuple code is at position i in sequence
- }
-}
-
-void FastPairwiseAlign::makeNPtrs(vector<int>& tptr,vector<int>& pl, const vector<int>* seq, int length)
-{
- int pot[] =
- {
- 0, 1, 4, 16, 64, 256, 1024, 4096
- };
- int i, j, limit, code;
- bool flag;
- int residue;
- int _ktup = userParameters->getKtup();
-
- limit = (int)pow((double)4, (double)_ktup);
-
- if(limit >= (int)pl.size())
- {
- pl.resize(limit + 1);
- }
- if(length >= (int)tptr.size())
- {
- tptr.resize(length + 1);
- }
-
- for (i = 1; i <= limit; ++i)
- {
- pl[i] = 0;
- }
- for (i = 1; i <= length; ++i)
- {
- tptr[i] = 0;
- }
-
- for (i = 1; i <= length - _ktup + 1; ++i)
- {
- code = 0;
- flag = false;
- for (j = 1; j <= _ktup; ++j)
- {
- residue = ziAA2DNA[(*seq)[i + j - 1]];
- if ((residue < 0) || (residue > ciMaxResID))
- {
- flag = true;
- break;
- }
- code += ((residue) * pot[j]);
- }
- if (flag)
- {
- continue;
- }
- ++code;
- if (pl[code] != 0)
- {
- tptr[i] = pl[code];
- }
- pl[code] = i;
- }
-
-}
-
-void FastPairwiseAlign::putFrag(int fs, int v1, int v2, int flen)
-{
- int end;
- accum[0][currFrag] = fs;
- accum[1][currFrag] = v1;
- accum[2][currFrag] = v2;
- accum[3][currFrag] = flen;
-
- if (!maxSoFar)
- {
- maxSoFar = 1;
- accum[4][currFrag] = 0;
- return ;
- }
-
- if (fs >= accum[0][maxSoFar])
- {
- accum[4][currFrag] = maxSoFar;
- maxSoFar = currFrag;
- return ;
- }
- else
- {
- next = maxSoFar;
- while (true)
- {
- end = next;
- next = accum[4][next];
- if (fs >= accum[0][next])
- {
- break;
- }
- }
- accum[4][currFrag] = next;
- accum[4][end] = currFrag;
- }
-}
-
-inline int FastPairwiseAlign::fragRelPos(int a1, int b1, int a2, int b2)
-{
- int ret;
- int _ktup = userParameters->getKtup();
-
- ret = false;
- if (a1 - b1 == a2 - b2)
- {
- if (a2 < a1)
- {
- ret = true;
- }
- }
- else
- {
- if (a2 + _ktup - 1 < a1 && b2 + _ktup - 1 < b1)
- {
- ret = true;
- }
- }
- return ret;
-}
-
-void FastPairwiseAlign::desQuickSort(vector<int>& array1, vector<int>& array2, int arraySize)
-{
- // Quicksort routine, adapted from chapter 4, page 115 of software tools
- // by Kernighan and Plauger, (1986)
- // Sort the elements of array1 and sort the
- // elements of array2 accordingly
- int temp1, temp2;
- int p, pivlin;
- int i, j;
- int lst[50], ust[50]; // the maximum no. of elements must be
- // < log(base2) of 50
-
- lst[1] = 1;
- ust[1] = arraySize - 1;
- p = 1;
-
- while (p > 0)
- {
- if (lst[p] >= ust[p])
- {
- p--;
- }
- else
- {
- i = lst[p] - 1;
- j = ust[p];
- pivlin = array1[j];
- while (i < j)
- {
- for (i = i + 1; array1[i] < pivlin; i++)
- ;
- for (j = j - 1; j > i; j--)
- if (array1[j] <= pivlin)
- {
- break;
- }
- if (i < j)
- {
- temp1 = array1[i];
- array1[i] = array1[j];
- array1[j] = temp1;
-
- temp2 = array2[i];
- array2[i] = array2[j];
- array2[j] = temp2;
- }
- }
-
- j = ust[p];
-
- temp1 = array1[i];
- array1[i] = array1[j];
- array1[j] = temp1;
-
- temp2 = array2[i];
- array2[i] = array2[j];
- array2[j] = temp2;
-
- if (i - lst[p] < ust[p] - i)
- {
- lst[p + 1] = lst[p];
- ust[p + 1] = i - 1;
- lst[p] = i + 1;
- }
- else
- {
- lst[p + 1] = i + 1;
- ust[p + 1] = ust[p];
- ust[p] = i - 1;
- }
- p = p + 1;
- }
- }
- return ;
-
-}
-
-}
git://source.jalview.org / jabaws.git / blobdiff