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 = 1.0000;
 117                         for (it = rscoreFreq.begin(); it != rscoreFreq.end(); it++) {
 118                                 rCumul[it->first] = rcumul;
 119                                 //get percentage of random trees with that info
 120                                 rscoreFreq[it->first] /= iters;
 121                                 rcumul-= it->second;
 122
 123                         }
 124
 125                         //save the signifigance of the users score for printing later
 126                         UWScoreSig.push_back(rCumul[userData[0]]);
 127
 128
 129                         //clear random data
 130                         rscoreFreq.clear();  //you clear this because in the summary file you want the unweighted signifinance to be relative to these 1000 trees.
 131                         rCumul.clear();
 132                 }
 133
 134                 float ucumul = 1.0000;
 135                 float rcumul = 1.0000;
 136                 //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
 137                 for (it = validScores.begin(); it != validScores.end(); it++) {
 138                         it2 = uscoreFreq.find(it->first);
 139                         //make uCumul map
 140                         uCumul[it->first] = ucumul;
 141                         //user data has that score
 142                         if (it2 != uscoreFreq.end()) { uscoreFreq[it->first] /= T.size(); ucumul-= it2->second;  }
 143                         else { uscoreFreq[it->first] = 0.0000; } //no user trees with that score
 144
 145                         //make rscoreFreq map and rCumul
 146                         it2 = totalrscoreFreq.find(it->first);
 147                         rCumul[it->first] = rcumul;
 148                         //get percentage of random trees with that info
 149                         if (it2 != totalrscoreFreq.end()) {  totalrscoreFreq[it->first] /= (iters*T.size()); rcumul-= it2->second;  }
 150                         else { totalrscoreFreq[it->first] = 0.0000; } //no random trees with that score
 151
 152                 }
 153
 154                 printUnweightedFile();
 155                 printUWSummaryFile();
 156
 157                 delete randT;
 158
 159                 return 0;
 160
 161         }
 162         catch(exception& e) {
 163                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 164                 exit(1);
 165         }
 166         catch(...) {
 167                 cout << "An unknown error has occurred in the UnifracUnweightedCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 168                 exit(1);
 169         }
 170 }
 171 /***********************************************************/
 172 void UnifracUnweightedCommand::printUnweightedFile() {
 173         try {
 174                 //column headers
 175
 176                 out << "Score" << '\t' << "UserFreq" << '\t' << "UserCumul" << '\t' << "RandFreq" << '\t' << "RandCumul" << endl;
 177
 178                 //format output
 179                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 180
 181                 //print each line
 182                 for (it = validScores.begin(); it != validScores.end(); it++) {
 183                         out << setprecision(6) << it->first << '\t' << '\t' << uscoreFreq[it->first] << '\t' << uCumul[it->first] << '\t' << totalrscoreFreq[it->first] << '\t' << rCumul[it->first] << endl;
 184                 }
 185
 186                 out.close();
 187
 188         }
 189         catch(exception& e) {
 190                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function printUnweightedFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 191                 exit(1);
 192         }
 193         catch(...) {
 194                 cout << "An unknown error has occurred in the UnifracUnweightedCommand class function printUnweightedFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 195                 exit(1);
 196         }
 197 }
 198
 199 /***********************************************************/
 200 void UnifracUnweightedCommand::printUWSummaryFile() {
 201         try {
 202                 //column headers
 203                 outSum << "Tree#" << '\t'  <<  "UWScore" << '\t' << '\t' << "UWSig" <<  endl;
 204
 205                 //format output
 206                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 207
 208                 //print each line
 209                 for (int i = 0; i< T.size(); i++) {
 210                         outSum << setprecision(6) << i+1 << '\t' << '\t' << utreeScores[i] << '\t' << UWScoreSig[i] << endl;
 211                 }
 212
 213                 outSum.close();
 214         }
 215         catch(exception& e) {
 216                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function printUWSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 217                 exit(1);
 218         }
 219         catch(...) {
 220                 cout << "An unknown error has occurred in the UnifracUnweightedCommand class function printUWSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 221                 exit(1);
 222         }
 223 }
 224 /***********************************************************/
 225 void UnifracUnweightedCommand::saveRandomScores() {
 226         try {
 227                 for (it = rscoreFreq.begin(); it != rscoreFreq.end(); it++) {
 228                         //does this score already exist in the total map
 229                         it2 = totalrscoreFreq.find(it->first);
 230                         //if yes then add them
 231                         if (it2 != totalrscoreFreq.end()) {
 232                                 totalrscoreFreq[it->first] += rscoreFreq[it->first];
 233                         }else{ //its a new score
 234                                 totalrscoreFreq[it->first] = rscoreFreq[it->first];
 235                         }
 236                 }
 237         }
 238         catch(exception& e) {
 239                 cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function saveRandomScores. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 240                 exit(1);
 241         }
 242         catch(...) {
 243                 cout << "An unknown error has occurred in the UnifracUnweightedCommand class function saveRandomScores. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 244                 exit(1);
 245         }
 246 }
 247
 248 /***********************************************************/