unifracweightedcommand.cpp

   1 /*
   2  *  unifracweightedcommand.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 "unifracweightedcommand.h"
  11
  12 /***********************************************************/
  13 UnifracWeightedCommand::UnifracWeightedCommand() {
  14         try {
  15                 globaldata = GlobalData::getInstance();
  16
  17                 T = globaldata->gTree;
  18                 tmap = globaldata->gTreemap;
  19                 sumFile = globaldata->getTreeFile() + ".wsummary";
  20                 openOutputFile(sumFile, outSum);
  21
  22                 util = new SharedUtil();
  23                 string s; //to make work with setgroups
  24                 util->setGroups(globaldata->Groups, tmap->namesOfGroups, s, numGroups, "weighted");     //sets the groups the user wants to analyze
  25                 util->getCombos(groupComb, globaldata->Groups, numComp);
  26                 globaldata->setGroups("");
  27
  28                 convert(globaldata->getIters(), iters);  //how many random trees to generate
  29                 weighted = new Weighted(tmap);
  30
  31         }
  32         catch(exception& e) {
  33                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracWeightedCommand class Function UnifracWeightedCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  34                 exit(1);
  35         }
  36         catch(...) {
  37                 cout << "An unknown error has occurred in the UnifracWeightedCommand class function UnifracWeightedCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  38                 exit(1);
  39         }
  40 }
  41 /***********************************************************/
  42 int UnifracWeightedCommand::execute() {
  43         try {
  44                 Progress* reading;
  45                 reading = new Progress("Comparing to random:", iters);
  46
  47                 //get weighted for users tree
  48                 userData.resize(numComp,0);  //data[0] = weightedscore AB, data[1] = weightedscore AC...
  49                 randomData.resize(numComp,0); //data[0] = weightedscore AB, data[1] = weightedscore AC...
  50
  51                 //create new tree with same num nodes and leaves as users
  52                 randT = new Tree();
  53
  54                 //get weighted scores for users trees
  55                 for (int i = 0; i < T.size(); i++) {
  56                         counter = 0;
  57                         rScores.resize(numComp);  //data[0] = weightedscore AB, data[1] = weightedscore AC...
  58                         uScores.resize(numComp);  //data[0] = weightedscore AB, data[1] = weightedscore AC...
  59
  60                         output = new ColumnFile(globaldata->getTreeFile()  + toString(i+1) + ".weighted");
  61
  62                         userData = weighted->getValues(T[i]);  //userData[0] = weightedscore
  63
  64                         //save users score
  65                         for (int s=0; s<numComp; s++) {
  66                                 //add users score to vector of user scores
  67                                 uScores[s].push_back(userData[s]);
  68
  69                                 //save users tree score for summary file
  70                                 utreeScores.push_back(userData[s]);
  71                         }
  72
  73                         //get scores for random trees
  74                         for (int j = 0; j < iters; j++) {
  75                                 int count = 0;
  76                                 for (int r=0; r<numGroups; r++) {
  77                                         for (int l = r+1; l < numGroups; l++) {
  78                                                 //copy T[i]'s info.
  79                                                 randT->getCopy(T[i]);
  80
  81                                                 //create a random tree with same topology as T[i], but different labels
  82                                                 randT->assembleRandomUnifracTree(globaldata->Groups[r], globaldata->Groups[l]);
  83                                                 //get wscore of random tree
  84                                                 randomData = weighted->getValues(randT, globaldata->Groups[r], globaldata->Groups[l]);
  85
  86                                                 //save scores
  87                                                 rScores[count].push_back(randomData[0]);
  88                                                 count++;
  89                                         }
  90                                 }
  91
  92                                 //update progress bar
  93                                 reading->update(j);
  94
  95                         }
  96
  97                         //removeValidScoresDuplicates();
  98                         //find the signifigance of the score for summary file
  99                         for (int f = 0; f < numComp; f++) {
 100                                 //sort random scores
 101                                 sort(rScores[f].begin(), rScores[f].end());
 102
 103                                 //the index of the score higher than yours is returned
 104                                 //so if you have 1000 random trees the index returned is 100
 105                                 //then there are 900 trees with a score greater then you.
 106                                 //giving you a signifigance of 0.900
 107                                 int index = findIndex(userData[f], f);    if (index == -1) { cout << "error in UnifracWeightedCommand" << endl; exit(1); } //error code
 108
 109                                 //the signifigance is the number of trees with the users score or higher
 110                                 WScoreSig.push_back((iters-index)/(float)iters);
 111                         }
 112
 113                         //out << "Tree# " << i << endl;
 114                         calculateFreqsCumuls();
 115                         printWeightedFile();
 116
 117                         delete output;
 118
 119                         //clear data
 120                         rScores.clear();
 121                         uScores.clear();
 122                         validScores.clear();
 123                 }
 124
 125                 //finish progress bar
 126                 reading->finish();
 127                 delete reading;
 128
 129                 printWSummaryFile();
 130
 131                 //clear out users groups
 132                 globaldata->Groups.clear();
 133
 134                 delete randT;
 135
 136                 return 0;
 137
 138         }
 139         catch(exception& e) {
 140                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracWeightedCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 141                 exit(1);
 142         }
 143         catch(...) {
 144                 cout << "An unknown error has occurred in the UnifracWeightedCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 145                 exit(1);
 146         }
 147 }
 148 /***********************************************************/
 149 void UnifracWeightedCommand::printWeightedFile() {
 150         try {
 151                 vector<double> data;
 152                 vector<string> tags;
 153                 tags.push_back("Score"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
 154
 155                 for(int a = 0; a < numComp; a++) {
 156                         output->initFile(groupComb[a], tags);
 157                         //print each line
 158                         for (it = validScores.begin(); it != validScores.end(); it++) {
 159                                 data.push_back(it->first);  data.push_back(rScoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
 160                                 output->output(data);
 161                                 data.clear();
 162                         }
 163                         output->resetFile();
 164                 }
 165         }
 166         catch(exception& e) {
 167                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracWeightedCommand class Function printWeightedFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 168                 exit(1);
 169         }
 170         catch(...) {
 171                 cout << "An unknown error has occurred in the UnifracWeightedCommand class function printWeightedFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 172                 exit(1);
 173         }
 174 }
 175
 176
 177 /***********************************************************/
 178 void UnifracWeightedCommand::printWSummaryFile() {
 179         try {
 180                 //column headers
 181                 outSum << "Tree#" << '\t' << "Groups" << '\t' << "WScore" << '\t' << "WSig" <<  endl;
 182                 cout << "Tree#" << '\t' << "Groups" << '\t' << "WScore" << '\t' << "WSig" <<  endl;
 183
 184                 //format output
 185                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 186
 187                 //print each line
 188                 int count = 0;
 189                 for (int i = 0; i < T.size(); i++) {
 190                         for (int j = 0; j < numComp; j++) {
 191                                 if (WScoreSig[count] > (1/(float)iters)) {
 192                                         outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(globaldata->getIters().length()) << WScoreSig[count] << endl;
 193                                         cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(globaldata->getIters().length()) << WScoreSig[count] << endl;
 194                                 }else{
 195                                         outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(globaldata->getIters().length()) << "<" << (1/float(iters)) << endl;
 196                                         cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << setprecision(globaldata->getIters().length()) << "<" << (1/float(iters)) << endl;
 197                                 }
 198                                 count++;
 199                         }
 200                 }
 201                 outSum.close();
 202         }
 203         catch(exception& e) {
 204                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracWeightedCommand class Function printWeightedFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 205                 exit(1);
 206         }
 207         catch(...) {
 208                 cout << "An unknown error has occurred in the UnifracWeightedCommand class function printWeightedFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 209                 exit(1);
 210         }
 211 }
 212
 213 /***********************************************************/
 214 int UnifracWeightedCommand::findIndex(float score, int index) {
 215         try{
 216                 for (int i = 0; i < rScores[index].size(); i++) {
 217                         if (rScores[index][i] >= score) {       return i;       }
 218                 }
 219                 return rScores[index].size();
 220         }
 221         catch(exception& e) {
 222                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracWeightedCommand class Function findIndex. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 223                 exit(1);
 224         }
 225         catch(...) {
 226                 cout << "An unknown error has occurred in the UnifracWeightedCommand class function findIndex. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 227                 exit(1);
 228         }
 229 }
 230
 231 /***********************************************************/
 232
 233 void UnifracWeightedCommand::calculateFreqsCumuls() {
 234         try {
 235                 //clear out old tree values
 236                 rScoreFreq.clear();
 237                 rScoreFreq.resize(numComp);
 238                 rCumul.clear();
 239                 rCumul.resize(numComp);
 240                 validScores.clear();
 241
 242                 //calculate frequency
 243                 for (int f = 0; f < numComp; f++) {
 244                         for (int i = 0; i < rScores[f].size(); i++) { //looks like 0,0,1,1,1,2,4,7...  you want to make a map that say rScoreFreq[0] = 2, rScoreFreq[1] = 3...
 245                                 validScores[rScores[f][i]] = rScores[f][i];
 246                                 it = rScoreFreq[f].find(rScores[f][i]);
 247                                 if (it != rScoreFreq[f].end()) {
 248                                         rScoreFreq[f][rScores[f][i]]++;
 249                                 }else{
 250                                         rScoreFreq[f][rScores[f][i]] = 1;
 251                                 }
 252                         }
 253                 }
 254
 255                 //calculate rcumul
 256                 for(int a = 0; a < numComp; a++) {
 257                         float rcumul = 1.0000;
 258                         //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
 259                         for (it = validScores.begin(); it != validScores.end(); it++) {
 260                                 //make rscoreFreq map and rCumul
 261                                 it2 = rScoreFreq[a].find(it->first);
 262                                 rCumul[a][it->first] = rcumul;
 263                                 //get percentage of random trees with that info
 264                                 if (it2 != rScoreFreq[a].end()) {  rScoreFreq[a][it->first] /= iters; rcumul-= it2->second;  }
 265                                 else { rScoreFreq[a][it->first] = 0.0000; } //no random trees with that score
 266                         }
 267                 }
 268
 269         }
 270         catch(exception& e) {
 271                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracWeightedCommand class Function calculateFreqsCums. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 272                 exit(1);
 273         }
 274         catch(...) {
 275                 cout << "An unknown error has occurred in the UnifracWeightedCommand class function calculateFreqsCums. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 276                 exit(1);
 277         }
 278
 279 }
 280
 281 /***********************************************************/
 282
 283
 284
 285
 286