parsimonycommand.cpp

   1 /*
   2  *  parsimonycommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 1/26/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "parsimonycommand.h"
  11
  12 //**********************************************************************************************************************
  13 vector<string> ParsimonyCommand::getValidParameters(){
  14         try {
  15                 string Array[] =  {"random","groups","iters","outputdir","inputdir"};
  16                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  17                 return myArray;
  18         }
  19         catch(exception& e) {
  20                 m->errorOut(e, "ParsimonyCommand", "getValidParameters");
  21                 exit(1);
  22         }
  23 }
  24 //**********************************************************************************************************************
  25 ParsimonyCommand::ParsimonyCommand(){
  26         try {
  27                 abort = true;
  28                 //initialize outputTypes
  29                 vector<string> tempOutNames;
  30                 outputTypes["parsimony"] = tempOutNames;
  31                 outputTypes["psummary"] = tempOutNames;
  32         }
  33         catch(exception& e) {
  34                 m->errorOut(e, "ParsimonyCommand", "ParsimonyCommand");
  35                 exit(1);
  36         }
  37 }
  38 //**********************************************************************************************************************
  39 vector<string> ParsimonyCommand::getRequiredParameters(){
  40         try {
  41                 vector<string> myArray;
  42                 return myArray;
  43         }
  44         catch(exception& e) {
  45                 m->errorOut(e, "ParsimonyCommand", "getRequiredParameters");
  46                 exit(1);
  47         }
  48 }
  49 //**********************************************************************************************************************
  50 vector<string> ParsimonyCommand::getRequiredFiles(){
  51         try {
  52                 vector<string> myArray;
  53                 return myArray;
  54         }
  55         catch(exception& e) {
  56                 m->errorOut(e, "ParsimonyCommand", "getRequiredFiles");
  57                 exit(1);
  58         }
  59 }
  60 /***********************************************************/
  61 ParsimonyCommand::ParsimonyCommand(string option)  {
  62         try {
  63                 globaldata = GlobalData::getInstance();
  64                 abort = false;
  65                 Groups.clear();
  66
  67                 //allow user to run help
  68                 if(option == "help") { help(); abort = true; }
  69
  70                 else {
  71                         //valid paramters for this command
  72                         string Array[] =  {"random","groups","iters","outputdir","inputdir"};
  73                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  74
  75                         OptionParser parser(option);
  76                         map<string, string> parameters = parser.getParameters();
  77
  78                         ValidParameters validParameter;
  79
  80                         //check to make sure all parameters are valid for command
  81                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  82                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  83                         }
  84
  85                         //initialize outputTypes
  86                         vector<string> tempOutNames;
  87                         outputTypes["parsimony"] = tempOutNames;
  88                         outputTypes["psummary"] = tempOutNames;
  89
  90                         randomtree = validParameter.validFile(parameters, "random", false);             if (randomtree == "not found") { randomtree = ""; }
  91
  92                         //are you trying to use parsimony without reading a tree or saying you want random distribution
  93                         if (randomtree == "")  {
  94                                 if (globaldata->gTree.size() == 0) {
  95                                         m->mothurOut("You must read a treefile and a groupfile or set the randomtree parameter to the output filename you wish, before you may execute the parsimony command."); m->mothurOutEndLine(); abort = true;  }
  96                         }
  97
  98                         //if the user changes the output directory command factory will send this info to us in the output parameter
  99                         string outputDir = validParameter.validFile(parameters, "outputdir", false);            if (outputDir == "not found"){  outputDir = ""; }
 100
 101                         //check for optional parameter and set defaults
 102                         // ...at some point should added some additional type checking...
 103                         groups = validParameter.validFile(parameters, "groups", false);
 104                         if (groups == "not found") { groups = ""; globaldata->Groups.clear(); }
 105                         else {
 106                                 m->splitAtDash(groups, Groups);
 107                                 globaldata->Groups = Groups;
 108                         }
 109
 110                         itersString = validParameter.validFile(parameters, "iters", false);                     if (itersString == "not found") { itersString = "1000"; }
 111                         convert(itersString, iters);
 112
 113                         if (abort == false) {
 114                                 //randomtree will tell us if user had their own treefile or if they just want the random distribution
 115                                 //user has entered their own tree
 116                                 if (randomtree == "") {
 117                                         T = globaldata->gTree;
 118                                         tmap = globaldata->gTreemap;
 119
 120                                         if(outputDir == "") { outputDir += m->hasPath(globaldata->getTreeFile()); }
 121                                         output = new ColumnFile(outputDir + m->getSimpleName(globaldata->getTreeFile())  +  ".parsimony", itersString);
 122                                         outputNames.push_back(outputDir + m->getSimpleName(globaldata->getTreeFile())  +  ".parsimony");
 123                                         outputTypes["parsimony"].push_back(outputDir + m->getSimpleName(globaldata->getTreeFile())  +  ".parsimony");
 124
 125                                         sumFile = outputDir + m->getSimpleName(globaldata->getTreeFile()) + ".psummary";
 126                                         m->openOutputFile(sumFile, outSum);
 127                                         outputNames.push_back(sumFile);
 128                                         outputTypes["psummary"].push_back(sumFile);
 129                                 }else { //user wants random distribution
 130                                         savetmap = globaldata->gTreemap;
 131                                         getUserInput();
 132
 133                                         if(outputDir == "") { outputDir += m->hasPath(randomtree); }
 134                                         output = new ColumnFile(outputDir+ m->getSimpleName(randomtree), itersString);
 135                                         outputNames.push_back(outputDir+ m->getSimpleName(randomtree));
 136                                         outputTypes["parsimony"].push_back(outputDir+ m->getSimpleName(randomtree));
 137                                 }
 138
 139                                 //set users groups to analyze
 140                                 util = new SharedUtil();
 141                                 util->setGroups(globaldata->Groups, tmap->namesOfGroups, allGroups, numGroups, "parsimony");    //sets the groups the user wants to analyze
 142                                 util->getCombos(groupComb, globaldata->Groups, numComp);
 143
 144                                 if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
 145
 146                                 pars = new Parsimony(tmap);
 147                                 counter = 0;
 148
 149                         }
 150
 151                 }
 152
 153         }
 154         catch(exception& e) {
 155                 m->errorOut(e, "ParsimonyCommand", "ParsimonyCommand");
 156                 exit(1);
 157         }
 158 }
 159
 160 //**********************************************************************************************************************
 161
 162 void ParsimonyCommand::help(){
 163         try {
 164                 m->mothurOut("The parsimony command can only be executed after a successful read.tree command, unless you use the random parameter.\n");
 165                 m->mothurOut("The parsimony command parameters are random, groups and iters.  No parameters are required.\n");
 166                 m->mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed.  You must enter at least 1 valid group.\n");
 167                 m->mothurOut("The group names are separated by dashes.  The iters parameter allows you to specify how many random trees you would like compared to your tree.\n");
 168                 m->mothurOut("The parsimony command should be in the following format: parsimony(random=yourOutputFilename, groups=yourGroups, iters=yourIters).\n");
 169                 m->mothurOut("Example parsimony(random=out, iters=500).\n");
 170                 m->mothurOut("The default value for random is "" (meaning you want to use the trees in your inputfile, randomtree=out means you just want the random distribution of trees outputted to out.rd_parsimony),\n");
 171                 m->mothurOut("and iters is 1000.  The parsimony command output two files: .parsimony and .psummary their descriptions are in the manual.\n");
 172                 m->mothurOut("Note: No spaces between parameter labels (i.e. random), '=' and parameters (i.e.yourOutputFilename).\n\n");
 173         }
 174         catch(exception& e) {
 175                 m->errorOut(e, "ParsimonyCommand", "help");
 176                 exit(1);
 177         }
 178 }
 179
 180
 181 /***********************************************************/
 182 int ParsimonyCommand::execute() {
 183         try {
 184
 185                 if (abort == true) { return 0; }
 186
 187                 Progress* reading;
 188                 reading = new Progress("Comparing to random:", iters);
 189
 190                 if (m->control_pressed) {
 191                         delete reading; delete pars; delete util; delete output;
 192                         if (randomtree == "") {  outSum.close();  }
 193                         for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } outputTypes.clear();
 194                         globaldata->Groups.clear();
 195                         return 0;
 196                 }
 197
 198
 199                 //get pscore for users tree
 200                 userData.resize(numComp,0);  //data = AB, AC, BC, ABC.
 201                 randomData.resize(numComp,0);  //data = AB, AC, BC, ABC.
 202                 rscoreFreq.resize(numComp);
 203                 uscoreFreq.resize(numComp);
 204                 rCumul.resize(numComp);
 205                 uCumul.resize(numComp);
 206                 userTreeScores.resize(numComp);
 207                 UScoreSig.resize(numComp);
 208
 209                 if (randomtree == "") {
 210                         //get pscores for users trees
 211                         for (int i = 0; i < T.size(); i++) {
 212                                 userData = pars->getValues(T[i]);  //data = AB, AC, BC, ABC.
 213
 214                                 if (m->control_pressed) {
 215                                         delete reading; delete pars; delete util; delete output;
 216                                         if (randomtree == "") {  outSum.close();  }
 217                                         for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } outputTypes.clear();
 218                                         globaldata->Groups.clear();
 219                                         return 0;
 220                                 }
 221
 222
 223                                 //output scores for each combination
 224                                 for(int k = 0; k < numComp; k++) {
 225
 226                                         //update uscoreFreq
 227                                         map<int,double>::iterator it = uscoreFreq[k].find(userData[k]);
 228                                         if (it == uscoreFreq[k].end()) {//new score
 229                                                 uscoreFreq[k][userData[k]] = 1;
 230                                         }else{ uscoreFreq[k][userData[k]]++; }
 231
 232                                         //add users score to valid scores
 233                                         validScores[userData[k]] = userData[k];
 234
 235                                         //save score for summary file
 236                                         userTreeScores[k].push_back(userData[k]);
 237                                 }
 238                         }
 239
 240                         //get pscores for random trees
 241                         for (int j = 0; j < iters; j++) {
 242
 243                                 //create new tree with same num nodes and leaves as users
 244                                 randT = new Tree();
 245
 246                                 //create random relationships between nodes
 247                                 randT->assembleRandomTree();
 248
 249                                 //get pscore of random tree
 250                                 randomData = pars->getValues(randT);
 251
 252                                 if (m->control_pressed) {
 253                                         delete reading; delete pars; delete util; delete output; delete randT;
 254                                         if (randomtree == "") {  outSum.close();  }
 255                                         for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } outputTypes.clear();
 256                                         globaldata->Groups.clear();
 257                                         return 0;
 258                                 }
 259
 260                                 for(int r = 0; r < numComp; r++) {
 261                                         //add trees pscore to map of scores
 262                                         map<int,double>::iterator it = rscoreFreq[r].find(randomData[r]);
 263                                         if (it != rscoreFreq[r].end()) {//already have that score
 264                                                 rscoreFreq[r][randomData[r]]++;
 265                                         }else{//first time we have seen this score
 266                                                 rscoreFreq[r][randomData[r]] = 1;
 267                                         }
 268
 269                                         //add randoms score to validscores
 270                                         validScores[randomData[r]] = randomData[r];
 271                                 }
 272
 273                                 //update progress bar
 274                                 reading->update(j);
 275
 276                                 delete randT;
 277                         }
 278
 279                 }else {
 280                         //get pscores for random trees
 281                         for (int j = 0; j < iters; j++) {
 282
 283                                 //create new tree with same num nodes and leaves as users
 284                                 randT = new Tree();
 285                                 //create random relationships between nodes
 286
 287                                 randT->assembleRandomTree();
 288
 289                                 if (m->control_pressed) {
 290                                         delete reading; delete pars; delete util; delete output; delete randT;
 291                                         globaldata->gTreemap = savetmap;
 292                                         for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } outputTypes.clear();
 293                                         globaldata->Groups.clear();
 294                                         return 0;
 295                                 }
 296
 297
 298                                 //get pscore of random tree
 299                                 randomData = pars->getValues(randT);
 300
 301                                 if (m->control_pressed) {
 302                                         delete reading; delete pars; delete util; delete output; delete randT;
 303                                         globaldata->gTreemap = savetmap;
 304                                         for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } outputTypes.clear();
 305                                         globaldata->Groups.clear();
 306                                         return 0;
 307                                 }
 308
 309                                 for(int r = 0; r < numComp; r++) {
 310                                         //add trees pscore to map of scores
 311                                         map<int,double>::iterator it = rscoreFreq[r].find(randomData[r]);
 312                                         if (it != rscoreFreq[r].end()) {//already have that score
 313                                                 rscoreFreq[r][randomData[r]]++;
 314                                         }else{//first time we have seen this score
 315                                                 rscoreFreq[r][randomData[r]] = 1;
 316                                         }
 317
 318                                         //add randoms score to validscores
 319                                         validScores[randomData[r]] = randomData[r];
 320                                 }
 321
 322                                 //update progress bar
 323                                 reading->update(j);
 324
 325                                 delete randT;
 326                         }
 327                 }
 328
 329                 for(int a = 0; a < numComp; a++) {
 330                         float rcumul = 0.0000;
 331                         float ucumul = 0.0000;
 332                         //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
 333                         for (map<int,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 334                                 if (randomtree == "") {
 335                                         map<int,double>::iterator it2 = uscoreFreq[a].find(it->first);
 336                                         //user data has that score
 337                                         if (it2 != uscoreFreq[a].end()) { uscoreFreq[a][it->first] /= T.size(); ucumul+= it2->second;  }
 338                                         else { uscoreFreq[a][it->first] = 0.0000; } //no user trees with that score
 339                                         //make uCumul map
 340                                         uCumul[a][it->first] = ucumul;
 341                                 }
 342
 343                                 //make rscoreFreq map and rCumul
 344                                 map<int,double>::iterator it2 = rscoreFreq[a].find(it->first);
 345                                 //get percentage of random trees with that info
 346                                 if (it2 != rscoreFreq[a].end()) {  rscoreFreq[a][it->first] /= iters; rcumul+= it2->second;  }
 347                                 else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
 348                                 rCumul[a][it->first] = rcumul;
 349                         }
 350
 351                         //find the signifigance of each user trees score when compared to the random trees and save for printing the summary file
 352                         for (int h = 0; h < userTreeScores[a].size(); h++) {
 353                                 UScoreSig[a].push_back(rCumul[a][userTreeScores[a][h]]);
 354                         }
 355                 }
 356
 357                 if (m->control_pressed) {
 358                                 delete reading; delete pars; delete util; delete output;
 359                                 if (randomtree == "") {  outSum.close();  }
 360                                 else { globaldata->gTreemap = savetmap; }
 361                                 for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } outputTypes.clear();
 362                                 globaldata->Groups.clear();
 363                                 return 0;
 364                 }
 365
 366                 //finish progress bar
 367                 reading->finish();
 368                 delete reading;
 369
 370
 371                 printParsimonyFile();
 372                 if (randomtree == "") { printUSummaryFile(); }
 373
 374                 //reset globaldata's treemap if you just did random distrib
 375                 if (randomtree != "") {
 376                         //memory leak prevention
 377                         //if (globaldata->gTreemap != NULL) { delete globaldata->gTreemap;  }
 378                         globaldata->gTreemap = savetmap;
 379                 }
 380
 381                 //reset groups parameter
 382                 globaldata->Groups.clear();
 383
 384                 if (m->control_pressed) {
 385                         delete pars; delete util; delete output;
 386                         for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } outputTypes.clear();
 387                         return 0;
 388                 }
 389
 390                 m->mothurOutEndLine();
 391                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 392                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 393                 m->mothurOutEndLine();
 394
 395
 396                 return 0;
 397
 398         }
 399         catch(exception& e) {
 400                 m->errorOut(e, "ParsimonyCommand", "execute");
 401                 exit(1);
 402         }
 403 }
 404
 405 /***********************************************************/
 406 void ParsimonyCommand::printParsimonyFile() {
 407         try {
 408                 vector<double> data;
 409                 vector<string> tags;
 410
 411                 if (randomtree == "") {
 412                         tags.push_back("Score"); tags.push_back("UserFreq"); tags.push_back("UserCumul"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
 413                 }else {
 414                         tags.push_back("Score"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
 415                 }
 416
 417                 for(int a = 0; a < numComp; a++) {
 418                         output->initFile(groupComb[a], tags);
 419                         //print each line
 420                         for (map<int,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 421                                 if (randomtree == "") {
 422                                         data.push_back(it->first);  data.push_back(uscoreFreq[a][it->first]); data.push_back(uCumul[a][it->first]); data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
 423                                 }else{
 424                                         data.push_back(it->first);  data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
 425                                 }
 426                                 output->output(data);
 427                                 data.clear();
 428                         }
 429                         output->resetFile();
 430                 }
 431         }
 432         catch(exception& e) {
 433                 m->errorOut(e, "ParsimonyCommand", "printParsimonyFile");
 434                 exit(1);
 435         }
 436 }
 437 /***********************************************************/
 438 int ParsimonyCommand::printUSummaryFile() {
 439         try {
 440                 //column headers
 441                 outSum << "Tree#" << '\t' << "Groups" << '\t'  <<  "ParsScore" << '\t' << "ParsSig" <<  endl;
 442                 m->mothurOut("Tree#\tGroups\tParsScore\tParsSig"); m->mothurOutEndLine();
 443
 444                 //format output
 445                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 446
 447
 448                 //print each line
 449                 for (int i = 0; i< T.size(); i++) {
 450                         for(int a = 0; a < numComp; a++) {
 451                                 if (m->control_pressed) {  outSum.close(); return 0; }
 452                                 if (UScoreSig[a][i] > (1/(float)iters)) {
 453                                         outSum << setprecision(6) << i+1 << '\t' << groupComb[a]  << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << UScoreSig[a][i] << endl;
 454                                         cout << setprecision(6) << i+1 << '\t' << groupComb[a]  << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << UScoreSig[a][i] << endl;
 455                                         m->mothurOutJustToLog(toString(i+1) + "\t" + groupComb[a] + "\t" + toString(userTreeScores[a][i]) + "\t" + toString(UScoreSig[a][i])); m->mothurOutEndLine();
 456                                 }else {
 457                                         outSum << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length())  << '\t' << "<" << (1/float(iters)) << endl;
 458                                         cout << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << "<" << (1/float(iters)) << endl;
 459                                         m->mothurOutJustToLog(toString(i+1) + "\t" + groupComb[a] + "\t" + toString(userTreeScores[a][i]) + "\t" + toString((1/float(iters)))); m->mothurOutEndLine();
 460                                 }
 461                         }
 462                 }
 463
 464                 outSum.close();
 465                 return 0;
 466         }
 467         catch(exception& e) {
 468                 m->errorOut(e, "ParsimonyCommand", "printUSummaryFile");
 469                 exit(1);
 470         }
 471 }
 472
 473 /***********************************************************/
 474 void ParsimonyCommand::getUserInput() {
 475         try {
 476
 477                 //create treemap
 478                 tmap = new TreeMap();
 479
 480                 m->mothurOut("Please enter the number of groups you would like to analyze: ");
 481                 cin >> numGroups;
 482                 m->mothurOutJustToLog(toString(numGroups)); m->mothurOutEndLine();
 483
 484                 int num, count;
 485                 count = 1;
 486                 numEachGroup.resize(numGroups, 0);
 487
 488                 for (int i = 1; i <= numGroups; i++) {
 489                         m->mothurOut("Please enter the number of sequences in group " + toString(i) +  ": ");
 490                         cin >> num;
 491                         m->mothurOutJustToLog(toString(num)); m->mothurOutEndLine();
 492
 493                         //set tmaps seqsPerGroup
 494                         tmap->seqsPerGroup[toString(i)] = num;
 495                         tmap->namesOfGroups.push_back(toString(i));
 496
 497                         //set tmaps namesOfSeqs
 498                         for (int j = 0; j < num; j++) {
 499                                 tmap->namesOfSeqs.push_back(toString(count));
 500                                 tmap->treemap[toString(count)].groupname = toString(i);
 501                                 count++;
 502                         }
 503                 }
 504
 505                 //clears buffer so next command doesn't have error
 506                 string s;
 507                 getline(cin, s);
 508
 509                 //save tmap for later
 510                 //memory leak prevention
 511                 //if (globaldata->gTreemap != NULL) { delete globaldata->gTreemap;  }
 512                 globaldata->gTreemap = tmap;
 513                 globaldata->Treenames = tmap->namesOfSeqs;
 514
 515         }
 516         catch(exception& e) {
 517                 m->errorOut(e, "ParsimonyCommand", "getUserInput");
 518                 exit(1);
 519         }
 520 }
 521
 522 /***********************************************************/
 523
 524