--- /dev/null
+/**
+ * Author: Mark Larkin
+ *
+ * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
+ */
+/**
+ * @author Mark Larkin, Conway Institute, UCD. mark.larkin@ucd.ie
+ * Changes:
+ *
+ * Mark: 23-01-2007: There was a problem with running an alignment with only
+ * one sequence in it. I needed to make a change to the multiSeqAlign function.
+ ****************************************************************************/
+#ifdef HAVE_CONFIG_H
+ #include "config.h"
+#endif
+#include "MSA.h"
+#include "MyersMillerProfileAlign.h"
+//#include "../phylogeneticTree/ClusterTree.h"
+#include "../general/debuglogObject.h"
+
+namespace clustalw
+{
+
+/**
+ *
+ * @param alnPtr
+ * @param distMat
+ * @param iStart
+ * @param phylipName
+ * @return
+ */
+int MSA::multiSeqAlign(Alignment* alnPtr, DistMatrix* distMat, vector<int>* seqWeight, AlignmentSteps* progSteps, int iStart)
+{
+
+ if(!progSteps)
+ {
+ return 0;
+ }
+
+ int* aligned;
+ vector<int> group;
+ int ix;
+
+ int* maxid;
+ int max = 0, sum = 0;
+ vector<int> treeWeight;
+ int i = 0, j = 0, set = 0, iseq = 0;
+ int entries = 0;
+ int score = 0;
+ int _numSteps = 0;
+ utilityObject->info("Start of Multiple Alignment\n");
+
+ int _numSeqs = alnPtr->getNumSeqs();
+
+ vector<int> newOutputIndex(_numSeqs);
+
+ alnPtr->addSeqWeight(seqWeight);
+
+ ProfileAlignAlgorithm* alignAlgorithm = new MyersMillerProfileAlign;
+ _numSteps = progSteps->getNumSteps();
+ // for each sequence, find the most closely related sequence
+
+ maxid = new int[_numSeqs + 1];
+
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ maxid[i] = -1;
+ for (j = 1; j <= _numSeqs; j++)
+ if (j != i && maxid[i] < (*distMat)(i, j))
+ {
+ maxid[i] = static_cast<int>((*distMat)(i, j));
+ }
+ }
+
+ // group the sequences according to their relative divergence
+
+ if (iStart == 0)
+ {
+
+ // start the multiple alignments.........
+
+ utilityObject->info("Aligning...");
+ // first pass, align closely related sequences first....
+
+
+ ix = 0;
+ aligned = new int[_numSeqs + 1];
+
+ for (i = 0; i <= _numSeqs; i++)
+ {
+ aligned[i] = 0;
+ }
+
+ const vector<vector<int> >* ptrToSets = progSteps->getSteps();
+
+
+ for (set = 1; set <= _numSteps; ++set)
+ {
+ entries = 0;
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ if (((*ptrToSets)[set][i] != 0) &&
+ (maxid[i] > userParameters->getDivergenceCutoff()))
+ {
+ entries++;
+ if (aligned[i] == 0)
+ {
+ if (userParameters->getOutputOrder() == INPUT)
+ {
+ ++ix;
+ newOutputIndex[i - 1] = i;
+ }
+ else
+ {
+ if(ix >= (int)newOutputIndex.size())
+ {
+ cerr << "ERROR: size = " << newOutputIndex.size()
+ << "ix = " << ix << "\n";
+ exit(1);
+ }
+ else
+ {
+ newOutputIndex[ix] = i;
+ ++ix;
+ }
+ }
+
+ aligned[i] = 1;
+ }
+ }
+ }
+
+ if (entries > 0)
+ {
+ #if DEBUGFULL
+ if(logObject && DEBUGLOG)
+ {
+ logObject->logMsg("Doing profile align");
+ }
+ #endif
+ score = alignAlgorithm->profileAlign(alnPtr, distMat, progSteps->getStep(set),
+ aligned);
+ }
+ else
+ {
+ score = 0;
+ }
+
+
+ // negative score means fatal error... exit now!
+
+ if (score < 0)
+ {
+ return (-1);
+ }
+ if(userParameters->getDisplayInfo())
+ {
+ if ((entries > 0) && (score > 0))
+ {
+ utilityObject->info("Group %d: Sequences:%4d Score:%d",
+ set, entries, score);
+ }
+ else
+ {
+ utilityObject->info("Group %d: Delayed", set);
+ }
+ }
+ }
+ }
+
+ else
+ {
+ aligned = new int[_numSeqs + 1];
+ ix = 0;
+ for (i = 1; i <= iStart + 1; i++)
+ {
+ aligned[i] = 1;
+ ++ix;
+ newOutputIndex[i - 1] = i;
+ }
+ for (i = iStart + 2; i <= _numSeqs; i++)
+ {
+ aligned[i] = 0;
+ }
+ }
+
+ // second pass - align remaining, more divergent sequences.....
+
+ // if not all sequences were aligned, for each unaligned sequence,
+ // find it's closest pair amongst the aligned sequences.
+
+ group.resize(_numSeqs + 1);
+ treeWeight.resize(_numSeqs);
+
+ for (i = 0; i < _numSeqs; i++)
+ {
+ treeWeight[i] = (*seqWeight)[i];
+ }
+
+ // if we haven't aligned any sequences, in the first pass - align the
+ // two most closely related sequences now
+ if (ix == 0)
+ {
+ max = -1;
+ iseq = 0;
+
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ for (j = i + 1; j <= _numSeqs; j++)
+ {
+ if (max < (*distMat)(i, j))
+ {
+ max = static_cast<int>((*distMat)(i, j)); // Mark change 17-5-07
+ iseq = i;
+ }
+ }
+ }
+ aligned[iseq] = 1;
+ if (userParameters->getOutputOrder() == INPUT)
+ {
+ ++ix;
+ newOutputIndex[iseq - 1] = iseq;
+ }
+ else
+ {
+ newOutputIndex[ix] = iseq;
+ ++ix;
+ }
+ }
+
+ while (ix < _numSeqs)
+ {
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ if (aligned[i] == 0)
+ {
+ maxid[i] = - 1;
+ for (j = 1; j <= _numSeqs; j++)
+ if ((maxid[i] < (*distMat)(i, j)) && (aligned[j] != 0))
+ {
+ maxid[i] = static_cast<int>((*distMat)(i, j));// Mark change 17-5-07
+ }
+
+ }
+ }
+ // find the most closely related sequence to those already aligned
+
+ max = - 1;
+ iseq = 0;
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ if ((aligned[i] == 0) && (maxid[i] > max))
+ {
+ max = maxid[i];
+ iseq = i;
+ }
+ }
+
+
+ // align this sequence to the existing alignment
+ // weight sequences with percent identity with profile
+ // OR...., multiply sequence weights from tree by percent identity with new sequence
+ if (userParameters->getNoWeights() == false)
+ {
+ for (j = 0; j < _numSeqs; j++)
+ if (aligned[j + 1] != 0)
+ {
+ (*seqWeight)[j] = static_cast<int>(treeWeight[j] * (*distMat)(j + 1, iseq));
+ }
+
+ // Normalise the weights, such that the sum of the weights = INT_SCALE_FACTOR
+
+ sum = 0;
+ for (j = 0; j < _numSeqs; j++)
+ {
+ if (aligned[j + 1] != 0)
+ {
+ sum += (*seqWeight)[j];
+ }
+ }
+
+ if (sum == 0)
+ {
+ for (j = 0; j < _numSeqs; j++)
+ {
+ (*seqWeight)[j] = 1;
+ }
+
+ sum = j;
+ }
+ for (j = 0; j < _numSeqs; j++)
+ {
+ if (aligned[j + 1] != 0)
+ {
+ (*seqWeight)[j] = ((*seqWeight)[j] * INT_SCALE_FACTOR) / sum;
+ if ((*seqWeight)[j] < 1)
+ {
+ (*seqWeight)[j] = 1;
+ }
+ }
+ }
+ }
+
+ entries = 0;
+ for (j = 1; j <= _numSeqs; j++)
+ {
+ if (aligned[j] != 0)
+ {
+ group[j] = 1;
+ entries++;
+ }
+ else if (iseq == j)
+ {
+ group[j] = 2;
+ entries++;
+ }
+ }
+
+ alnPtr->addSeqWeight(seqWeight);
+ aligned[iseq] = 1;
+
+ score = alignAlgorithm->profileAlign(alnPtr, distMat, &group, aligned);
+
+ if (userParameters->getOutputOrder() == INPUT)
+ {
+ ++ix;
+ newOutputIndex[iseq - 1] = iseq;
+ }
+ else
+ {
+ newOutputIndex[ix] = iseq;
+ ++ix;
+ }
+ }
+
+ alnPtr->addOutputIndex(&newOutputIndex);
+
+ if(userParameters->getDisplayInfo())
+ {
+ int alignmentScore = alnPtr->alignScore(); // ?? check, FS, 2009-05-18
+ }
+
+ delete alignAlgorithm;
+ delete [] aligned;
+ delete [] maxid;
+ return (_numSeqs);
+}
+
+
+/**
+ *
+ * @param alnPtr
+ * @param distMat
+ * @param iStart
+ * @param phylipName
+ * @return
+ */
+int MSA::seqsAlignToProfile(Alignment* alnPtr, DistMatrix* distMat, vector<int>* seqWeight, int iStart,
+ string phylipName)
+{
+ int *aligned;
+ vector<int> treeWeight;
+ vector<int> group;
+ int ix;
+
+ int *maxid;
+ int max = 0;
+ int i = 0, j = 0, iseq = 0;
+ int sum = 0, entries = 0;
+ int score = 0;
+ int _numSeqs = alnPtr->getNumSeqs();
+
+ utilityObject->info("Start of Multiple Alignment\n");
+
+ ProfileAlignAlgorithm* alignAlgorithm = new MyersMillerProfileAlign;
+
+ // calculate sequence weights according to branch lengths of the tree -
+ // weights in global variable seq_weight normalised to sum to 100
+ vector<int> newOutputIndex(_numSeqs);
+ //groupTree.calcSeqWeights(0, _numSeqs, &seqWeight);
+
+ treeWeight.resize(_numSeqs);
+ for (i = 0; i < _numSeqs; i++)
+ {
+ treeWeight[i] = (*seqWeight)[i];
+ }
+
+ // for each sequence, find the most closely related sequence
+
+ maxid = new int[_numSeqs + 1];
+
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ maxid[i] = - 1;
+ for (j = 1; j <= _numSeqs; j++)
+ {
+ if (maxid[i] < (*distMat)(i, j))
+ {
+ maxid[i] = static_cast<int>((*distMat)(i, j)); // Mark change 17-5-07
+ }
+ }
+ }
+
+ aligned = new int[_numSeqs + 1];
+ ix = 0;
+ for (i = 1; i <= iStart + 1; i++)
+ {
+ aligned[i] = 1;
+ ++ix;
+ newOutputIndex[i - 1] = i;
+ }
+ for (i = iStart + 2; i <= _numSeqs; i++)
+ {
+ aligned[i] = 0;
+ }
+
+ // for each unaligned sequence, find it's closest pair amongst the
+ // aligned sequences.
+
+ group.resize(_numSeqs + 1);
+
+ while (ix < _numSeqs)
+ {
+ if (ix > 0)
+ {
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ if (aligned[i] == 0)
+ {
+ maxid[i] = - 1;
+ for (j = 1; j <= _numSeqs; j++)
+ {
+ if ((maxid[i] < (*distMat)(i, j)) && (aligned[j] != 0))
+ {
+ maxid[i] = static_cast<int>((*distMat)(i, j));
+ }
+ }
+ }
+ }
+ }
+
+ // find the most closely related sequence to those already aligned
+
+ max = -1;
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ if ((aligned[i] == 0) && (maxid[i] > max))
+ {
+ max = maxid[i];
+ iseq = i;
+ }
+ }
+
+ // align this sequence to the existing alignment
+ entries = 0;
+ for (j = 1; j <= _numSeqs; j++)
+ {
+ if (aligned[j] != 0)
+ {
+ group[j] = 1;
+ entries++;
+ }
+ else if (iseq == j)
+ {
+ group[j] = 2;
+ entries++;
+ }
+ }
+
+ aligned[iseq] = 1;
+
+ // multiply sequence weights from tree by percent
+ // identity with new sequence
+
+ for (j = 0; j < _numSeqs; j++)
+ {
+ (*seqWeight)[j] = static_cast<int>(treeWeight[j] * (*distMat)(j + 1, iseq));
+ }
+
+ //
+ // Normalise the weights, such that the sum of the weights = INT_SCALE_FACTOR
+ //
+
+ sum = 0;
+ for (j = 0; j < _numSeqs; j++)
+ {
+ if (group[j + 1] == 1)
+ {
+ sum += (*seqWeight)[j];
+ }
+ }
+
+ if (sum == 0)
+ {
+ for (j = 0; j < _numSeqs; j++)
+ {
+ (*seqWeight)[j] = 1;
+ }
+
+ sum = j;
+ }
+ for (j = 0; j < _numSeqs; j++)
+ {
+ (*seqWeight)[j] = ((*seqWeight)[j] * INT_SCALE_FACTOR) / sum;
+ if ((*seqWeight)[j] < 1)
+ {
+ (*seqWeight)[j] = 1;
+ }
+ }
+ // Add the seqWeights to the Alignment object!!!!
+ alnPtr->addSeqWeight(seqWeight);
+
+ score = alignAlgorithm->profileAlign(alnPtr, distMat, &group, aligned);
+ utilityObject->info("Sequence:%d Score:%d", iseq, score);
+ if (userParameters->getOutputOrder() == INPUT)
+ {
+ ++ix;
+ newOutputIndex[iseq - 1] = iseq;
+ }
+ else
+ {
+ newOutputIndex[ix] = iseq;
+ ++ix;
+ }
+ }
+
+ delete [] aligned;
+ delete [] maxid;
+ delete alignAlgorithm;
+ alnPtr->addOutputIndex(&newOutputIndex);
+
+ if(userParameters->getDisplayInfo())
+ {
+ alnPtr->alignScore();
+ }
+
+ return (_numSeqs);
+}
+
+
+/**
+ *
+ * @param alnPtr
+ * @param distMat
+ * @return
+ */
+int MSA::calcPairwiseForProfileAlign(Alignment* alnPtr, DistMatrix* distMat)
+{
+ //Tree groupTree;
+ int i, j, temp;
+ int entries;
+ int* aligned;
+ vector<int> group;
+ vector<int> seqWeight;
+ float dscore;
+ int score;
+ int _numSeqs = alnPtr->getNumSeqs();
+
+ seqWeight.resize(_numSeqs);
+ ProfileAlignAlgorithm* alignAlg = new MyersMillerProfileAlign;
+
+ utilityObject->info("Start of Initial Alignment");
+ /* calculate sequence weights according to branch lengths of the tree -
+ * weights in global variable seq_weight normalised to sum to INT_SCALE_FACTOR */
+
+ temp = INT_SCALE_FACTOR / _numSeqs;
+ for (i = 0; i < _numSeqs; i++)
+ {
+ seqWeight[i] = temp;
+ }
+
+ userParameters->setDistanceTree(false);
+
+ /* do the initial alignment......... */
+
+ group.resize(_numSeqs + 1);
+
+ for (i = 1; i <= alnPtr->getProfile1NumSeqs(); ++i)
+ {
+ group[i] = 1;
+ }
+
+ for (i = alnPtr->getProfile1NumSeqs() + 1; i <= _numSeqs; ++i)
+ {
+ group[i] = 2;
+ }
+
+ entries = _numSeqs;
+
+ aligned = new int[_numSeqs + 1];
+
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ aligned[i] = 1;
+ }
+
+ alnPtr->addSeqWeight(&seqWeight);
+
+ score = alignAlg->profileAlign(alnPtr, distMat, &group, aligned);
+ utilityObject->info("Sequences:%d Score:%d", entries, score);
+ delete [] aligned;
+
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ for (j = i + 1; j <= _numSeqs; j++)
+ {
+ dscore = alnPtr->countid(i, j);
+ (*distMat)(i, j) = ((double)100.0 - (double)dscore) / (double)100.0;
+ (*distMat)(j, i) = (*distMat)(i, j);
+ }
+ }
+ delete alignAlg;
+ return (_numSeqs);
+
+}
+
+
+/**
+ *
+ * @param alnPtr
+ * @param distMat
+ * @param p1TreeName
+ * @param p2TreeName
+ * @return
+ */
+int MSA::doProfileAlign(Alignment* alnPtr, DistMatrix* distMat, vector<int>* prof1Weight, vector<int>* prof2Weight)
+{
+ //Tree groupTree1, groupTree2;
+ int i, j, sum, entries;
+ int score;
+ int *aligned;
+ vector<int> group;
+ int *maxid;
+ //vector<int> prof1Weight, prof2Weight;
+ int _profile1NumSeqs = alnPtr->getProfile1NumSeqs();
+ int _numSeqs = alnPtr->getNumSeqs();
+ vector<int> _seqWeight, _outputIndex;
+
+ utilityObject->info("Start of Multiple Alignment\n");
+
+ _seqWeight.resize(_numSeqs + 1);
+ _outputIndex.resize(_numSeqs);
+ ProfileAlignAlgorithm* alignAlgorithm = new MyersMillerProfileAlign;
+
+ // weight sequences with max percent identity with other profile
+
+ maxid = new int[_numSeqs + 1];
+ for (i = 0; i < _profile1NumSeqs; i++)
+ {
+ maxid[i] = 0;
+ for (j = _profile1NumSeqs + 1; j <= _numSeqs; j++)
+ {
+ if (maxid[i] < (*distMat)(i + 1, j))
+ {
+ maxid[i] = static_cast<int>((*distMat)(i + 1, j)); // Mark change 17-5-07
+ }
+ }
+ _seqWeight[i] = maxid[i] * (*prof1Weight)[i];
+ }
+
+ for (i = _profile1NumSeqs; i < _numSeqs; i++)
+ {
+ maxid[i] = - 1;
+ for (j = 1; j <= _profile1NumSeqs; j++)
+ {
+ if (maxid[i] < (*distMat)(i + 1, j))
+ {
+ maxid[i] = static_cast<int>((*distMat)(i + 1, j));// Mark change 17-5-07
+ }
+ }
+ _seqWeight[i] = maxid[i] * (*prof2Weight)[i];
+ }
+ //
+ // Normalise the weights, such that the sum of the weights = INT_SCALE_FACTOR
+ //
+
+ sum = 0;
+ for (j = 0; j < _numSeqs; j++)
+ {
+ sum += _seqWeight[j];
+ }
+
+ if (sum == 0)
+ {
+ for (j = 0; j < _numSeqs; j++)
+ {
+ _seqWeight[j] = 1;
+ }
+
+ sum = j;
+ }
+ for (j = 0; j < _numSeqs; j++)
+ {
+ _seqWeight[j] = (_seqWeight[j] * INT_SCALE_FACTOR) / sum;
+ if (_seqWeight[j] < 1)
+ {
+ _seqWeight[j] = 1;
+ }
+ }
+
+ // do the alignment......... /
+
+ utilityObject->info("Aligning...");
+ group.resize(_numSeqs + 1);
+
+ for (i = 1; i <= _profile1NumSeqs; ++i)
+ {
+ group[i] = 1;
+ }
+
+ for (i = _profile1NumSeqs + 1; i <= _numSeqs; ++i)
+ {
+ group[i] = 2;
+ }
+
+ entries = _numSeqs;
+
+ aligned = new int[_numSeqs + 1];
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ aligned[i] = 1;
+ }
+
+ alnPtr->addSeqWeight(&_seqWeight);
+
+ score = alignAlgorithm->profileAlign(alnPtr, distMat, &group, aligned);
+
+ utilityObject->info("Sequences:%d Score:%d", entries, score);
+
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ _outputIndex[i - 1] = i;
+ }
+
+ alnPtr->addOutputIndex(&_outputIndex);
+
+ delete alignAlgorithm;
+ delete [] aligned;
+ delete [] maxid;
+ return (_numSeqs);
+ return 1;
+}
+
+}
git://source.jalview.org / jabaws.git / blobdiff