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