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