unifracunweightedcommand.cpp

   1 /*
   2  *  unifracunweightedcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 2/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "unifracunweightedcommand.h"
  11
  12 /***********************************************************/
  13 UnifracUnweightedCommand::UnifracUnweightedCommand() {
  14         try {
  15                 globaldata = GlobalData::getInstance();
  16
  17                 T = globaldata->gTree;
  18                 tmap = globaldata->gTreemap;
  19                 unweightedFile = globaldata->getTreeFile() + ".unweighted";
  20                 openOutputFile(unweightedFile, out);
  21                 sumFile = globaldata->getTreeFile() + ".uwsummary";
  22                 openOutputFile(sumFile, outSum);
  23                 distFile = globaldata->getTreeFile() + ".uwdistrib";
  24                 openOutputFile(distFile, outDist);
  25
  26                 convert(globaldata->getIters(), iters);  //how many random trees to generate
  27                 unweighted = new Unweighted(tmap);
  28
  29         }
  30         catch(exception& e) {
  31                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function UnifracUnweightedCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  32                 exit(1);
  33         }
  34         catch(...) {
  35                 cout << "An unknown error has occurred in the UnifracUnweightedCommand class function UnifracUnweightedCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  36                 exit(1);
  37         }
  38 }
  39 /***********************************************************/
  40 int UnifracUnweightedCommand::execute() {
  41         try {
  42
  43                 //get unweighted for users tree
  44                 userData.resize(1,0);  //data[0] = unweightedscore
  45                 randomData.resize(1,0); //data[0] = unweightedscore
  46
  47                 //format output
  48                 outDist.setf(ios::fixed, ios::floatfield); outDist.setf(ios::showpoint);
  49                 outDist << "Tree#" << '\t' << "Iter" << '\t' << "UWScore" << endl;
  50
  51                 //create new tree with same num nodes and leaves as users
  52                 randT = new Tree();
  53
  54                 //get pscores for users trees
  55                 for (int i = 0; i < T.size(); i++) {
  56                         cout << "Processing tree " << i+1 << endl;
  57                         userData = unweighted->getValues(T[i]);  //userData[0] = unweightedscore
  58
  59                         //update uscoreFreq
  60                         it = uscoreFreq.find(userData[0]);
  61                         if (it == uscoreFreq.end()) {//new score
  62                                 uscoreFreq[userData[0]] = 1;
  63                         }else{ uscoreFreq[userData[0]]++; }
  64
  65                         //add users score to valid scores
  66                         validScores[userData[0]] = userData[0];
  67
  68                         //saves users score
  69                         utreeScores.push_back(userData[0]);
  70
  71                         //copy T[i]'s info.
  72                         randT->getCopy(T[i]);
  73
  74                         //get unweighted scores for random trees
  75                         for (int j = 0; j < iters; j++) {
  76                                 //create a random tree with same topology as T[i], but different labels
  77                                 randT->assembleRandomUnifracTree();
  78                                 //get pscore of random tree
  79                                 randomData = unweighted->getValues(randT);
  80
  81                                 //add trees unweighted score to map of scores
  82                                 it2 = rscoreFreq.find(randomData[0]);
  83                                 if (it2 != rscoreFreq.end()) {//already have that score
  84                                         rscoreFreq[randomData[0]]++;
  85                                 }else{//first time we have seen this score
  86                                         rscoreFreq[randomData[0]] = 1;
  87                                 }
  88
  89                                 //add randoms score to validscores
  90                                 validScores[randomData[0]] = randomData[0];
  91
  92                                 //output info to uwdistrib file
  93                                 outDist << i+1 << '\t' << '\t'<< j+1 << '\t' << '\t' << randomData[0] << endl;
  94                         }
  95
  96                         //find the signifigance of the score
  97                         float rcumul = 0.0000;
  98                         for (it = rscoreFreq.begin(); it != rscoreFreq.end(); it++) {
  99                                 //get percentage of random trees with that info
 100                                 rscoreFreq[it->first] /= iters;
 101                                 rcumul+= it->second;
 102                                 rCumul[it->first] = rcumul;
 103                         }
 104
 105                         //save the signifigance of the users score for printing later
 106                         UWScoreSig.push_back(rCumul[userData[0]]);
 107
 108                         saveRandomScores(); //save all random scores for unweighted file
 109
 110                         //clear random data
 111                         rscoreFreq.clear();  //you clear this because in the summary file you want the unweighted signifinance to be relative to these 1000 trees.
 112                         rCumul.clear();
 113                 }
 114
 115                 float ucumul = 0.0000;
 116                 float rcumul = 0.0000;
 117                 //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
 118                 for (it = validScores.begin(); it != validScores.end(); it++) {
 119                         it2 = uscoreFreq.find(it->first);
 120                         //user data has that score
 121                         if (it2 != uscoreFreq.end()) { uscoreFreq[it->first] /= T.size(); ucumul+= it2->second;  }
 122                         else { uscoreFreq[it->first] = 0.0000; } //no user trees with that score
 123                         //make uCumul map
 124                         uCumul[it->first] = ucumul;
 125
 126                         //make rscoreFreq map and rCumul
 127                         it2 = totalrscoreFreq.find(it->first);
 128                         //get percentage of random trees with that info
 129                         if (it2 != totalrscoreFreq.end()) {  totalrscoreFreq[it->first] /= (iters*T.size()); rcumul+= it2->second;  }
 130                         else { totalrscoreFreq[it->first] = 0.0000; } //no random trees with that score
 131                         rCumul[it->first] = rcumul;
 132                 }
 133
 134                 printUnweightedFile();
 135                 printUWSummaryFile();
 136
 137                 //reset randomTree parameter to 0
 138                 globaldata->setRandomTree("0");
 139
 140                 delete randT;
 141
 142                 return 0;
 143
 144         }
 145         catch(exception& e) {
 146                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 147                 exit(1);
 148         }
 149         catch(...) {
 150                 cout << "An unknown error has occurred in the UnifracUnweightedCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 151                 exit(1);
 152         }
 153 }
 154 /***********************************************************/
 155 void UnifracUnweightedCommand::printUnweightedFile() {
 156         try {
 157                 //column headers
 158
 159                 out << "Score" << '\t' << "UserFreq" << '\t' << "UserCumul" << '\t' << "RandFreq" << '\t' << "RandCumul" << endl;
 160
 161                 //format output
 162                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 163
 164                 //print each line
 165                 for (it = validScores.begin(); it != validScores.end(); it++) {
 166                         out << setprecision(6) << it->first << '\t' << '\t' << uscoreFreq[it->first] << '\t' << uCumul[it->first] << '\t' << totalrscoreFreq[it->first] << '\t' << rCumul[it->first] << endl;
 167                 }
 168
 169                 out.close();
 170
 171         }
 172         catch(exception& e) {
 173                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function printUnweightedFile. 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 UnifracUnweightedCommand class function printUnweightedFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 178                 exit(1);
 179         }
 180 }
 181
 182 /***********************************************************/
 183 void UnifracUnweightedCommand::printUWSummaryFile() {
 184         try {
 185                 //column headers
 186                 outSum << "Tree#" << '\t'  <<  "UWScore" << '\t' << '\t' << "UWSig" <<  endl;
 187
 188                 //format output
 189                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 190
 191                 //print each line
 192                 for (int i = 0; i< T.size(); i++) {
 193                         outSum << setprecision(6) << i+1 << '\t' << '\t' << utreeScores[i] << '\t' << UWScoreSig[i] << endl;
 194                 }
 195
 196                 outSum.close();
 197         }
 198         catch(exception& e) {
 199                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function printUWSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 200                 exit(1);
 201         }
 202         catch(...) {
 203                 cout << "An unknown error has occurred in the UnifracUnweightedCommand class function printUWSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 204                 exit(1);
 205         }
 206 }
 207 /***********************************************************/
 208 void UnifracUnweightedCommand::saveRandomScores() {
 209         try {
 210                 //update total map with new random scores
 211                 for (it = rscoreFreq.begin(); it != rscoreFreq.end(); it++) {
 212                         //does this score already exist in the total map
 213                         it2 = totalrscoreFreq.find(it->first);
 214                         //if yes then add them
 215                         if (it2 != totalrscoreFreq.end()) {
 216                                 it2->second += it->second;
 217                         }else{ //its a new score
 218                                 totalrscoreFreq[it->first] = 1;
 219                         }
 220                 }
 221         }
 222         catch(exception& e) {
 223                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function saveRandomScores. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 224                 exit(1);
 225         }
 226         catch(...) {
 227                 cout << "An unknown error has occurred in the UnifracUnweightedCommand class function saveRandomScores. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 228                 exit(1);
 229         }
 230 }
 231
 232 /***********************************************************/