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["coremicrobiome"] = 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             vector<string> tempOutNames;
 115             outputTypes["coremicrobiome"] = tempOutNames;
 116
 117                         //check for parameters
 118             sharedfile = validParameter.validFile(parameters, "shared", true);
 119                         if (sharedfile == "not open") { abort = true; }
 120                         else if (sharedfile == "not found") { sharedfile = ""; }
 121                         else { inputFileName = sharedfile; format = "sharedfile"; m->setSharedFile(sharedfile); }
 122
 123                         relabundfile = validParameter.validFile(parameters, "relabund", true);
 124                         if (relabundfile == "not open") { abort = true; }
 125                         else if (relabundfile == "not found") { relabundfile = ""; }
 126                         else { inputFileName = relabundfile; format = "relabund"; m->setRelAbundFile(relabundfile); }
 127
 128             if ((relabundfile == "") && (sharedfile == "")) {
 129                                 //is there are current file available for either of these?
 130                                 //give priority to shared, then relabund
 131                                 sharedfile = m->getSharedFile();
 132                                 if (sharedfile != "") {  inputFileName = sharedfile; format="sharedfile"; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 133                                 else {
 134                                         relabundfile = m->getRelAbundFile();
 135                                         if (relabundfile != "") {  inputFileName = relabundfile; format="relabund"; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
 136                                         else {
 137                                                 m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine();
 138                                                 abort = true;
 139                                         }
 140                                 }
 141                         }
 142
 143             //if the user changes the output directory command factory will send this info to us in the output parameter
 144                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 145                                 outputDir = m->hasPath(inputFileName); //if user entered a file with a path then preserve it
 146                         }
 147
 148             string groups = validParameter.validFile(parameters, "groups", false);
 149                         if (groups == "not found") { groups = ""; }
 150                         else { m->splitAtDash(groups, Groups); }
 151                         m->setGroups(Groups);
 152
 153             string label = validParameter.validFile(parameters, "label", false);
 154                         if (label == "not found") { label = ""; }
 155                         else {
 156                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 157                                 else { allLines = 1;  }
 158                         }
 159
 160             output = validParameter.validFile(parameters, "output", false);             if(output == "not found"){      output = "fraction"; }
 161
 162                         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"; }
 163
 164             string temp = validParameter.validFile(parameters, "abundance", false);     if (temp == "not found"){       temp = "-1";    }
 165                         m->mothurConvert(temp, abund);
 166
 167             if (abund != -1) {
 168                 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"); }
 169             }
 170
 171             temp = validParameter.validFile(parameters, "samples", false);      if (temp == "not found"){       temp = "-1";    }
 172                         m->mothurConvert(temp, samples);
 173
 174                 }
 175
 176         }
 177         catch(exception& e) {
 178                 m->errorOut(e, "GetCoreMicroBiomeCommand", "GetCoreMicroBiomeCommand");
 179                 exit(1);
 180         }
 181 }
 182 //**********************************************************************************************************************
 183
 184 int GetCoreMicroBiomeCommand::execute(){
 185         try {
 186
 187                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 188
 189         InputData input(inputFileName, format);
 190         vector<SharedRAbundFloatVector*> lookup = input.getSharedRAbundFloatVectors();
 191         string lastLabel = lookup[0]->getLabel();
 192
 193         if (samples != -1) {
 194             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(); }
 195         }
 196
 197
 198         //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 199         set<string> processedLabels;
 200         set<string> userLabels = labels;
 201
 202         //as long as you are not at the end of the file or done wih the lines you want
 203         while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 204
 205             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; }
 206
 207             if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 208
 209                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 210
 211                 createTable(lookup);
 212
 213                 processedLabels.insert(lookup[0]->getLabel());
 214                 userLabels.erase(lookup[0]->getLabel());
 215             }
 216
 217             if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 218                 string saveLabel = lookup[0]->getLabel();
 219
 220                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 221                 lookup = input.getSharedRAbundFloatVectors(lastLabel);
 222                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 223
 224                 createTable(lookup);
 225
 226                 processedLabels.insert(lookup[0]->getLabel());
 227                 userLabels.erase(lookup[0]->getLabel());
 228
 229                 //restore real lastlabel to save below
 230                 lookup[0]->setLabel(saveLabel);
 231             }
 232
 233             lastLabel = lookup[0]->getLabel();
 234             //prevent memory leak
 235             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 236
 237             if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); }  return 0; }
 238
 239             //get next line to process
 240             lookup = input.getSharedRAbundFloatVectors();
 241         }
 242
 243         if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); }  return 0; }
 244
 245         //output error messages about any remaining user labels
 246         set<string>::iterator it;
 247         bool needToRun = false;
 248         for (it = userLabels.begin(); it != userLabels.end(); it++) {
 249             m->mothurOut("Your file does not include the label " + *it);
 250             if (processedLabels.count(lastLabel) != 1) {
 251                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 252                 needToRun = true;
 253             }else {
 254                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 255             }
 256         }
 257
 258         //run last label if you need to
 259         if (needToRun == true)  {
 260             for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
 261             lookup = input.getSharedRAbundFloatVectors(lastLabel);
 262
 263             m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 264
 265             createTable(lookup);
 266
 267             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 268         }
 269
 270         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); }  return 0; }
 271
 272         //output files created by command
 273                 m->mothurOutEndLine();
 274                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 275                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 276                 m->mothurOutEndLine();
 277         return 0;
 278
 279     }
 280         catch(exception& e) {
 281                 m->errorOut(e, "GetCoreMicroBiomeCommand", "execute");
 282                 exit(1);
 283         }
 284 }
 285 //**********************************************************************************************************************
 286
 287 int GetCoreMicroBiomeCommand::createTable(vector<SharedRAbundFloatVector*>& lookup){
 288         try {
 289
 290         string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + ".core.microbiome";
 291         outputNames.push_back(outputFileName);  outputTypes["coremicrobiome"].push_back(outputFileName);
 292                 ofstream out;
 293                 m->openOutputFile(outputFileName, out);
 294
 295         int numSamples = lookup.size();
 296         int numOtus = lookup[0]->getNumBins();
 297
 298         //table is 100 by numsamples
 299         //question we are answering is: what fraction of OTUs in a study have a relative abundance at or above %X
 300         //in at least %Y samples. x goes from 0 to 100, y from 1 to numSamples
 301         vector< vector<double> > table; table.resize(101);
 302         for (int i = 0; i < table.size(); i++) { table[i].resize(numSamples, 0.0); }
 303
 304         map<int, vector<string> > otuNames;
 305         if ((abund != -1) && (samples == -1)) { //fill with all samples
 306             for (int i = 0; i < numSamples; i++) {
 307                 vector<string> temp;
 308                 otuNames[i+1] = temp;
 309             }
 310         }else if ((abund == -1) && (samples != -1)) { //fill with all relabund
 311             for (int i = 0; i < 101; i++) {
 312                 vector<string> temp;
 313                 otuNames[i] = temp;
 314             }
 315         }else if ((abund != -1) && (samples != -1)) { //only one line is wanted
 316             vector<string> temp;
 317             otuNames[abund] = temp;
 318         }
 319
 320         for (int i = 0; i < numOtus; i++) {
 321
 322             if (m->control_pressed) { break; }
 323
 324             //count number of samples in this otu with a relabund >= spot in count
 325             vector<int> counts; counts.resize(101, 0);
 326
 327             for (int j = 0; j < lookup.size(); j++) {
 328                 double relabund = lookup[j]->getAbundance(i);
 329                 int wholeRelabund = (int) (floor(relabund*100));
 330                 for (int k = 0; k < wholeRelabund+1; k++) { counts[k]++; }
 331             }
 332
 333             //add this otus info to table
 334             for (int j = 0; j < table.size(); j++) {
 335                 for (int k = 0; k < counts[j]; k++) { table[j][k]++; }
 336
 337                 if ((abund == -1) && (samples != -1)) { //we want all OTUs with this number of samples
 338                     if (counts[j] >= samples) { otuNames[j].push_back(m->currentBinLabels[i]); }
 339                 }else if ((abund != -1) && (samples == -1)) { //we want all OTUs with this relabund
 340                     if (j == abund) {
 341                         for (int k = 0; k < counts[j]; k++) {  otuNames[k+1].push_back(m->currentBinLabels[i]); }
 342                     }
 343                 }else if ((abund != -1) && (samples != -1)) { //we want only OTUs with this relabund for this number of samples
 344                     if ((j == abund) && (counts[j] >= samples)) {
 345                         otuNames[j].push_back(m->currentBinLabels[i]);
 346                     }
 347                 }
 348             }
 349         }
 350
 351         //format output
 352                 if (output == "fraction") { out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); }
 353         out << "NumSamples\t";
 354
 355         //convert table counts to percents
 356         for (int i = 0; i < table.size(); i++) {
 357             out << "Relabund-" << i << "%\t";
 358             if (m->control_pressed) { break; }
 359             for (int j = 0; j < table[i].size(); j++) {  if (output == "fraction") { table[i][j] /= (double) numOtus; } }
 360         }
 361         out << endl;
 362
 363         for (int i = 0; i < numSamples; i++) {
 364             if (m->control_pressed) { break; }
 365             out << i+1 << '\t';
 366             for (int j = 0; j < table.size(); j++) {  out << table[j][i] << '\t'; }
 367             out << endl;
 368         }
 369
 370         out.close();
 371
 372         if (m->control_pressed) { return 0; }
 373
 374         if ((samples != -1) || (abund != -1))  {
 375             string outputFileName2 = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + ".core.microbiomelist";
 376             outputNames.push_back(outputFileName2);  outputTypes["coremicrobiome"].push_back(outputFileName2);
 377             ofstream out2;
 378             m->openOutputFile(outputFileName2, out2);
 379
 380             if ((abund == -1) && (samples != -1)) { //we want all OTUs with this number of samples
 381                 out2 << "Relabund%\tOTUList_for_samples=" << samples << "\n";
 382             }else if ((abund != -1) && (samples == -1)) { //we want all OTUs with this relabund
 383                 out2 << "Samples\tOTUList_for_abund=" << abund << "\n";
 384             }else if ((abund != -1) && (samples != -1)) { //we want only OTUs with this relabund for this number of samples
 385                 out2 << "Relabund%\tOTUList_for_samples=" << samples << "\n";
 386             }
 387
 388             for (map<int, vector<string> >::iterator it = otuNames.begin(); it != otuNames.end(); it++) {
 389                 if (m->control_pressed) { break; }
 390
 391                 vector<string> temp = it->second;
 392                 string list = m->makeList(temp);
 393
 394                 out2 << it->first << '\t' << list << endl;
 395             }
 396
 397             out2.close();
 398         }
 399
 400         return 0;
 401     }
 402         catch(exception& e) {
 403                 m->errorOut(e, "GetCoreMicroBiomeCommand", "createTable");
 404                 exit(1);
 405         }
 406 }
 407
 408 //**********************************************************************************************************************
 409
 410