--- /dev/null
+/**
+ * Author: Mark Larkin
+ *
+ * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
+ */
+#ifdef HAVE_CONFIG_H
+ #include "config.h"
+#endif
+#include "UnRootedClusterTree.h"
+#include "../general/utils.h"
+#include "RandomGenerator.h"
+#include <math.h>
+#include <sstream>
+#include "../general/OutputFile.h"
+
+namespace clustalw
+{
+UnRootedClusterTree::UnRootedClusterTree()
+{
+
+}
+/*
+ * The function treeFromAlignment is called to generate a tree from an alignment
+ * without a distance matrix. It calculates the distance matrix from the alignment
+ * quickly.
+ */
+
+// Note this function was called phylogenetic_tree before
+void UnRootedClusterTree::treeFromAlignment(TreeNames* treeNames, Alignment *alignPtr)
+{
+ try
+ {
+ OutputFile phylipPhyTreeFile;
+ OutputFile clustalPhyTreeFile;
+ OutputFile distancesPhyTreeFile;
+ OutputFile nexusPhyTreeFile;
+ OutputFile pimFile;
+
+ string path;
+ int i, j;
+ int overspill = 0;
+ int totalDists;
+ numSeqs = alignPtr->getNumSeqs(); // NOTE class variable
+
+ /**
+ * Check if numSeqs is ok
+ */
+ if(!checkIfConditionsMet(numSeqs, 2))
+ {
+ return;
+ }
+
+ firstSeq = 1;
+ lastSeq = numSeqs;
+ phyloTree = new PhyloTree;
+
+ // The SeqInfo struct is passed to reduce the number of parameters passed!
+ SeqInfo info;
+ info.firstSeq = firstSeq;
+ info.lastSeq = lastSeq;
+ info.numSeqs = numSeqs;
+
+ outputTree = new ClusterTreeOutput(&info, 0); // No bootstrap!
+
+ phyloTree->treeDesc.resize(numSeqs + 1, vector<int>(numSeqs + 1));
+
+ TreeGroups saveTree(numSeqs + 1, vector<int>(numSeqs + 1));
+
+ // NOTE at the moment there is only one type of clustering algorithm, but there will
+ // be more!
+ clusAlgorithm = new NJTree();
+
+ utilityObject->getPath(userParameters->getSeqName(), &path);
+
+ /**
+ * Open the required output files.
+ */
+ if(!openFilesForTreeFromAlignment(&clustalPhyTreeFile, &phylipPhyTreeFile,
+ &distancesPhyTreeFile, &nexusPhyTreeFile, &pimFile, treeNames, &path))
+ {
+ return; // Problem opeing one of the files, cannot continue!
+ }
+
+ int _lenFirstSeq = alignPtr->getSeqLength(firstSeq);
+
+ bootPositions.resize(_lenFirstSeq + 2);
+
+ for (j = 1; j <= _lenFirstSeq; ++j)
+ {
+ bootPositions[j] = j;
+ }
+
+ /**
+ * Calculate quickDist and overspill
+ */
+ overspill = calcQuickDistMatForAll(clustalPhyTreeFile.getPtrToFile(),
+ phylipPhyTreeFile.getPtrToFile(), nexusPhyTreeFile.getPtrToFile(),
+ pimFile.getPtrToFile(), distancesPhyTreeFile.getPtrToFile(), alignPtr);
+
+ // check if any distances overflowed the distance corrections
+ if (overspill > 0)
+ {
+ totalDists = (numSeqs *(numSeqs - 1)) / 2;
+ overspillMessage(overspill, totalDists);
+ }
+
+ if (userParameters->getOutputTreeClustal())
+ {
+ verbose = true;
+ }
+ // Turn on file output
+
+
+ if (userParameters->getOutputTreeClustal() || userParameters->getOutputTreePhylip()
+ || userParameters->getOutputTreeNexus())
+ {
+
+ clusAlgorithm->setVerbose(true);
+ clusAlgorithm->generateTree(phyloTree, quickDistMat.get(), &info,
+ clustalPhyTreeFile.getPtrToFile());
+ clusAlgorithm->setVerbose(false);
+ }
+
+ for (i = 1; i < numSeqs + 1; i++)
+ for (j = 1; j < numSeqs + 1; j++)
+ {
+ saveTree[i][j] = phyloTree->treeDesc[i][j];
+ }
+
+ if (userParameters->getOutputTreePhylip())
+ {
+ outputTree->printPhylipTree(phyloTree, phylipPhyTreeFile.getPtrToFile(), alignPtr,
+ quickDistMat.get(), &bootTotals);
+ }
+
+ for (i = 1; i < numSeqs + 1; i++)
+ for (j = 1; j < numSeqs + 1; j++)
+ {
+ phyloTree->treeDesc[i][j] = saveTree[i][j];
+ }
+
+ if (userParameters->getOutputTreeNexus())
+ {
+ outputTree->printNexusTree(phyloTree, nexusPhyTreeFile.getPtrToFile(), alignPtr,
+ quickDistMat.get(), &bootTotals);
+ }
+
+ /** Free up resources!!!!! */
+
+ treeGaps.clear();
+ bootPositions.clear();
+
+ delete clusAlgorithm;
+ delete phyloTree;
+ delete outputTree;
+ }
+ catch(const exception& ex)
+ {
+ cerr << ex.what() << endl;
+ utilityObject->error("Terminating program. Cannot continue\n");
+ exit(1);
+ }
+}
+
+/*
+ * Routine for producing unrooted NJ trees from seperately aligned
+ * pairwise distances. This produces the GUIDE DENDROGRAMS in
+ * PHYLIP format.
+ */
+void UnRootedClusterTree::treeFromDistMatrix(DistMatrix* distMat,
+ Alignment *alignPtr, int seq1, int nSeqs,
+ string& phylipName)
+{
+ OutputFile phylipPhyTreeFile;
+
+ try
+ {
+ // Test to see if the inputs are valid
+ if(seq1 < 1 || nSeqs < 1)
+ {
+ cerr << "Invalid inputs into treeFromDistMatrix \n"
+ << "seq1 = " << seq1 << " nSeqs = " << nSeqs << "\n"
+ << "Need to end program!\n";
+ exit(1);
+ return;
+ }
+ PhyloTree phyloTree;
+ float dist;
+ string path;
+ verbose = false;
+ firstSeq = seq1;
+ lastSeq = firstSeq + nSeqs - 1;
+
+ SeqInfo info;
+ info.firstSeq = firstSeq;
+ info.lastSeq = lastSeq;
+ info.numSeqs = nSeqs;
+ // Open up the phylipPhyTreeFile now!
+
+ // NOTE that this is not exactly correct. This may cause a problem when outputing
+ // a tree for each of the profiles. But then we can pass it in, maybe.
+ utilityObject->getPath(userParameters->getSeqName(), &path);
+
+ if(nSeqs >= 2)
+ {
+ string name = phylipName;
+ if(!phylipPhyTreeFile.openFile(&name,
+ "\nEnter name for new GUIDE TREE file ", &path, "dnd",
+ "Guide tree"))
+ {
+ return;
+ }
+ phylipName = name;
+ }
+ else
+ {
+ return;
+ }
+
+ // Not sure about bootstrapping here!
+ clusAlgorithm = new NJTree();
+ outputTree = new ClusterTreeOutput(&info, 0);
+
+ ofstream* ptrToFile = phylipPhyTreeFile.getPtrToFile();
+
+ if (nSeqs == 2)
+ {
+ dist = (*distMat)(firstSeq, firstSeq + 1) / 2.0;
+ if(ptrToFile->is_open())
+ {
+ (*ptrToFile) << "(" << alignPtr->getName(firstSeq) << ":"
+ << setprecision(5)
+ << dist << "," << alignPtr->getName(firstSeq + 1) << ":"
+ << setprecision(5) << dist <<");\n";
+ }
+ }
+ else
+ {
+ int dimensions = lastSeq - firstSeq + 2;
+ phyloTree.treeDesc.resize(dimensions, vector<int>(dimensions));
+/* debugging, also used when -OUTPUTTREE is used */
+#if 0
+ ofstream debuglog("debug.treeFromDistMatrix.txt", ios::out);
+ clusAlgorithm->setVerbose(true);
+ clusAlgorithm->generateTree(&phyloTree, distMat, &info, &debuglog);
+#else
+ clusAlgorithm->generateTree(&phyloTree, distMat, &info);
+#endif
+
+ outputTree->printPhylipTree(&phyloTree, ptrToFile, alignPtr,
+ distMat, &bootTotals);
+ }
+ delete clusAlgorithm;
+ delete outputTree;
+
+ //phylipPhyTreeFile.close();
+ }
+ catch(const exception &ex)
+ {
+ cerr << "ERROR: Error has occured in treeFromDistMatrix. "
+ << "Need to terminate program.\n"
+ << ex.what();
+ exit(1);
+ }
+ catch(...)
+ {
+ cerr << "ERROR: Error has occured in treeFromDistMatrix. "
+ << "Need to terminate program.\n";
+ exit(1);
+ }
+}
+
+/*
+ * Note before I had this function was accepting a distance matrix. But I think it
+ * just uses the quickly generated one. So it doesnt need it anymore. But be careful !!!!
+ */
+void UnRootedClusterTree::bootstrapTree(TreeNames* treeNames, Alignment *alignPtr)
+{
+ int i, j;
+ int ranno;
+ string path;
+
+ OutputFile clustalPhyTreeFile;
+ ofstream* ptrToClustalFile;
+ OutputFile phylipPhyTreeFile;
+ OutputFile nexusPhyTreeFile;
+
+ try
+ {
+ phyloTree = new PhyloTree;
+ PhyloTree sampleTree;
+ PhyloTree standardTree;
+ PhyloTree saveTree;
+ int totalDists, overspill = 0, totalOverspill = 0;
+ int nfails = 0;
+ numSeqs = alignPtr->getNumSeqs();
+ firstSeq = 1;
+ lastSeq = numSeqs;
+ clusAlgorithm = new NJTree();
+
+ SeqInfo info;
+ info.firstSeq = firstSeq;
+ info.lastSeq = lastSeq;
+ info.numSeqs = numSeqs;
+
+ /**
+ * Check if numSeqs is ok
+ */
+ if(!checkIfConditionsMet(numSeqs, 4))
+ {
+ return;
+ }
+
+ if (!userParameters->getOutputTreeClustal() && !userParameters->getOutputTreePhylip()
+ && !userParameters->getOutputTreeNexus())
+ {
+ utilityObject->error("You must select either clustal or phylip or nexus tree output format");
+ return;
+ }
+
+ utilityObject->getPath(userParameters->getSeqName(), &path);
+
+ if(!openFilesForBootstrap(&clustalPhyTreeFile, &phylipPhyTreeFile, &nexusPhyTreeFile,
+ treeNames, &path))
+ {
+ return; // There was a problem opening the output files.
+ }
+
+ int _lenFirstSeq = alignPtr->getSeqLength(firstSeq);
+ bootTotals.clear();
+ bootTotals.resize(numSeqs + 1);
+ bootPositions.clear();
+ bootPositions.resize(_lenFirstSeq + 2);
+
+ for (j = 1; j <= _lenFirstSeq; ++j)
+ // First select all positions for
+ {
+ bootPositions[j] = j;
+ }
+
+ // the "standard" tree
+ overspill = calcQuickDistMatForSubSet(clustalPhyTreeFile.getPtrToFile(),
+ phylipPhyTreeFile.getPtrToFile(), nexusPhyTreeFile.getPtrToFile(), alignPtr);
+
+ // check if any distances overflowed the distance corrections
+ if (overspill > 0)
+ {
+ totalDists = (numSeqs *(numSeqs - 1)) / 2;
+ overspillMessage(overspill, totalDists);
+ }
+
+ treeGaps.clear();
+
+ if (userParameters->getOutputTreeClustal())
+ {
+ verbose = true;
+ }
+ // Turn on screen output
+ phyloTree->treeDesc.resize(numSeqs + 1, vector<int>(numSeqs + 1));
+
+ // compute the standard tree
+
+ if (userParameters->getOutputTreeClustal() || userParameters->getOutputTreePhylip() ||
+ userParameters->getOutputTreeNexus())
+ {
+ clusAlgorithm->setVerbose(true);
+ clusAlgorithm->generateTree(phyloTree, quickDistMat.get(), &info,
+ clustalPhyTreeFile.getPtrToFile());
+ }
+
+ ptrToClustalFile = clustalPhyTreeFile.getPtrToFile();
+
+ promptForBoolSeedAndNumTrials();
+
+ RandomGenerator randGenerator(userParameters->getBootRanSeed());
+
+ /**
+ * Print bootstrap information to top of clustal bootstrap file!
+ */
+ if (userParameters->getOutputTreeClustal())
+ {
+ printBootstrapHeaderToClustalFile(ptrToClustalFile);
+ }
+
+ verbose = false; // Turn OFF screen output
+ clusAlgorithm->setVerbose(false);
+
+ sampleTree.treeDesc.resize(numSeqs + 1, vector<int>(numSeqs + 1));
+
+ if (userParameters->getMenuFlag())
+ {
+ cout << "\n\nEach dot represents 10 trials\n\n";
+ }
+
+ totalOverspill = 0;
+ nfails = 0;
+ int lenSeq1 = alignPtr->getSeqLength(1);
+ for (i = 1; i <= userParameters->getBootNumTrials(); ++i)
+ {
+ for (j = 1; j <= alignPtr->getSeqLength(firstSeq); ++j)
+ {
+ // select alignment positions for
+ ranno = randGenerator.addRand((unsigned long)lenSeq1) + 1;
+ bootPositions[j] = ranno; // bootstrap sample
+ }
+
+ overspill = calcQuickDistMatForSubSet(clustalPhyTreeFile.getPtrToFile(),
+ phylipPhyTreeFile.getPtrToFile(), nexusPhyTreeFile.getPtrToFile(), alignPtr,
+ true);
+
+ if (overspill > 0)
+ {
+ totalOverspill = totalOverspill + overspill;
+ nfails++;
+ }
+
+ treeGaps.clear();
+
+ if (userParameters->getOutputTreeClustal() ||
+ userParameters->getOutputTreePhylip() || userParameters->getOutputTreeNexus())
+ {
+ // NOTE this is bad to pass in clustalPhyTreeFile
+ clusAlgorithm->generateTree(&sampleTree, quickDistMat.get(), &info,
+ clustalPhyTreeFile.getPtrToFile());
+ }
+
+ sampleTree.leftBranch.clear();
+ sampleTree.rightBranch.clear();
+
+ compareTree(phyloTree, &sampleTree, &bootTotals, lastSeq - firstSeq + 1);
+
+ if (userParameters->getMenuFlag())
+ {
+ if (i % 10 == 0)
+ {
+ cout << ".";
+ }
+ if (i % 100 == 0)
+ {
+ cout << "\n";
+ }
+ }
+ }
+
+ // check if any distances overflowed the distance corrections
+ if (nfails > 0)
+ {
+ totalDists = (numSeqs *(numSeqs - 1)) / 2;
+ printErrorMessageForBootstrap(totalOverspill, totalDists, nfails);
+ }
+
+ bootPositions.clear();
+
+ outputTree = new ClusterTreeOutput(&info, userParameters->getBootstrapFormat());
+
+ /**
+ * Print ClustalTree with bootTotals
+ */
+ if (userParameters->getOutputTreeClustal())
+ {
+ outputTree->printTree(phyloTree, clustalPhyTreeFile.getPtrToFile(), &bootTotals);
+ }
+
+ /**
+ * Print phylip tree with boottotals.
+ */
+ if (userParameters->getOutputTreePhylip())
+ {
+ // Save the old tree!
+ saveTree.treeDesc.resize(numSeqs + 1, vector<int>(numSeqs + 1));
+ saveTree.treeDesc.assign(phyloTree->treeDesc.begin(), phyloTree->treeDesc.end());
+
+ outputTree->printPhylipTree(phyloTree, phylipPhyTreeFile.getPtrToFile(),
+ alignPtr, quickDistMat.get(), &bootTotals);
+ // reassign the old values!
+ phyloTree->treeDesc.assign(saveTree.treeDesc.begin(), saveTree.treeDesc.end());
+ }
+
+ /**
+ * print nexus tree with boottotals
+ */
+ if (userParameters->getOutputTreeNexus())
+ {
+ outputTree->printNexusTree(phyloTree, nexusPhyTreeFile.getPtrToFile(), alignPtr,
+ quickDistMat.get(), &bootTotals);
+ }
+
+ delete phyloTree;
+ delete clusAlgorithm;
+ delete outputTree;
+ }
+ catch(const exception &ex)
+ {
+
+ }
+}
+
+}
git://source.jalview.org / jabaws.git / blobdiff