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