4 * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
9 #include "RootedClusterTree.h"
10 #include "../../general/OutputFile.h"
11 #include "UPGMAAlgorithm.h"
12 #include "RootedTreeOutput.h"
13 //#include "../RandomGenerator.h"
17 auto_ptr<AlignmentSteps> RootedClusterTree::treeFromDistMatrix(RootedGuideTree* phyloTree,DistMatrix* distMat, Alignment *alignPtr,
18 int seq1, int nSeqs, string& phylipName)
20 OutputFile phylipPhyTreeFile;
21 auto_ptr<AlignmentSteps> progSteps;
24 // Test to see if the inputs are valid
25 if(seq1 < 1 || nSeqs < 1)
27 cerr << "Invalid inputs into treeFromDistMatrix \n"
28 << "seq1 = " << seq1 << " nSeqs = " << nSeqs << "\n"
29 << "Need to end program!\n";
38 lastSeq = firstSeq + nSeqs - 1;
41 info.firstSeq = firstSeq;
42 info.lastSeq = lastSeq;
45 utilityObject->getPath(userParameters->getSeqName(), &path);
49 string name = phylipName;
50 if(!phylipPhyTreeFile.openFile(&name,
51 "\nEnter name for new GUIDE TREE file ", &path, "dnd",
63 RootedTreeOutput outputTree(&info);
65 ofstream* ptrToFile = phylipPhyTreeFile.getPtrToFile();
69 dist = (*distMat)(firstSeq, firstSeq + 1) / 2.0;
70 if(ptrToFile->is_open())
72 (*ptrToFile) << "(" << alignPtr->getName(firstSeq) << ":"
74 << dist << "," << alignPtr->getName(firstSeq + 1) << ":"
75 << setprecision(5) << dist <<");\n";
77 progSteps.reset(new AlignmentSteps);
79 groups.resize(nSeqs + 1, 0);
85 UPGMAAlgorithm clusAlgorithm;
86 progSteps = clusAlgorithm.generateTree(phyloTree, distMat, &info, false);
87 outputTree.printPhylipTree(phyloTree, ptrToFile, alignPtr, distMat);
91 catch(const exception &ex)
93 cerr << "ERROR: Error has occured in treeFromDistMatrix. "
94 << "Need to terminate program.\n"
100 cerr << "ERROR: Error has occured in treeFromDistMatrix. "
101 << "Need to terminate program.\n";
106 void RootedClusterTree::treeFromAlignment(TreeNames* treeNames, Alignment *alignPtr)
110 OutputFile phylipPhyTreeFile;
111 OutputFile clustalPhyTreeFile;
112 OutputFile distancesPhyTreeFile;
113 OutputFile nexusPhyTreeFile;
116 RootedGuideTree phyloTree;
122 numSeqs = alignPtr->getNumSeqs(); // NOTE class variable
125 * Check if numSeqs is ok
127 if(!checkIfConditionsMet(numSeqs, 2))
135 // The SeqInfo struct is passed to reduce the number of parameters passed!
137 info.firstSeq = firstSeq;
138 info.lastSeq = lastSeq;
139 info.numSeqs = numSeqs;
141 RootedTreeOutput outputTree(&info); // No bootstrap!
143 utilityObject->getPath(userParameters->getSeqName(), &path);
146 * Open the required output files.
148 if(!openFilesForTreeFromAlignment(&clustalPhyTreeFile, &phylipPhyTreeFile,
149 &distancesPhyTreeFile, &nexusPhyTreeFile, &pimFile, treeNames, &path))
151 return; // Problem opeing one of the files, cannot continue!
154 int _lenFirstSeq = alignPtr->getSeqLength(firstSeq);
156 bootPositions.clear();
157 bootPositions.resize(_lenFirstSeq + 2);
159 for (j = 1; j <= _lenFirstSeq; ++j)
161 bootPositions[j] = j;
165 * Calculate quickDist and overspill
167 overspill = calcQuickDistMatForAll(clustalPhyTreeFile.getPtrToFile(),
168 phylipPhyTreeFile.getPtrToFile(), nexusPhyTreeFile.getPtrToFile(),
169 pimFile.getPtrToFile(), distancesPhyTreeFile.getPtrToFile(), alignPtr);
171 // check if any distances overflowed the distance corrections
174 totalDists = (numSeqs *(numSeqs - 1)) / 2;
175 overspillMessage(overspill, totalDists);
178 if (userParameters->getOutputTreeClustal())
182 // Turn on file output
185 if (userParameters->getOutputTreeClustal() ||
186 userParameters->getOutputTreePhylip()
187 || userParameters->getOutputTreeNexus())
189 UPGMAAlgorithm clusAlgorithm;
190 clusAlgorithm.setVerbose(true);
191 clusAlgorithm.generateTree(&phyloTree, quickDistMat.get(), &info, false,
192 clustalPhyTreeFile.getPtrToFile());
193 clusAlgorithm.setVerbose(false);
196 if (userParameters->getOutputTreePhylip())
198 outputTree.printPhylipTree(&phyloTree, phylipPhyTreeFile.getPtrToFile(), alignPtr,
202 if (userParameters->getOutputTreeNexus())
204 outputTree.printNexusTree(&phyloTree, nexusPhyTreeFile.getPtrToFile(), alignPtr,
208 /** Free up resources!!!!! */
211 bootPositions.clear();
214 catch(const exception& ex)
216 cerr << ex.what() << endl;
217 utilityObject->error("Terminating program. Cannot continue\n");