getcoremicrobiomecommand.cpp

   1 //
   2 //  GetCoreMicroBiomeCommand.cpp
   3 //  Mothur
   4 //
   5 //  Created by John Westcott on 5/8/12.
   6 //  Copyright (c) 2012 Schloss Lab. All rights reserved.
   7 //
   8
   9 #include "GetCoreMicroBiomeCommand.h"
  10
  11
  12 //**********************************************************************************************************************
  13 vector<string> GetCoreMicroBiomeCommand::setParameters(){
  14         try {
  15         CommandParameter pshared("shared", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(pshared);
  16                 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(prelabund);
  17         CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
  18                 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
  19                 CommandParameter poutput("output", "Multiple", "fraction-count", "fraction", "", "", "",false,false); parameters.push_back(poutput);
  20         CommandParameter pabund("abundance", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pabund);
  21                 CommandParameter psamples("samples", "Number", "", "-1", "", "", "",false,false); parameters.push_back(psamples);
  22                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  23                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  24
  25                 vector<string> myArray;
  26                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  27                 return myArray;
  28         }
  29         catch(exception& e) {
  30                 m->errorOut(e, "GetCoreMicroBiomeCommand", "setParameters");
  31                 exit(1);
  32         }
  33 }
  34 //**********************************************************************************************************************
  35 string GetCoreMicroBiomeCommand::getHelpString(){
  36         try {
  37                 string helpString = "";
  38                 helpString += "The get.coremicrobiome determines the fraction of OTUs that are found in varying numbers of samples for different minimum relative abundances.\n";
  39                 helpString += "The get.coremicrobiome parameters are: shared, relabund, groups, label, output, abundance and samples. Shared or relabund is required.\n";
  40                 helpString += "The label parameter is used to analyze specific labels in your input.\n";
  41                 helpString += "The groups parameter allows you to specify which of the groups you would like analyzed.\n";
  42         helpString += "The output parameter is used to specify whether you would like the fraction of OTU's or OTU count outputted. Options are fraction or count. Default=fraction.\n";
  43                 helpString += "The abundance parameter allows you to specify an abundance you would like the OTU names outputted for. Must be an integer between 0 and 100, indicating the relative abundance. \n";
  44         helpString += "The samples parameter allows you to specify the minimum number of samples you would like the OTU names outputted for. Must be an interger between 1 and number of samples in your file.\n";
  45                 helpString += "The new command should be in the following format: get.coremicrobiome(shared=yourSharedFile)\n";
  46                 helpString += "get.coremicrobiom(shared=final.an.shared, abund=30)\n";
  47                 return helpString;
  48         }
  49         catch(exception& e) {
  50                 m->errorOut(e, "GetCoreMicroBiomeCommand", "getHelpString");
  51                 exit(1);
  52         }
  53 }
  54 //**********************************************************************************************************************
  55 GetCoreMicroBiomeCommand::GetCoreMicroBiomeCommand(){
  56         try {
  57                 abort = true; calledHelp = true;
  58                 setParameters();
  59         vector<string> tempOutNames;
  60                 outputTypes["coremicrobiom"] = tempOutNames;
  61         }
  62         catch(exception& e) {
  63                 m->errorOut(e, "GetCoreMicroBiomeCommand", "GetCoreMicroBiomeCommand");
  64                 exit(1);
  65         }
  66 }
  67 //**********************************************************************************************************************
  68 GetCoreMicroBiomeCommand::GetCoreMicroBiomeCommand(string option)  {
  69         try {
  70                 abort = false; calledHelp = false;   allLines = 1;
  71
  72                 //allow user to run help
  73                 if(option == "help") { help(); abort = true; calledHelp = true; }
  74                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  75
  76                 else {
  77                         //valid paramters for this command
  78                         vector<string> myArray = setParameters();
  79
  80                         OptionParser parser(option);
  81                         map<string,string> parameters = parser.getParameters();
  82
  83                         ValidParameters validParameter;
  84                         map<string,string>::iterator it;
  85                         //check to make sure all parameters are valid for command
  86                         for (it = parameters.begin(); it != parameters.end(); it++) {
  87                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  88                         }
  89
  90
  91                         //if the user changes the input directory command factory will send this info to us in the output parameter
  92                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  93                         if (inputDir == "not found"){   inputDir = "";          }
  94                         else {
  95
  96                                 string path;
  97                                 it = parameters.find("relabund");
  98                                 //user has given a template file
  99                                 if(it != parameters.end()){
 100                                         path = m->hasPath(it->second);
 101                                         //if the user has not given a path then, add inputdir. else leave path alone.
 102                                         if (path == "") {       parameters["relabund"] = inputDir + it->second;         }
 103                                 }
 104
 105                 it = parameters.find("shared");
 106                                 //user has given a template file
 107                                 if(it != parameters.end()){
 108                                         path = m->hasPath(it->second);
 109                                         //if the user has not given a path then, add inputdir. else leave path alone.
 110                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 111                                 }
 112             }
 113
 114
 115                         //check for parameters
 116             sharedfile = validParameter.validFile(parameters, "shared", true);
 117                         if (sharedfile == "not open") { abort = true; }
 118                         else if (sharedfile == "not found") { sharedfile = ""; }
 119                         else { inputFileName = sharedfile; format = "sharedfile"; m->setSharedFile(sharedfile); }
 120
 121                         relabundfile = validParameter.validFile(parameters, "relabund", true);
 122                         if (relabundfile == "not open") { abort = true; }
 123                         else if (relabundfile == "not found") { relabundfile = ""; }
 124                         else { inputFileName = relabundfile; format = "relabund"; m->setRelAbundFile(relabundfile); }
 125
 126             if ((relabundfile == "") && (sharedfile == "")) {
 127                                 //is there are current file available for either of these?
 128                                 //give priority to shared, then relabund
 129                                 sharedfile = m->getSharedFile();
 130                                 if (sharedfile != "") {  inputFileName = sharedfile; format="sharedfile"; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 131                                 else {
 132                                         relabundfile = m->getRelAbundFile();
 133                                         if (relabundfile != "") {  inputFileName = relabundfile; format="relabund"; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
 134                                         else {
 135                                                 m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine();
 136                                                 abort = true;
 137                                         }
 138                                 }
 139                         }
 140
 141             //if the user changes the output directory command factory will send this info to us in the output parameter
 142                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 143                                 outputDir = m->hasPath(inputFileName); //if user entered a file with a path then preserve it
 144                         }
 145
 146             string groups = validParameter.validFile(parameters, "groups", false);
 147                         if (groups == "not found") { groups = ""; }
 148                         else { m->splitAtDash(groups, Groups); }
 149                         m->setGroups(Groups);
 150
 151             string label = validParameter.validFile(parameters, "label", false);
 152                         if (label == "not found") { label = ""; }
 153                         else {
 154                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 155                                 else { allLines = 1;  }
 156                         }
 157
 158             output = validParameter.validFile(parameters, "output", false);             if(output == "not found"){      output = "fraction"; }
 159
 160                         if ((output != "fraction") && (output != "count")) { m->mothurOut(output + " is not a valid output form. Options are fraction and count. I will use fraction."); m->mothurOutEndLine(); output = "fraction"; }
 161
 162             string temp = validParameter.validFile(parameters, "abundance", false);     if (temp == "not found"){       temp = "-1";    }
 163                         m->mothurConvert(temp, abund);
 164
 165             if (abund != -1) {
 166                 if ((abund < 0) || (abund > 100)) { m->mothurOut(toString(abund) + " is not a valid number for abund. Must be an integer between 0 and 100.\n"); }
 167             }
 168
 169             temp = validParameter.validFile(parameters, "samples", false);      if (temp == "not found"){       temp = "-1";    }
 170                         m->mothurConvert(temp, samples);
 171
 172                 }
 173
 174         }
 175         catch(exception& e) {
 176                 m->errorOut(e, "GetCoreMicroBiomeCommand", "GetCoreMicroBiomeCommand");
 177                 exit(1);
 178         }
 179 }
 180 //**********************************************************************************************************************
 181
 182 int GetCoreMicroBiomeCommand::execute(){
 183         try {
 184
 185                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 186
 187         InputData input(inputFileName, format);
 188         vector<SharedRAbundFloatVector*> lookup = input.getSharedRAbundFloatVectors();
 189         string lastLabel = lookup[0]->getLabel();
 190
 191         if (samples != -1) {
 192             if ((samples < 1) || (samples > lookup.size())) { m->mothurOut(toString(samples) + " is not a valid number for samples. Must be an integer between 1 and the number of samples in your file. Your file contains " + toString(lookup.size()) + " samples, so I will use that.\n"); samples = lookup.size(); }
 193         }
 194
 195
 196         //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 197         set<string> processedLabels;
 198         set<string> userLabels = labels;
 199
 200         //as long as you are not at the end of the file or done wih the lines you want
 201         while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 202
 203             if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }  for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
 204
 205             if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 206
 207                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 208
 209                 createTable(lookup);
 210
 211                 processedLabels.insert(lookup[0]->getLabel());
 212                 userLabels.erase(lookup[0]->getLabel());
 213             }
 214
 215             if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 216                 string saveLabel = lookup[0]->getLabel();
 217
 218                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 219                 lookup = input.getSharedRAbundFloatVectors(lastLabel);
 220                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 221
 222                 createTable(lookup);
 223
 224                 processedLabels.insert(lookup[0]->getLabel());
 225                 userLabels.erase(lookup[0]->getLabel());
 226
 227                 //restore real lastlabel to save below
 228                 lookup[0]->setLabel(saveLabel);
 229             }
 230
 231             lastLabel = lookup[0]->getLabel();
 232             //prevent memory leak
 233             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 234
 235             if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); }  return 0; }
 236
 237             //get next line to process
 238             lookup = input.getSharedRAbundFloatVectors();
 239         }
 240
 241         if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); }  return 0; }
 242
 243         //output error messages about any remaining user labels
 244         set<string>::iterator it;
 245         bool needToRun = false;
 246         for (it = userLabels.begin(); it != userLabels.end(); it++) {
 247             m->mothurOut("Your file does not include the label " + *it);
 248             if (processedLabels.count(lastLabel) != 1) {
 249                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 250                 needToRun = true;
 251             }else {
 252                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 253             }
 254         }
 255
 256         //run last label if you need to
 257         if (needToRun == true)  {
 258             for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
 259             lookup = input.getSharedRAbundFloatVectors(lastLabel);
 260
 261             m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 262
 263             createTable(lookup);
 264
 265             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 266         }
 267
 268         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); }  return 0; }
 269
 270         //output files created by command
 271                 m->mothurOutEndLine();
 272                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 273                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 274                 m->mothurOutEndLine();
 275         return 0;
 276
 277     }
 278         catch(exception& e) {
 279                 m->errorOut(e, "GetCoreMicroBiomeCommand", "execute");
 280                 exit(1);
 281         }
 282 }
 283 //**********************************************************************************************************************
 284
 285 int GetCoreMicroBiomeCommand::createTable(vector<SharedRAbundFloatVector*>& lookup){
 286         try {
 287
 288         string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + ".core.microbiom";
 289         outputNames.push_back(outputFileName);  outputTypes["coremicrobiom"].push_back(outputFileName);
 290                 ofstream out;
 291                 m->openOutputFile(outputFileName, out);
 292
 293         int numSamples = lookup.size();
 294         int numOtus = lookup[0]->getNumBins();
 295
 296         //table is 100 by numsamples
 297         //question we are answering is: what fraction of OTUs in a study have a relative abundance at or above %X
 298         //in at least %Y samples. x goes from 0 to 100, y from 1 to numSamples
 299         vector< vector<double> > table; table.resize(101);
 300         for (int i = 0; i < table.size(); i++) { table[i].resize(numSamples, 0.0); }
 301
 302         map<int, vector<string> > otuNames;
 303         if ((abund != -1) && (samples == -1)) { //fill with all samples
 304             for (int i = 0; i < numSamples; i++) {
 305                 vector<string> temp;
 306                 otuNames[i+1] = temp;
 307             }
 308         }else if ((abund == -1) && (samples != -1)) { //fill with all relabund
 309             for (int i = 0; i < 101; i++) {
 310                 vector<string> temp;
 311                 otuNames[i] = temp;
 312             }
 313         }else if ((abund != -1) && (samples != -1)) { //only one line is wanted
 314             vector<string> temp;
 315             otuNames[abund] = temp;
 316         }
 317
 318         for (int i = 0; i < numOtus; i++) {
 319
 320             if (m->control_pressed) { break; }
 321
 322             //count number of samples in this otu with a relabund >= spot in count
 323             vector<int> counts; counts.resize(101, 0);
 324
 325             for (int j = 0; j < lookup.size(); j++) {
 326                 double relabund = lookup[j]->getAbundance(i);
 327                 int wholeRelabund = (int) (floor(relabund*100));
 328                 for (int k = 0; k < wholeRelabund+1; k++) { counts[k]++; }
 329             }
 330
 331             //add this otus info to table
 332             for (int j = 0; j < table.size(); j++) {
 333                 for (int k = 0; k < counts[j]; k++) { table[j][k]++; }
 334
 335                 if ((abund == -1) && (samples != -1)) { //we want all OTUs with this number of samples
 336                     if (counts[j] >= samples) { otuNames[j].push_back(m->currentBinLabels[i]); }
 337                 }else if ((abund != -1) && (samples == -1)) { //we want all OTUs with this relabund
 338                     if (j == abund) {
 339                         for (int k = 0; k < counts[j]; k++) {  otuNames[k+1].push_back(m->currentBinLabels[i]); }
 340                     }
 341                 }else if ((abund != -1) && (samples != -1)) { //we want only OTUs with this relabund for this number of samples
 342                     if ((j == abund) && (counts[j] >= samples)) {
 343                         otuNames[j].push_back(m->currentBinLabels[i]);
 344                     }
 345                 }
 346             }
 347         }
 348
 349         //format output
 350                 if (output == "fraction") { out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); }
 351         out << "NumSamples\t";
 352
 353         //convert table counts to percents
 354         for (int i = 0; i < table.size(); i++) {
 355             out << "Relabund-" << i << "%\t";
 356             if (m->control_pressed) { break; }
 357             for (int j = 0; j < table[i].size(); j++) {  if (output == "fraction") { table[i][j] /= (double) numOtus; } }
 358         }
 359         out << endl;
 360
 361         for (int i = 0; i < numSamples; i++) {
 362             if (m->control_pressed) { break; }
 363             out << i+1 << '\t';
 364             for (int j = 0; j < table.size(); j++) {  out << table[j][i] << '\t'; }
 365             out << endl;
 366         }
 367
 368         out.close();
 369
 370         if (m->control_pressed) { return 0; }
 371
 372         if ((samples != -1) || (abund != -1))  {
 373             string outputFileName2 = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + ".core.microbiomlist";
 374             outputNames.push_back(outputFileName2);  outputTypes["coremicrobiom"].push_back(outputFileName2);
 375             ofstream out2;
 376             m->openOutputFile(outputFileName2, out2);
 377
 378             if ((abund == -1) && (samples != -1)) { //we want all OTUs with this number of samples
 379                 out2 << "Relabund%\tOTUList_for_samples=" << samples << "\n";
 380             }else if ((abund != -1) && (samples == -1)) { //we want all OTUs with this relabund
 381                 out2 << "Samples\tOTUList_for_abund=" << abund << "\n";
 382             }else if ((abund != -1) && (samples != -1)) { //we want only OTUs with this relabund for this number of samples
 383                 out2 << "Relabund%\tOTUList_for_samples=" << samples << "\n";
 384             }
 385
 386             for (map<int, vector<string> >::iterator it = otuNames.begin(); it != otuNames.end(); it++) {
 387                 if (m->control_pressed) { break; }
 388
 389                 vector<string> temp = it->second;
 390                 string list = m->makeList(temp);
 391
 392                 out2 << it->first << '\t' << list << endl;
 393             }
 394
 395             out2.close();
 396         }
 397
 398         return 0;
 399     }
 400         catch(exception& e) {
 401                 m->errorOut(e, "GetCoreMicroBiomeCommand", "createTable");
 402                 exit(1);
 403         }
 404 }
 405
 406 //**********************************************************************************************************************
 407
 408