5 * Created by Sarah Westcott on 4/29/09.
6 * Copyright 2009 Schloss Lab UMASS AMherst. All rights reserved.
10 #include "consensuscommand.h"
12 //**********************************************************************************************************************
14 ConcensusCommand::ConcensusCommand(string fileroot){
16 globaldata = GlobalData::getInstance();
20 //allow user to run help
21 //if(option == "help") { help(); abort = true; }
24 //if (option != "") { mothurOut("There are no valid parameters for the consensus command."); mothurOutEndLine(); abort = true; }
26 // //no trees were read
27 // if (globaldata->gTree.size() == 0) { mothurOut("You must execute the read.tree command, before you may use the consensus command."); mothurOutEndLine(); abort = true; }
29 t = globaldata->gTree;
34 errorOut(e, "ConcensusCommand", "ConcensusCommand");
39 //**********************************************************************************************************************
41 void ConcensusCommand::help(){
43 mothurOut("The consensus command can only be executed after a successful read.tree command.\n");
44 mothurOut("The consensus command has no parameters.\n");
45 mothurOut("The consensus command should be in the following format: consensus().\n");
46 mothurOut("The consensus command output two files: .consensus.tre and .consensuspairs.\n");
47 mothurOut("The .consensus.tre file contains the consensus tree of the trees in your input file.\n");
48 mothurOut("The branch lengths are the percentage of trees in your input file that had the given pair.\n");
49 mothurOut("The .consensuspairs file contains a list of the internal nodes in your tree. For each node, the pair that was used in the consensus tree \n");
50 mothurOut("is reported with its percentage, as well as the other pairs that were seen for that node but not used and their percentages.\n\n");
53 errorOut(e, "ConcensusCommand", "help");
58 //**********************************************************************************************************************
60 ConcensusCommand::~ConcensusCommand(){}
62 //**********************************************************************************************************************
64 int ConcensusCommand::execute(){
67 if (abort == true) { return 0; }
69 numNodes = t[0]->getNumNodes();
70 numLeaves = t[0]->getNumLeaves();
73 //get the possible pairings
76 //open file for pairing not included in the tree
77 notIncluded = filename + ".cons.pairs";
78 openOutputFile(notIncluded, out2);
80 consensusTree = new Tree();
82 it2 = nodePairs.find(treeSet);
84 nodePairsInTree[treeSet] = it2->second;
86 //erase treeset because you are adding it
87 nodePairs.erase(treeSet);
89 //set count to numLeaves;
92 buildConcensusTree(treeSet);
94 consensusTree->assembleTree();
96 //output species in order
97 out2 << "Species in Order: " << endl << endl;
98 for (int i = 0; i < treeSet.size(); i++) { out2 << i+1 << ". " << treeSet[i] << endl; }
101 //output sets included
102 out2 << endl << "Sets included in the consensus tree:" << endl << endl;
103 for (it2 = nodePairsInTree.begin(); it2 != nodePairsInTree.end(); it2++) {
104 //only output pairs not leaves
105 if (it2->first.size() > 1) {
107 //initialize temp to all "."
108 temp.resize(treeSet.size(), ".");
110 //set the spot in temp that represents it2->first[i] to a "*"
111 for (int i = 0; i < it2->first.size(); i++) {
113 int index = findSpot(it2->first[i]);
118 for (int j = 0; j < temp.size(); j++) {
121 out2 << '\t' << it2->second << endl;
125 //output sets not included
126 out2 << endl << "Sets NOT included in the consensus tree:" << endl << endl;
127 for (it2 = nodePairs.begin(); it2 != nodePairs.end(); it2++) {
129 //initialize temp to all "."
130 temp.resize(treeSet.size(), ".");
132 //set the spot in temp that represents it2->first[i] to a "*"
133 for (int i = 0; i < it2->first.size(); i++) {
135 int index = findSpot(it2->first[i]);
140 for (int j = 0; j < temp.size(); j++) {
143 out2 << '\t' << it2->second << endl;
146 outputFile = filename + ".cons.tre";
147 openOutputFile(outputFile, out);
149 consensusTree->printForBoot(out);
151 out.close(); out2.close();
153 delete consensusTree;
157 catch(exception& e) {
158 errorOut(e, "ConcensusCommand", "execute");
163 //**********************************************************************************************************************
164 int ConcensusCommand::buildConcensusTree(vector<string> nodeSet) {
166 vector<string> leftChildSet;
167 vector<string> rightChildSet;
169 //if you are at a leaf
170 if (nodeSet.size() == 1) {
171 //return the vector index of the leaf you are at
172 return consensusTree->getIndex(nodeSet[0]);
173 //terminate recursion
174 }else if (count == numNodes) { return 0; }
176 leftChildSet = getNextAvailableSet(nodeSet);
177 rightChildSet = getRestSet(nodeSet, leftChildSet);
178 int left = buildConcensusTree(leftChildSet);
179 int right = buildConcensusTree(rightChildSet);
180 consensusTree->tree[count].setChildren(left, right);
181 consensusTree->tree[count].setLabel(nodePairsInTree[nodeSet]);
182 consensusTree->tree[left].setParent(count);
183 consensusTree->tree[right].setParent(count);
189 catch(exception& e) {
190 errorOut(e, "ConcensusCommand", "buildConcensusTree");
195 //**********************************************************************************************************************
196 void ConcensusCommand::getSets() {
201 //for each tree add the possible pairs you find
202 for (int i = 0; i < t.size(); i++) {
204 //for each non-leaf node get descendant info.
205 for (int j = numLeaves; j < numNodes; j++) {
207 //go through pcounts and pull out descendants
208 for (it = t[i]->tree[j].pcount.begin(); it != t[i]->tree[j].pcount.end(); it++) {
209 temp.push_back(it->first);
213 sort(temp.begin(), temp.end());
215 it2 = nodePairs.find(temp);
216 if (it2 != nodePairs.end()) {
224 //add each leaf to terminate recursion in consensus
225 //you want the leaves in there but with insignifigant sightings value so it is added last
226 //for each leaf node get descendant info.
227 for (int j = 0; j < numLeaves; j++) {
229 //only need the first one since leaves have no descendants but themselves
230 it = t[0]->tree[j].pcount.begin();
231 temp.clear(); temp.push_back(it->first);
234 treeSet.push_back(it->first);
236 //add leaf to list but with sighting value less then all non leaf pairs
240 sort(treeSet.begin(), treeSet.end());
242 catch(exception& e) {
243 errorOut(e, "ConcensusCommand", "getSets");
248 //**********************************************************************************************************************
249 vector<string> ConcensusCommand::getNextAvailableSet(vector<string> bigset) {
251 vector<string> largest; largest.clear();
252 int largestSighting = -1;
254 //go through the sets
255 for (it2 = nodePairs.begin(); it2 != nodePairs.end(); it2++) {
256 //are you a subset of bigset
257 if (isSubset(bigset, it2->first) == true) {
259 //are you the largest. if you are the same size as current largest refer to sighting
260 if (it2->first.size() > largest.size()) { largest = it2->first; largestSighting = it2->second; }
261 else if (it2->first.size() == largest.size()) {
262 if (it2->second > largestSighting) { largest = it2->first; largestSighting = it2->second; }
268 //save for printing out later and for branch lengths
269 nodePairsInTree[largest] = nodePairs[largest];
271 //delete whatever set you return because it is no longer available
272 nodePairs.erase(largest);
277 catch(exception& e) {
278 errorOut(e, "ConcensusCommand", "getNextAvailableSet");
283 //**********************************************************************************************************************
284 vector<string> ConcensusCommand::getRestSet(vector<string> bigset, vector<string> subset) {
288 for (int i = 0; i < bigset.size(); i++) {
289 bool inSubset = false;
290 for (int j = 0; j < subset.size(); j++) {
291 if (bigset[i] == subset[j]) { inSubset = true; break; }
294 //its not in the subset so put it in the rest
295 if (inSubset == false) { rest.push_back(bigset[i]); }
298 //save for printing out later and for branch lengths
299 nodePairsInTree[rest] = nodePairs[rest];
301 //delete whatever set you return because it is no longer available
302 nodePairs.erase(rest);
307 catch(exception& e) {
308 errorOut(e, "ConcensusCommand", "getRestSet");
313 //**********************************************************************************************************************
314 bool ConcensusCommand::isSubset(vector<string> bigset, vector<string> subset) {
317 //check if each guy in suset is also in bigset
318 for (int i = 0; i < subset.size(); i++) {
320 for (int j = 0; j < bigset.size(); j++) {
321 if (subset[i] == bigset[j]) { match = true; break; }
324 //you have a guy in subset that had no match in bigset
325 if (match == false) { return false; }
331 catch(exception& e) {
332 errorOut(e, "ConcensusCommand", "isSubset");
336 //**********************************************************************************************************************
337 int ConcensusCommand::findSpot(string node) {
341 //check if each guy in suset is also in bigset
342 for (int i = 0; i < treeSet.size(); i++) {
343 if (treeSet[i] == node) { spot = i; break; }
349 catch(exception& e) {
350 errorOut(e, "ConcensusCommand", "findSpot");
354 //**********************************************************************************************************************