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(){
  15         try {
  16                 globaldata = GlobalData::getInstance();
  17                 t = globaldata->gTree;
  18         }
  19         catch(exception& e) {
  20                 cout << "Standard Error: " << e.what() << " has occurred in the ConcensusCommand class Function ConcensusCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  21                 exit(1);
  22         }
  23         catch(...) {
  24                 cout << "An unknown error has occurred in the ConcensusCommand class function ConcensusCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  25                 exit(1);
  26         }
  27 }
  28 //**********************************************************************************************************************
  29
  30 ConcensusCommand::~ConcensusCommand(){}
  31
  32 //**********************************************************************************************************************
  33
  34 int ConcensusCommand::execute(){
  35         try {
  36
  37                 if (t.size() == 0) { return 0; }
  38                 else {
  39                         numNodes = t[0]->getNumNodes();
  40                         numLeaves = t[0]->getNumLeaves();
  41                 }
  42
  43                 //initialize nodepairs
  44                 nodePairs.resize(numNodes-numLeaves);
  45
  46                 //process each tree and fill counts
  47                 for (int i = 0; i < t.size(); i++) {
  48                         //process each nonleaf node
  49                         for (int j = numLeaves; j < numNodes; j++) {
  50                                 //get pairing
  51                                 int leftChild = t[i]->tree[j].getLChild();
  52                                 int rightChild = t[i]->tree[j].getRChild();
  53                                 string pair = toString(leftChild) + "-" + toString(rightChild);
  54
  55                                 //if this is an existing pairing for this node then increment the count otherwise add new pairing.
  56                                 it = nodePairs[j-numLeaves].find(pair);
  57                                 if (it != nodePairs[j-numLeaves].end()) {//already have that score
  58                                         nodePairs[j-numLeaves][pair]++;
  59                                 }else{//first time we have seen this score
  60                                         nodePairs[j-numLeaves][pair] = 1;
  61                                 }
  62                         }
  63                 }
  64
  65
  66                 //print out pairings for testing
  67                 /*for (int i = 0; i < nodePairs.size(); i++) {
  68                         cout << "pairs for node " << i+numLeaves << endl;
  69                         for (it = nodePairs[i].begin(); it != nodePairs[i].end(); it++) {
  70                                 cout << " pair = " << it->first <<  " count = " << it->second << endl;
  71                         }
  72                 }*/
  73
  74                 //open file for pairing not included in the tree
  75                 notIncluded = getRootName(globaldata->inputFileName) + "concensuspairs";
  76                 openOutputFile(notIncluded, out2);
  77
  78                 concensusTree = new Tree();
  79
  80                 //set relationships for nonleaf nodes
  81                 for (int j = numLeaves; j < numNodes; j++) {
  82
  83                         //find that nodes pairing with the highest count
  84                         int large = 0;
  85                         for (it = nodePairs[j-numLeaves].begin(); it != nodePairs[j-numLeaves].end(); it++) {
  86                                 if (it->second > large) { large = it->second;  it2 = it; }
  87                         }
  88
  89                         string pair = it2->first;
  90                         int pos = pair.find_first_of('-');
  91                         string lefts =  pair.substr(0, pos);
  92                         string rights =  pair.substr(pos+1, pair.length());
  93
  94                         //converts string to int
  95                         int left = atoi(lefts.c_str());
  96                         int right = atoi(rights.c_str());
  97
  98                         //set parents and children
  99                         concensusTree->tree[j].setChildren(left, right);
 100                         concensusTree->tree[left].setParent(j);
 101                         concensusTree->tree[right].setParent(j);
 102
 103                         // set Branchlengths //
 104                         //if your children are leaves remove their branchlengths
 105                         if (concensusTree->tree[left].getLChild() == -1) {  concensusTree->tree[left].setBranchLength(1.0); }
 106                         if (concensusTree->tree[right].getLChild() == -1) {  concensusTree->tree[right].setBranchLength(1.0); }
 107
 108                         //set your branch length to the percentage of trees with this pairing
 109                         float bLength = (float) it2->second / (float) t.size();
 110                         concensusTree->tree[j].setBranchLength(bLength);
 111
 112                         //print out set used
 113                         string leftName, rightName;
 114                         getNames(it2->first, leftName, rightName);
 115
 116                         out2 << "Node " << j+1 << " in concensus tree: " << leftName << "-" << rightName << '\t' << (float)it2->second / (float) t.size() << endl;
 117                         out2 << "Node " << j+1 << " sets not included in concensus tree: " << endl;
 118
 119                         //print out sets not included
 120                         for (it = nodePairs[j-numLeaves].begin(); it != nodePairs[j-numLeaves].end(); it++) {
 121                                 if (it != it2) {
 122                                         getNames(it->first, leftName, rightName);
 123                                         out2 << leftName << "-" << rightName << '\t' << (float)it->second / (float) t.size() << endl;
 124                                 }
 125                         }
 126                         out2 << endl;
 127
 128                 }
 129
 130                 concensusTree->assembleTree();
 131
 132                 outputFile = getRootName(globaldata->inputFileName) + "concensus.tre";
 133                 openOutputFile(outputFile, out);
 134
 135                 concensusTree->print(out);
 136
 137                 out.close(); out2.close();
 138
 139                 delete concensusTree;
 140
 141                 return 0;
 142         }
 143         catch(exception& e) {
 144                 cout << "Standard Error: " << e.what() << " has occurred in the ConcensusCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 145                 exit(1);
 146         }
 147         catch(...) {
 148                 cout << "An unknown error has occurred in the ConcensusCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 149                 exit(1);
 150         }
 151 }
 152 //**********************************************************************************************************************
 153
 154 void ConcensusCommand::getNames(string pair, string& leftName, string& rightName) {
 155         try {
 156                 int pos = pair.find_first_of('-');
 157                 string lefts =  pair.substr(0, pos);
 158                 string rights =  pair.substr(pos+1, pair.length());
 159
 160                 //converts string to int
 161                 int     left = atoi(lefts.c_str());
 162                 int     right = atoi(rights.c_str());
 163
 164                 //get name
 165                 leftName = concensusTree->tree[left].getName();
 166                 rightName = concensusTree->tree[right].getName();
 167
 168                 //if you are not a leaf use node number as name
 169                 if (leftName == "") {  leftName = toString(left+1); }
 170                 if (rightName == "") {  rightName = toString(right+1); }
 171         }
 172         catch(exception& e) {
 173                 cout << "Standard Error: " << e.what() << " has occurred in the ConcensusCommand class Function getNames. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 174                 exit(1);
 175         }
 176         catch(...) {
 177                 cout << "An unknown error has occurred in the ConcensusCommand class function getNames. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 178                 exit(1);
 179         }
 180 }
 181
 182 //**********************************************************************************************************************
 183
 184