concensuscommand.cpp

   1 /*
   2  *  concensuscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 4/29/09.
   6  *  Copyright 2009 Schloss Lab UMASS AMherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "concensuscommand.h"
  11
  12 //**********************************************************************************************************************
  13
  14 ConcensusCommand::ConcensusCommand(string fileroot){
  15         try {
  16                 globaldata = GlobalData::getInstance();
  17                 abort = false;
  18
  19                 filename = fileroot;
  20                 //allow user to run help
  21                 //if(option == "help") { help(); abort = true; }
  22
  23                 //else {
  24                         //if (option != "") { mothurOut("There are no valid parameters for the concensus command."); mothurOutEndLine(); abort = true; }
  25
  26                 //      //no trees were read
  27                 //      if (globaldata->gTree.size() == 0) {  mothurOut("You must execute the read.tree command, before you may use the concensus command."); mothurOutEndLine(); abort = true;  }
  28                         //else {
  29                          t = globaldata->gTree;
  30                          //     }
  31                 //}
  32         }
  33         catch(exception& e) {
  34                 errorOut(e, "ConcensusCommand", "ConcensusCommand");
  35                 exit(1);
  36         }
  37 }
  38
  39 //**********************************************************************************************************************
  40
  41 void ConcensusCommand::help(){
  42         try {
  43                 mothurOut("The concensus command can only be executed after a successful read.tree command.\n");
  44                 mothurOut("The concensus command has no parameters.\n");
  45                 mothurOut("The concensus command should be in the following format: concensus().\n");
  46                 mothurOut("The concensus command output two files: .concensus.tre and .concensuspairs.\n");
  47                 mothurOut("The .concensus.tre file contains the concensus 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 .concensuspairs file contains a list of the internal nodes in your tree.  For each node, the pair that was used in the concensus 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");
  51         }
  52         catch(exception& e) {
  53                 errorOut(e, "ConcensusCommand", "help");
  54                 exit(1);
  55         }
  56 }
  57
  58 //**********************************************************************************************************************
  59
  60 ConcensusCommand::~ConcensusCommand(){}
  61
  62 //**********************************************************************************************************************
  63
  64 int ConcensusCommand::execute(){
  65         try {
  66
  67                 if (abort == true) { return 0; }
  68                 else {
  69                         numNodes = t[0]->getNumNodes();
  70                         numLeaves = t[0]->getNumLeaves();
  71                 }
  72
  73                 //get the possible pairings
  74                 getSets();
  75
  76                 //open file for pairing not included in the tree
  77                 notIncluded = filename + ".concensuspairs";
  78                 openOutputFile(notIncluded, out2);
  79
  80                 concensusTree = new Tree();
  81
  82                 it2 = nodePairs.find(treeSet);
  83
  84                 nodePairsInTree[treeSet] = it2->second;
  85
  86                 //erase treeset because you are adding it
  87                 nodePairs.erase(treeSet);
  88
  89                 //set count to numLeaves;
  90                 count = numLeaves;
  91
  92                 buildConcensusTree(treeSet);
  93
  94                 concensusTree->assembleTree();
  95
  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; }
  99
 100                 vector<string> temp;
 101                 //output sets included
 102                 out2 << endl << "Sets included in the concensus tree:" << endl << endl;
 103                 for (it2 = nodePairsInTree.begin(); it2 != nodePairsInTree.end(); it2++) {
 104                         //only output pairs not leaves
 105                         if (it2->first.size() > 1) {
 106                                 temp.clear();
 107                                 //initialize temp to all "."
 108                                 temp.resize(treeSet.size(), ".");
 109
 110                                 //set the spot in temp that represents it2->first[i] to a "*"
 111                                 for (int i = 0; i < it2->first.size(); i++) {
 112                                         //find spot
 113                                         int index = findSpot(it2->first[i]);
 114                                         temp[index] = "*";
 115                                 }
 116
 117                                 //output temp
 118                                 for (int j = 0; j < temp.size(); j++) {
 119                                         out2 << temp[j];
 120                                 }
 121                                 out2 << '\t' << it2->second << endl;
 122                         }
 123                 }
 124
 125                 //output sets not included
 126                 out2 << endl << "Sets NOT included in the concensus tree:" << endl << endl;
 127                 for (it2 = nodePairs.begin(); it2 != nodePairs.end(); it2++) {
 128                         temp.clear();
 129                         //initialize temp to all "."
 130                         temp.resize(treeSet.size(), ".");
 131
 132                         //set the spot in temp that represents it2->first[i] to a "*"
 133                         for (int i = 0; i < it2->first.size(); i++) {
 134                                 //find spot
 135                                 int index = findSpot(it2->first[i]);
 136                                 temp[index] = "*";
 137                         }
 138
 139                         //output temp
 140                         for (int j = 0; j < temp.size(); j++) {
 141                                 out2 << temp[j];
 142                         }
 143                         out2 << '\t' << it2->second << endl;
 144                 }
 145
 146                 outputFile = filename + ".cons.tre";
 147                 openOutputFile(outputFile, out);
 148
 149                 concensusTree->printForBoot(out);
 150
 151                 out.close(); out2.close();
 152
 153                 delete concensusTree;
 154
 155                 return 0;
 156         }
 157         catch(exception& e) {
 158                 errorOut(e, "ConcensusCommand", "execute");
 159                 exit(1);
 160         }
 161 }
 162
 163 //**********************************************************************************************************************
 164 int ConcensusCommand::buildConcensusTree(vector<string> nodeSet) {
 165         try {
 166                 vector<string> leftChildSet;
 167                 vector<string> rightChildSet;
 168
 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 concensusTree->getIndex(nodeSet[0]);
 173                 //terminate recursion
 174                 }else if (count == numNodes) { return 0; }
 175                 else {
 176                         leftChildSet = getNextAvailableSet(nodeSet);
 177                         rightChildSet = getRestSet(nodeSet, leftChildSet);
 178                         int left = buildConcensusTree(leftChildSet);
 179                         int right = buildConcensusTree(rightChildSet);
 180                         concensusTree->tree[count].setChildren(left, right);
 181                         concensusTree->tree[count].setLabel(nodePairsInTree[nodeSet]);
 182                         concensusTree->tree[left].setParent(count);
 183                         concensusTree->tree[right].setParent(count);
 184                         count++;
 185                         return (count-1);
 186                 }
 187
 188         }
 189         catch(exception& e) {
 190                 errorOut(e, "ConcensusCommand", "buildConcensusTree");
 191                 exit(1);
 192         }
 193 }
 194
 195 //**********************************************************************************************************************
 196 void ConcensusCommand::getSets() {
 197         try {
 198                 vector<string> temp;
 199                 treeSet.clear();
 200
 201                 //for each tree add the possible pairs you find
 202                 for (int i = 0; i < t.size(); i++) {
 203
 204                         //for each non-leaf node get descendant info.
 205                         for (int j = numLeaves; j < numNodes; j++) {
 206                                 temp.clear();
 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);
 210                                 }
 211
 212                                 //sort temp
 213                                 sort(temp.begin(), temp.end());
 214
 215                                 it2 = nodePairs.find(temp);
 216                                 if (it2 != nodePairs.end()) {
 217                                         nodePairs[temp]++;
 218                                 }else{
 219                                         nodePairs[temp] = 1;
 220                                 }
 221                         }
 222                 }
 223
 224                 //add each leaf to terminate recursion in concensus
 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++) {
 228
 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);
 232
 233                         //fill treeSet
 234                         treeSet.push_back(it->first);
 235
 236                         //add leaf to list but with sighting value less then all non leaf pairs
 237                         nodePairs[temp] = 0;
 238                 }
 239
 240                 sort(treeSet.begin(), treeSet.end());
 241         }
 242         catch(exception& e) {
 243                 errorOut(e, "ConcensusCommand", "getSets");
 244                 exit(1);
 245         }
 246 }
 247
 248 //**********************************************************************************************************************
 249 vector<string> ConcensusCommand::getNextAvailableSet(vector<string> bigset) {
 250         try {
 251                 vector<string> largest; largest.clear();
 252                 int largestSighting = -1;
 253
 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) {
 258
 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; }
 263                                 }
 264
 265                         }
 266                 }
 267
 268                 //save for printing out later and for branch lengths
 269                 nodePairsInTree[largest] = nodePairs[largest];
 270
 271                 //delete whatever set you return because it is no longer available
 272                 nodePairs.erase(largest);
 273
 274                 return largest;
 275
 276         }
 277         catch(exception& e) {
 278                 errorOut(e, "ConcensusCommand", "getNextAvailableSet");
 279                 exit(1);
 280         }
 281 }
 282
 283 //**********************************************************************************************************************
 284 vector<string> ConcensusCommand::getRestSet(vector<string> bigset, vector<string> subset) {
 285         try {
 286                 vector<string> rest;
 287
 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; }
 292                         }
 293
 294                         //its not in the subset so put it in the rest
 295                         if (inSubset == false) { rest.push_back(bigset[i]); }
 296                 }
 297
 298                 //save for printing out later and for branch lengths
 299                 nodePairsInTree[rest] = nodePairs[rest];
 300
 301                 //delete whatever set you return because it is no longer available
 302                 nodePairs.erase(rest);
 303
 304                 return rest;
 305
 306         }
 307         catch(exception& e) {
 308                 errorOut(e, "ConcensusCommand", "getRestSet");
 309                 exit(1);
 310         }
 311 }
 312
 313 //**********************************************************************************************************************
 314 bool ConcensusCommand::isSubset(vector<string> bigset, vector<string> subset) {
 315         try {
 316
 317                 //check if each guy in suset is also in bigset
 318                 for (int i = 0; i < subset.size(); i++) {
 319                         bool match = false;
 320                         for (int j = 0; j < bigset.size(); j++) {
 321                                 if (subset[i] == bigset[j]) { match = true; break; }
 322                         }
 323
 324                         //you have a guy in subset that had no match in bigset
 325                         if (match == false) { return false; }
 326                 }
 327
 328                 return true;
 329
 330         }
 331         catch(exception& e) {
 332                 errorOut(e, "ConcensusCommand", "isSubset");
 333                 exit(1);
 334         }
 335 }
 336 //**********************************************************************************************************************
 337 int ConcensusCommand::findSpot(string node) {
 338         try {
 339                 int spot;
 340
 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; }
 344                 }
 345
 346                 return spot;
 347
 348         }
 349         catch(exception& e) {
 350                 errorOut(e, "ConcensusCommand", "findSpot");
 351                 exit(1);
 352         }
 353 }
 354 //**********************************************************************************************************************
 355
 356
 357
 358