1 /**
2 * Author: Mark Larkin
3 *
4 * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
5 */
6 #ifdef HAVE_CONFIG_H
7 #include "config.h"
8 #endif
9 #include "RootedClusterTree.h"
10 #include "../../general/OutputFile.h"
11 #include "UPGMAAlgorithm.h"
12 #include "RootedTreeOutput.h"
13 //#include "../RandomGenerator.h"
14 namespace clustalw
15 {
16
treeFromDistMatrix(RootedGuideTree * phyloTree,DistMatrix * distMat,Alignment * alignPtr,int seq1,int nSeqs,string & phylipName)17 auto_ptr<AlignmentSteps> RootedClusterTree::treeFromDistMatrix(RootedGuideTree* phyloTree,DistMatrix* distMat, Alignment *alignPtr,
18 int seq1, int nSeqs, string& phylipName)
19 {
20 OutputFile phylipPhyTreeFile;
21 auto_ptr<AlignmentSteps> progSteps;
22 try
23 {
24 // Test to see if the inputs are valid
25 if(seq1 < 1 || nSeqs < 1)
26 {
27 cerr << "Invalid inputs into treeFromDistMatrix \n"
28 << "seq1 = " << seq1 << " nSeqs = " << nSeqs << "\n"
29 << "Need to end program!\n";
30 exit(1);
31 return progSteps;
32 }
33
34 float dist;
35 string path;
36 verbose = false;
37 firstSeq = seq1;
38 lastSeq = firstSeq + nSeqs - 1;
39
40 SeqInfo info;
41 info.firstSeq = firstSeq;
42 info.lastSeq = lastSeq;
43 info.numSeqs = nSeqs;
44
45 utilityObject->getPath(userParameters->getSeqName(), &path);
46
47 if(nSeqs >= 2)
48 {
49 string name = phylipName;
50 if(!phylipPhyTreeFile.openFile(&name,
51 "\nEnter name for new GUIDE TREE file ", &path, "dnd",
52 "Guide tree"))
53 {
54 return progSteps;
55 }
56 phylipName = name;
57 }
58 else
59 {
60 return progSteps;
61 }
62
63 RootedTreeOutput outputTree(&info);
64
65 ofstream* ptrToFile = phylipPhyTreeFile.getPtrToFile();
66
67 if (nSeqs == 2)
68 {
69 dist = (*distMat)(firstSeq, firstSeq + 1) / 2.0;
70 if(ptrToFile->is_open())
71 {
72 (*ptrToFile) << "(" << alignPtr->getName(firstSeq) << ":"
73 << setprecision(5)
74 << dist << "," << alignPtr->getName(firstSeq + 1) << ":"
75 << setprecision(5) << dist <<");\n";
76 }
77 progSteps.reset(new AlignmentSteps);
78 vector<int> groups;
79 groups.resize(nSeqs + 1, 0);
80 groups[1] = 1;
81 groups[2] = 2;
82 }
83 else
84 {
85 UPGMAAlgorithm clusAlgorithm;
86 progSteps = clusAlgorithm.generateTree(phyloTree, distMat, &info, false);
87 outputTree.printPhylipTree(phyloTree, ptrToFile, alignPtr, distMat);
88 }
89 return progSteps;
90 }
91 catch(const exception &ex)
92 {
93 cerr << "ERROR: Error has occured in treeFromDistMatrix. "
94 << "Need to terminate program.\n"
95 << ex.what();
96 exit(1);
97 }
98 catch(...)
99 {
100 cerr << "ERROR: Error has occured in treeFromDistMatrix. "
101 << "Need to terminate program.\n";
102 exit(1);
103 }
104 }
105
treeFromAlignment(TreeNames * treeNames,Alignment * alignPtr)106 void RootedClusterTree::treeFromAlignment(TreeNames* treeNames, Alignment *alignPtr)
107 {
108 try
109 {
110 OutputFile phylipPhyTreeFile;
111 OutputFile clustalPhyTreeFile;
112 OutputFile distancesPhyTreeFile;
113 OutputFile nexusPhyTreeFile;
114 OutputFile pimFile;
115
116 RootedGuideTree phyloTree;
117
118 string path;
119 int j;
120 int overspill = 0;
121 int totalDists;
122 numSeqs = alignPtr->getNumSeqs(); // NOTE class variable
123
124 /**
125 * Check if numSeqs is ok
126 */
127 if(!checkIfConditionsMet(numSeqs, 2))
128 {
129 return;
130 }
131
132 firstSeq = 1;
133 lastSeq = numSeqs;
134
135 // The SeqInfo struct is passed to reduce the number of parameters passed!
136 SeqInfo info;
137 info.firstSeq = firstSeq;
138 info.lastSeq = lastSeq;
139 info.numSeqs = numSeqs;
140
141 RootedTreeOutput outputTree(&info); // No bootstrap!
142
143 utilityObject->getPath(userParameters->getSeqName(), &path);
144
145 /**
146 * Open the required output files.
147 */
148 if(!openFilesForTreeFromAlignment(&clustalPhyTreeFile, &phylipPhyTreeFile,
149 &distancesPhyTreeFile, &nexusPhyTreeFile, &pimFile, treeNames, &path))
150 {
151 return; // Problem opeing one of the files, cannot continue!
152 }
153
154 int _lenFirstSeq = alignPtr->getSeqLength(firstSeq);
155
156 bootPositions.clear();
157 bootPositions.resize(_lenFirstSeq + 2);
158
159 for (j = 1; j <= _lenFirstSeq; ++j)
160 {
161 bootPositions[j] = j;
162 }
163
164 /**
165 * Calculate quickDist and overspill
166 */
167 overspill = calcQuickDistMatForAll(clustalPhyTreeFile.getPtrToFile(),
168 phylipPhyTreeFile.getPtrToFile(), nexusPhyTreeFile.getPtrToFile(),
169 pimFile.getPtrToFile(), distancesPhyTreeFile.getPtrToFile(), alignPtr);
170
171 // check if any distances overflowed the distance corrections
172 if (overspill > 0)
173 {
174 totalDists = (numSeqs *(numSeqs - 1)) / 2;
175 overspillMessage(overspill, totalDists);
176 }
177
178 if (userParameters->getOutputTreeClustal())
179 {
180 verbose = true;
181 }
182 // Turn on file output
183
184
185 if (userParameters->getOutputTreeClustal() ||
186 userParameters->getOutputTreePhylip()
187 || userParameters->getOutputTreeNexus())
188 {
189 UPGMAAlgorithm clusAlgorithm;
190 clusAlgorithm.setVerbose(true);
191 clusAlgorithm.generateTree(&phyloTree, quickDistMat.get(), &info, false,
192 clustalPhyTreeFile.getPtrToFile());
193 clusAlgorithm.setVerbose(false);
194 }
195
196 if (userParameters->getOutputTreePhylip())
197 {
198 outputTree.printPhylipTree(&phyloTree, phylipPhyTreeFile.getPtrToFile(), alignPtr,
199 quickDistMat.get());
200 }
201
202 if (userParameters->getOutputTreeNexus())
203 {
204 outputTree.printNexusTree(&phyloTree, nexusPhyTreeFile.getPtrToFile(), alignPtr,
205 quickDistMat.get());
206 }
207
208 /** Free up resources!!!!! */
209
210 treeGaps.clear();
211 bootPositions.clear();
212
213 }
214 catch(const exception& ex)
215 {
216 cerr << ex.what() << endl;
217 utilityObject->error("Terminating program. Cannot continue\n");
218 exit(1);
219 }
220 }
221
222 }
223