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