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");