--- /dev/null
+/**
+ * Author: Mark Larkin
+ *
+ * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
+ */
+#ifdef HAVE_CONFIG_H
+ #include "config.h"
+#endif
+#include "ClusterTree.h"
+#include "../general/utils.h"
+#include "dayhoff.h"
+#include "RandomGenerator.h"
+#include <math.h>
+#include <sstream>
+#include "../general/OutputFile.h"
+
+namespace clustalw
+{
+
+ClusterTree::ClusterTree()
+ : numSeqs(0),
+ firstSeq(0),
+ lastSeq(0)
+{
+ bootstrapPrompt = "\nEnter name for bootstrap output file ";
+ bootstrapFileTypeMsg = "Bootstrap output";
+}
+
+
+
+/** ****************************************************************************************
+ * The rest are private functions. *
+ * *
+ * *
+ * *
+ *******************************************************************************************/
+
+
+
+void ClusterTree::distanceMatrixOutput(ofstream* outFile, clustalw::DistMatrix* matToPrint,
+ clustalw::Alignment *alignPtr)
+{
+ if(outFile == NULL || !outFile->is_open())
+ {
+ clustalw::utilityObject->error("Cannot output the distance matrix, file is not open\n");
+ return;
+ }
+
+ int i, j;
+ int _maxNames = alignPtr->getMaxNames();
+ (*outFile) << setw(6) << lastSeq - firstSeq + 1;
+
+ for (i = 1; i <= lastSeq - firstSeq + 1; i++)
+ {
+ (*outFile) << "\n" << left << setw(_maxNames) << alignPtr->getName(i) << " ";
+ for (j = 1; j <= lastSeq - firstSeq + 1; j++)
+ {
+ (*outFile) << " " << setw(6) << setprecision(3) << fixed << (*matToPrint)(i, j);
+ if (j % 8 == 0)
+ {
+ if (j != lastSeq - firstSeq + 1)
+ {
+ (*outFile) << "\n";
+ }
+ if (j != lastSeq - firstSeq + 1)
+ {
+ (*outFile) << " ";
+ }
+ }
+ }
+ }
+}
+
+void ClusterTree::overspillMessage(int overspill,int totalDists)
+{
+ std::stringstream ssOverSpill;
+ ssOverSpill << overspill << " of the distances out of a total of " << totalDists
+ << "\n were out of range for the distance correction." << "\n"
+ << "\n SUGGESTIONS: 1) remove the most distant sequences"
+ << "\n or 2) use the PHYLIP package"
+ << "\n or 3) turn off the correction."
+ << "\n Note: Use option 3 with caution! With this degree"
+ << "\n of divergence you will have great difficulty"
+ << "\n getting robust and reliable trees." << "\n\n";
+ string message;
+ ssOverSpill >> message;
+ clustalw::utilityObject->warning(message.c_str());
+}
+
+/*
+ * The function treeGapDelete flags all the positions that have a gap in any sequence.
+ *
+ */
+// NOTE there is something wrong with using _lenFirstSeq. But this is what old clustal does.
+void ClusterTree::treeGapDelete(clustalw::Alignment *alignPtr)
+{
+ int seqn;
+ int posn;
+ int _maxAlnLength = alignPtr->getMaxAlnLength();
+ int _lenFirstSeq = alignPtr->getSeqLength(firstSeq);
+ int _gapPos1 = clustalw::userParameters->getGapPos1();
+ int _gapPos2 = clustalw::userParameters->getGapPos2();
+
+ treeGaps.resize(_maxAlnLength + 1);
+
+ for (posn = 1; posn <= _lenFirstSeq; ++posn)
+ {
+ treeGaps[posn] = 0;
+ for (seqn = 1; seqn <= lastSeq - firstSeq + 1; ++seqn)
+ {
+ const vector<int>* _seqM = alignPtr->getSequence(seqn + firstSeq - 1);
+
+ if(posn > alignPtr->getSeqLength(seqn + firstSeq - 1))
+ {
+ break; // Dont read locations that cannot be read!
+ }
+ if (((*_seqM)[posn] == _gapPos1) || ((*_seqM)[posn] == _gapPos2))
+ {
+ treeGaps[posn] = 1;
+ break;
+ }
+ }
+ }
+}
+
+int ClusterTree::dnaDistanceMatrix(ofstream* treeFile, clustalw::Alignment *alignPtr)
+{
+ int m, n;
+ int j, i;
+ int res1, res2;
+ int overspill = 0;
+ double p, q, e, a, b, k;
+
+ treeGapDelete(alignPtr); // flag positions with gaps (tree_gaps[i] = 1 )
+
+ if (verbose)
+ {
+ (*treeFile) << "\n";
+ (*treeFile) << "\n DIST = percentage divergence (/100)";
+ (*treeFile) << "\n p = rate of transition (A <-> G; C <-> T)";
+ (*treeFile) << "\n q = rate of transversion";
+ (*treeFile) << "\n Length = number of sites used in comparison";
+ (*treeFile) << "\n";
+ if (clustalw::userParameters->getTossGaps())
+ {
+ (*treeFile) << "\n All sites with gaps (in any sequence) deleted!";
+ (*treeFile) << "\n";
+ }
+ if (clustalw::userParameters->getKimura())
+ {
+ (*treeFile) << "\n Distances corrected by Kimura's 2 parameter model:";
+ (*treeFile) << "\n\n Kimura, M. (1980)";
+ (*treeFile) << " A simple method for estimating evolutionary ";
+ (*treeFile) << "rates of base";
+ (*treeFile) << "\n substitutions through comparative studies of ";
+ (*treeFile) << "nucleotide sequences.";
+ (*treeFile) << "\n J. Mol. Evol., 16, 111-120.";
+ (*treeFile) << "\n\n";
+ }
+ }
+
+ int _numSeqs = alignPtr->getNumSeqs();
+ quickDistMat.reset(new clustalw::DistMatrix(_numSeqs + 1));
+ int _lenFirstSeq = alignPtr->getSeqLength(firstSeq);
+ int _gapPos1 = clustalw::userParameters->getGapPos1();
+ int _gapPos2 = clustalw::userParameters->getGapPos2();
+ int lenSeqM, lenSeqN;
+ // for every pair of sequence
+ for (m = 1; m < lastSeq - firstSeq + 1; ++m)
+ {
+ const vector<int>* _seqM = alignPtr->getSequence(m + firstSeq - 1);
+ lenSeqM = alignPtr->getSeqLength(m + firstSeq - 1);
+
+ for (n = m + 1; n <= lastSeq - firstSeq + 1; ++n)
+ {
+ const vector<int>* _seqN = alignPtr->getSequence(n + firstSeq - 1);
+ lenSeqN = alignPtr->getSeqLength(n + firstSeq - 1);
+
+ p = q = e = 0.0;
+ quickDistMat->SetAt(m, n, 0.0);
+ quickDistMat->SetAt(n ,m, 0.0);
+
+ for (i = 1; i <= _lenFirstSeq; ++i)
+ {
+ j = bootPositions[i];
+ if (clustalw::userParameters->getTossGaps() && (treeGaps[j] > 0))
+ {
+ goto skip;
+ }
+ // gap position
+/** *******************************************************************************
+ * BUG!!!!!!! NOTE this was found for protein. Presuming the same here *
+ * NOTE: the following if statements were coded in so as to produce *
+ * the same distance results as the old clustal. Old clustal compares *
+ * up to the length of the first sequence. If this is longer than the *
+ * other sequences, then the -3 and 0's are compared at the end of the *
+ * array. These should not be compared, but I need to stick to this to *
+ * produce the same results as the old version for testing! *
+ **********************************************************************************/
+ if(j > lenSeqM)
+ {
+ if(j == lenSeqM + 1)
+ {
+ res1 = -3;
+ }
+ else
+ {
+ res1 = 0;
+ }
+ }
+ else
+ {
+ res1 = (*_seqM)[j];
+ }
+ if(j > lenSeqN)
+ {
+ if(j == lenSeqN + 1)
+ {
+ res2 = -3;
+ }
+ else
+ {
+ res2 = 0;
+ }
+ }
+ else
+ {
+ res2 = (*_seqN)[j];
+ }
+ if ((res1 == _gapPos1) || (res1 == _gapPos2) || (res2 == _gapPos1)
+ || (res2 == _gapPos2))
+ {
+ goto skip;
+ }
+ // gap in a seq
+ if (!clustalw::userParameters->getUseAmbiguities())
+ {
+ if (isAmbiguity(res1) || isAmbiguity(res2))
+ {
+ goto skip;
+ }
+ }
+ // ambiguity code in a seq
+ e = e+1.0;
+ if (res1 != res2)
+ {
+ if (transition(res1, res2))
+ {
+ p = p + 1.0;
+ }
+ else
+ {
+ q = q + 1.0;
+ }
+ }
+ skip: ;
+ }
+
+
+ // Kimura's 2 parameter correction for multiple substitutions
+
+ if (!clustalw::userParameters->getKimura())
+ {
+ if (e == 0)
+ {
+ cerr << "\n WARNING: sequences " << m << " and " << n
+ << " are non-overlapping\n";
+ k = 0.0;
+ p = 0.0;
+ q = 0.0;
+ }
+ else
+ {
+ k = (p + q) / e;
+ if (p > 0.0)
+ {
+ p = p / e;
+ }
+ else
+ {
+ p = 0.0;
+ }
+ if (q > 0.0)
+ {
+ q = q / e;
+ }
+ else
+ {
+ q = 0.0;
+ }
+ }
+ quickDistMat->SetAt(m, n, k);
+ quickDistMat->SetAt(n ,m, k);
+ if (verbose)
+ {
+ (*treeFile) << setw(4) << m << " vs." << setw(4) << n << ": DIST = "
+ << setw(4) << fixed << setprecision(4)
+ << k << "; p = " << fixed << setprecision(4) << p << "; q = "
+ << fixed << setprecision(4) << q << "; length = " << setw(6)
+ << fixed << setprecision(0) << e << "\n";
+ }
+ }
+ else
+ {
+ if (e == 0)
+ {
+ cerr << "\n WARNING: sequences " << m << " and " << n
+ << " are non-overlapping\n";;
+ p = 0.0;
+ q = 0.0;
+ }
+ else
+ {
+ if (p > 0.0)
+ {
+ p = p / e;
+ }
+ else
+ {
+ p = 0.0;
+ }
+ if (q > 0.0)
+ {
+ q = q / e;
+ }
+ else
+ {
+ q = 0.0;
+ }
+ }
+
+ if (((2.0 *p) + q) == 1.0)
+ {
+ a = 0.0;
+ }
+ else
+ {
+ a = 1.0 / (1.0 - (2.0 *p) - q);
+ }
+
+ if (q == 0.5)
+ {
+ b = 0.0;
+ }
+ else
+ {
+ b = 1.0 / (1.0 - (2.0 *q));
+ }
+
+ // watch for values going off the scale for the correction.
+ if ((a <= 0.0) || (b <= 0.0))
+ {
+ overspill++;
+ k = 3.5; // arbitrary high score
+ }
+ else
+ {
+ k = 0.5 * log(a) + 0.25 * log(b);
+ }
+ quickDistMat->SetAt(m, n, k);
+ quickDistMat->SetAt(n ,m, k);
+ if (verbose)
+ // if screen output
+ {
+ (*treeFile) << setw(4) << m << " vs." << setw(4) << n
+ << ": DIST = " << fixed << setprecision(4)
+ << k << "; p = " << fixed << setprecision(4) << p << "; q = "
+ << fixed << setprecision(4) << q << "; length = " << setw(6)
+ << fixed << setprecision(0) << e << "\n";
+ }
+
+ }
+ }
+ }
+ return overspill; // return the number of off-scale values
+}
+
+int ClusterTree::protDistanceMatrix(ofstream* treeFile, clustalw::Alignment *alignPtr)
+{
+ int m, n;
+ int j, i;
+ int res1, res2;
+ int overspill = 0;
+ double p, e, k, tableEntry;
+
+ treeGapDelete(alignPtr); // flag positions with gaps (tree_gaps[i] = 1 )
+
+ if (verbose)
+ {
+ (*treeFile) << "\n";
+ (*treeFile) << "\n DIST = percentage divergence (/100)";
+ (*treeFile) << "\n Length = number of sites used in comparison";
+ (*treeFile) << "\n\n";
+ if (clustalw::userParameters->getTossGaps())
+ {
+ (*treeFile) << "\n All sites with gaps (in any sequence) deleted";
+ (*treeFile) << "\n";
+ }
+ if (clustalw::userParameters->getKimura())
+ {
+ (*treeFile) <<
+ "\n Distances up to 0.75 corrected by Kimura's empirical method:";
+ (*treeFile) << "\n\n Kimura, M. (1983)";
+ (*treeFile) << " The Neutral Theory of Molecular Evolution.";
+ (*treeFile) <<
+ "\n Page 75. Cambridge University Press, Cambridge, England.";
+ (*treeFile) << "\n\n";
+ }
+ }
+ int _numSeqs = alignPtr->getNumSeqs();
+ int _lenSeq1 = alignPtr->getSeqLength(1);
+ quickDistMat.reset(new clustalw::DistMatrix(_numSeqs + 1));
+ int _gapPos1 = clustalw::userParameters->getGapPos1();
+ int _gapPos2 = clustalw::userParameters->getGapPos2();
+ int lenSeqM, lenSeqN;
+ // for every pair of sequence
+ for (m = 1; m < _numSeqs; ++m)
+ {
+ const vector<int>* _seqM = alignPtr->getSequence(m);
+ lenSeqM = alignPtr->getSeqLength(m);
+ for (n = m + 1; n <= _numSeqs; ++n)
+ {
+ const vector<int>* _seqN = alignPtr->getSequence(n);
+ lenSeqN = alignPtr->getSeqLength(n);
+ p = e = 0.0;
+ quickDistMat->SetAt(m, n, 0.0);
+ quickDistMat->SetAt(n ,m, 0.0);
+ for (i = 1; i <= _lenSeq1; ++i) // It may be this here!
+ {
+ j = bootPositions[i];
+ if (clustalw::userParameters->getTossGaps() && (treeGaps[j] > 0))
+ {
+ goto skip;
+ }
+ // gap position
+/** *******************************************************************************
+ * BUG!!!!!!! *
+ * NOTE: the following if statements were coded in so as to produce *
+ * the same distance results as the old clustal. Old clustal compares *
+ * up to the length of the first sequence. If this is longer than the *
+ * other sequences, then the -3 and 0's are compared at the end of the *
+ * array. These should not be compared, but I need to stick to this to *
+ * produce the same results as the old version for testing! *
+ **********************************************************************************/
+ if(j > lenSeqM)
+ {
+ if(j == lenSeqM + 1)
+ {
+ res1 = -3;
+ }
+ else
+ {
+ res1 = 0;
+ }
+ }
+ else
+ {
+ res1 = (*_seqM)[j];
+ }
+ if(j > lenSeqN)
+ {
+ if(j == lenSeqN + 1)
+ {
+ res2 = -3;
+ }
+ else
+ {
+ res2 = 0;
+ }
+ }
+ else
+ {
+ res2 = (*_seqN)[j];
+ }
+ if ((res1 == _gapPos1) || (res1 == _gapPos2) || (res2 == _gapPos1)
+ || (res2 == _gapPos2))
+ {
+ goto skip;
+ }
+ // gap in a seq
+ e = e + 1.0;
+ if (res1 != res2)
+ {
+ p = p + 1.0;
+ }
+ skip: ;
+ }
+
+ if (p <= 0.0)
+ {
+ k = 0.0;
+ }
+ else
+ {
+ k = p / e;
+ }
+
+ if (clustalw::userParameters->getKimura())
+ {
+ if (k < 0.75)
+ {
+ // use Kimura's formula
+ if (k > 0.0)
+ {
+ k = - log(1.0 - k - (k * k / 5.0));
+ }
+ }
+ else
+ {
+ if (k > 0.930)
+ {
+ overspill++;
+ k = 10.0; // arbitrarily set to 1000%
+ }
+ else // dayhoff_pams is from dayhoff.h file
+ {
+ tableEntry = (k * 1000.0) - 750.0;
+ k = (double)dayhoff_pams[(int)tableEntry];
+ k = k / 100.0;
+ }
+ }
+ }
+
+ quickDistMat->SetAt(m, n, k);
+ quickDistMat->SetAt(n ,m, k);
+ if (verbose)
+ {
+ (*treeFile) << setw(4) << m << " vs." << setw(4) << n
+ << " DIST = " << fixed << setprecision(4)
+ << k << "; length = " << setw(6)
+ << setprecision(0) << e << "\n";
+ }
+ }
+ }
+ return overspill;
+}
+
+bool ClusterTree::isAmbiguity(int c)
+{
+ int i;
+ char codes[] = "ACGTU";
+
+ if (clustalw::userParameters->getUseAmbiguities() == true)
+ {
+ return false;
+ }
+
+ for (i = 0; i < 5; i++)
+ if (clustalw::userParameters->getAminoAcidCode(c) == codes[i])
+ {
+ return false;
+ }
+
+ return true;
+}
+
+/*
+ * The function calcPercIdentity calculates the percent identity of the sequences
+ * and outputs it to a the file pfile. NOTE this is not used at the moment. It was in
+ * the old code, but there was no way to access it from the menu. This may change.
+ */
+void ClusterTree::calcPercIdentity(ofstream* pfile, clustalw::Alignment *alignPtr)
+{
+ clustalw::DistMatrix percentMat;
+
+ float ident;
+ int nmatch;
+
+ int val1, val2;
+
+ int i,j,k, length_longest;
+ int length_shortest;
+
+ int rs = 0, rl = 0;
+ // findout sequence length, longest and shortest;
+ length_longest = 0;
+ length_shortest = 0;
+
+ int _numSeqs = alignPtr->getNumSeqs();
+ int _seqLength;
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ _seqLength = alignPtr->getSeqLength(i);
+ if (length_longest < _seqLength)
+ {
+ length_longest = _seqLength;
+ rs = i;
+ }
+ if (length_shortest > _seqLength)
+ {
+ length_shortest = _seqLength;
+ rl = i;
+ }
+ }
+
+ percentMat.ResizeRect(_numSeqs + 1);
+ nmatch = 0;
+ int _lenSeqI, _lenSeqJ;
+ int _maxAA = clustalw::userParameters->getMaxAA();
+
+ for (i = 1; i <= _numSeqs; i++)
+ {
+ const vector<int>* _seqI = alignPtr->getSequence(i);
+ _lenSeqI = alignPtr->getSeqLength(i);
+
+ for (j = i; j <= _numSeqs; j++)
+ {
+ const vector<int>* _seqJ = alignPtr->getSequence(j);
+ _lenSeqJ = alignPtr->getSeqLength(j);
+
+ cout << "\n " << alignPtr->getName(j) << " ";
+ ident = 0;
+ nmatch = 0;
+ for(k = 1; k <= length_longest; k++)
+ {
+ if((k > _lenSeqI) || (k > _lenSeqJ))
+ {
+ break;
+ }
+ val1 = (*_seqI)[k];
+ val2 = (*_seqJ)[k];
+
+ if (((val1 < 0) || (val1 > _maxAA)) || ((val2 < 0) || (val2 > _maxAA)))
+ {
+ continue; // residue = '-';
+ }
+ if (val1 == val2)
+ {
+ ident++ ;
+ nmatch++;
+ }
+ else
+ {
+ nmatch++ ;
+ }
+ }
+ ident = ident/nmatch * 100.0 ;
+ percentMat.SetAt(i, j, ident);
+ percentMat.SetAt(j, i, ident);
+ }
+ }
+
+ int _maxNameSize = alignPtr->getMaxNames();
+
+ (*pfile) << "#\n#\n# Percent Identity Matrix - created by Clustal"
+ << clustalw::userParameters->getRevisionLevel() << " \n#\n#\n";
+ for(i = 1; i <= _numSeqs; i++)
+ {
+ (*pfile) << "\n " << right << setw(5) << i << ": ";
+ (*pfile) << left << setw(_maxNameSize) << alignPtr->getName(i);
+
+ for(j = 1; j <= _numSeqs; j++)
+ {
+ (*pfile) << setw(8) << right << fixed << setprecision(0) << percentMat(i, j);
+ }
+ }
+ (*pfile) << "\n";
+
+}
+
+void ClusterTree::compareTree(PhyloTree* tree1, PhyloTree* tree2, vector<int>* hits, int n)
+{
+ int i, j, k;
+ int nhits1, nhits2;
+
+ for (i = 1; i <= n - 3; i++)
+ {
+ for (j = 1; j <= n - 3; j++)
+ {
+ nhits1 = 0;
+ nhits2 = 0;
+ for (k = 1; k <= n; k++)
+ {
+ if (tree1->treeDesc[i][k] == tree2->treeDesc[j][k])
+ {
+ nhits1++;
+ }
+ if (tree1->treeDesc[i][k] != tree2->treeDesc[j][k])
+ {
+ nhits2++;
+ }
+ }
+ if ((nhits1 == lastSeq - firstSeq + 1) || (nhits2 == lastSeq -
+ firstSeq + 1))
+ {
+ (*hits)[i]++;
+ }
+ }
+ }
+}
+
+/**
+ * NOTE this will go into the OutputFile class and will not be needed here anymore.
+ *
+ */
+/*string ClusterTree::getOutputFileName(const string prompt, string path,
+ const string fileExtension)
+{
+ string temp;
+ string _fileName; // Will return this name.
+ string message;
+ _fileName = path + fileExtension;
+
+ if(_fileName.compare(clustalw::userParameters->getSeqName()) == 0)
+ {
+ cout << "Output file name is the same as input file.\n";
+ if (clustalw::userParameters->getMenuFlag())
+ {
+ message = "\n\nEnter new name to avoid overwriting [" + _fileName + "]: ";
+ clustalw::utilityObject->getStr(message, temp);
+ if(temp != "")
+ {
+ _fileName = temp;
+ }
+ }
+ }
+ else if (clustalw::userParameters->getMenuFlag())
+ {
+
+ message = prompt + " [" + _fileName + "]";
+ clustalw::utilityObject->getStr(message, temp);
+ if(temp != "")
+ {
+ _fileName = temp;
+ }
+ }
+ return _fileName;
+
+}*/
+
+bool ClusterTree::transition(int base1, int base2)
+{
+// assumes that the bases of DNA sequences have been translated as
+// a,A = 0; c,C = 1; g,G = 2; t,T,u,U = 3; N = 4;
+// a,A = 0; c,C = 2; g,G = 6; t,T,u,U =17;
+// A <--> G and T <--> C are transitions; all others are transversions.
+ if (((base1 == 0) && (base2 == 6)) || ((base1 == 6) && (base2 == 0)))
+ {
+ return true;
+ }
+ // A <--> G
+ if (((base1 == 17) && (base2 == 2)) || ((base1 == 2) && (base2 == 17)))
+ {
+ return true;
+ }
+ // T <--> C
+ return false;
+}
+
+/**
+ * This function is used to open all the bootstrap tree files. It opens them with the
+ * correct message prompt.
+ */
+bool ClusterTree::openFilesForBootstrap(clustalw::OutputFile* clustalFile, clustalw::OutputFile* phylipFile,
+ clustalw::OutputFile* nexusFile, clustalw::TreeNames* treeNames, string* path)
+{
+ if (clustalw::userParameters->getOutputTreeClustal())
+ {
+ if(!clustalFile || !clustalFile->openFile(&(treeNames->clustalName),
+ bootstrapPrompt, path, "njb", bootstrapFileTypeMsg))
+ {
+ return false;
+ }
+ }
+
+ if (clustalw::userParameters->getOutputTreePhylip())
+ {
+ if(!phylipFile || !phylipFile->openFile(&(treeNames->phylipName),
+ bootstrapPrompt, path, "phb", bootstrapFileTypeMsg))
+ {
+ return false;
+ }
+ }
+
+ if (clustalw::userParameters->getOutputTreeNexus())
+ {
+ if(!nexusFile || !nexusFile->openFile(&(treeNames->nexusName),
+ bootstrapPrompt, path, "treb", bootstrapFileTypeMsg))
+ {
+ return false;
+ }
+ }
+ return true;
+}
+
+bool ClusterTree::openFilesForTreeFromAlignment(clustalw::OutputFile* clustalFile,
+ clustalw::OutputFile* phylipFile, clustalw::OutputFile* distFile, clustalw::OutputFile* nexusFile, clustalw::OutputFile* pimFile,
+ clustalw::TreeNames* treeNames, string* path)
+{
+ if (clustalw::userParameters->getOutputTreeClustal())
+ {
+ if(!clustalFile || !clustalFile->openFile(&(treeNames->clustalName),
+ "\nEnter name for CLUSTAL tree output file ",
+ path, "nj", "Phylogenetic tree"))
+ {
+ return false;
+ }
+ }
+
+ if (clustalw::userParameters->getOutputTreePhylip())
+ {
+ if(!phylipFile || !phylipFile->openFile(&(treeNames->phylipName),
+ "\nEnter name for PHYLIP tree output file ", path, "ph",
+ "Phylogenetic tree"))
+ {
+ return false;
+ }
+ }
+
+ if (clustalw::userParameters->getOutputTreeDistances())
+ {
+ if(!distFile || !distFile->openFile(&(treeNames->distName),
+ "\nEnter name for distance matrix output file ",
+ path, "dst", "Distance matrix"))
+ {
+ return false;
+ }
+ }
+
+ if (clustalw::userParameters->getOutputTreeNexus())
+ {
+ if(!nexusFile || !nexusFile->openFile(&(treeNames->nexusName),
+ "\nEnter name for NEXUS tree output file ", path,
+ "tre", "Nexus tree"))
+ {
+ return false;
+ }
+ }
+
+ if (clustalw::userParameters->getOutputPim())
+ {
+ if(!pimFile || !pimFile->openFile(&(treeNames->pimName),
+ "\nEnter name for % Identity matrix output file ", path, "pim",
+ "perc identity"))
+ {
+ return false;
+ }
+ }
+ return true;
+}
+
+int ClusterTree::calcQuickDistMatForAll(ofstream* clustalFile, ofstream* phylipFile,
+ ofstream* nexusFile, ofstream* pimFile, ofstream* distFile,
+ clustalw::Alignment* alignPtr)
+{
+ int overspill = 0;
+ bool _DNAFlag = clustalw::userParameters->getDNAFlag();
+
+ overspill = calcQuickDistMatForSubSet(clustalFile, phylipFile, nexusFile, alignPtr);
+
+ if (pimFile && clustalw::userParameters->getOutputPim())
+ {
+ verbose = false; // Turn off file output
+ if (_DNAFlag)
+ {
+ calcPercIdentity(pimFile, alignPtr);
+ }
+ else
+ {
+ calcPercIdentity(pimFile, alignPtr);
+ }
+ }
+
+
+ if (distFile && clustalw::userParameters->getOutputTreeDistances())
+ {
+ verbose = false; // Turn off file output
+ if (_DNAFlag)
+ {
+ overspill = dnaDistanceMatrix(distFile, alignPtr);
+ }
+ else
+ {
+ overspill = protDistanceMatrix(distFile, alignPtr);
+ }
+ distanceMatrixOutput(distFile, quickDistMat.get(),
+ alignPtr);
+ }
+ return overspill;
+}
+
+int ClusterTree::calcQuickDistMatForSubSet(ofstream* clustalFile, ofstream* phylipFile,
+ ofstream* nexusFile, clustalw::Alignment* alignPtr,
+ bool inBootLoop)
+{
+ int overspill = 0;
+ bool _DNAFlag = clustalw::userParameters->getDNAFlag();
+
+ if (clustalFile && clustalw::userParameters->getOutputTreeClustal())
+ {
+ if(!inBootLoop)
+ {
+ verbose = true; // Turn on file output
+ }
+ else
+ {
+ verbose = false; // Turn off when we are in the loop in bootstrap!
+ }
+ if (_DNAFlag)
+ {
+ overspill = dnaDistanceMatrix(clustalFile, alignPtr);
+ }
+ else
+ {
+ overspill = protDistanceMatrix(clustalFile, alignPtr);
+ }
+ }
+
+ if (phylipFile && clustalw::userParameters->getOutputTreePhylip())
+ {
+ verbose = false; // Turn off file output
+ if (_DNAFlag)
+ {
+ overspill = dnaDistanceMatrix(phylipFile, alignPtr);
+ }
+ else
+ {
+ overspill = protDistanceMatrix(phylipFile, alignPtr);
+ }
+ }
+
+ if (nexusFile && clustalw::userParameters->getOutputTreeNexus())
+ {
+ verbose = false; // Turn off file output
+ if (_DNAFlag)
+ {
+ overspill = dnaDistanceMatrix(nexusFile, alignPtr);
+ }
+ else
+ {
+ overspill = protDistanceMatrix(nexusFile, alignPtr);
+ }
+ }
+ return overspill;
+}
+
+void ClusterTree::printBootstrapHeaderToClustalFile(ofstream* clustalFile)
+{
+ if(clustalFile)
+ {
+ (*clustalFile) << "\n\n\t\t\tBootstrap Confidence Limits\n\n";
+ (*clustalFile) << "\n Random number generator seed = "
+ << setw(7)
+ << clustalw::userParameters->getBootRanSeed() << "\n";
+ (*clustalFile) << "\n Number of bootstrap trials = " << setw(7)
+ << clustalw::userParameters->getBootNumTrials() << "\n";
+ (*clustalFile) << "\n\n Diagrammatic representation of the above tree: \n";
+ (*clustalFile) << "\n Each row represents 1 tree cycle;";
+ (*clustalFile) << " defining 2 groups.\n";
+ (*clustalFile) << "\n Each column is 1 sequence; ";
+ (*clustalFile) << "the stars in each line show 1 group; ";
+ (*clustalFile) << "\n the dots show the other\n";
+ (*clustalFile) << "\n Numbers show occurences in bootstrap samples.";
+ }
+}
+
+void ClusterTree::promptForBoolSeedAndNumTrials()
+{
+ if (clustalw::userParameters->getMenuFlag())
+ {
+ unsigned int tempSeed;
+ tempSeed = clustalw::utilityObject->getInt(
+ "\n\nEnter seed no. for random number generator ", 1, 1000,
+ clustalw::userParameters->getBootRanSeed());
+ clustalw::userParameters->setBootRanSeed(tempSeed);
+
+ clustalw::userParameters->setBootNumTrials(
+ clustalw::utilityObject->getInt("\n\nEnter number of bootstrap trials ",
+ 1, 10000, clustalw::userParameters->getBootNumTrials()));
+ }
+}
+
+void ClusterTree::printErrorMessageForBootstrap(int totalOverspill, int totalDists,
+ int nfails)
+{
+ cerr << "\n";
+ cerr << "\n WARNING: " << totalOverspill
+ << " of the distances out of a total of "
+ << totalDists << " times" << clustalw::userParameters->getBootNumTrials();
+ cerr << "\n were out of range for the distance correction.";
+ cerr << "\n This affected " << nfails << " out of "
+ << clustalw::userParameters->getBootNumTrials() << " bootstrap trials.";
+ cerr << "\n This may not be fatal but you have been warned!" << "\n";
+ cerr << "\n SUGGESTIONS: 1) turn off the correction";
+ cerr << "\n or 2) remove the most distant sequences";
+ cerr << "\n or 3) use the PHYLIP package.\n\n";
+
+ if (clustalw::userParameters->getMenuFlag())
+ {
+ string dummy;
+ clustalw::utilityObject->getStr(string("Press [RETURN] to continue"), dummy);
+ }
+}
+
+bool ClusterTree::checkIfConditionsMet(int numSeqs, int min)
+{
+ if (clustalw::userParameters->getEmpty())
+ {
+ clustalw::utilityObject->error("You must load an alignment first");
+ return false;
+ }
+
+ if (numSeqs < min)
+ {
+ clustalw::utilityObject->error("Alignment has only %d sequences", numSeqs);
+ return false;
+ }
+
+ return true;
+}
+
+}
git://source.jalview.org / jabaws.git / blobdiff