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