classifysharedcommand.cpp

   1 //
   2 //  classifysharedcommand.cpp
   3 //  Mothur
   4 //
   5 //  Created by Abu Zaher Md. Faridee on 8/13/12.
   6 //  Copyright (c) 2012 Schloss Lab. All rights reserved.
   7 //
   8
   9 #include "classifysharedcommand.h"
  10 #include "randomforest.hpp"
  11 #include "decisiontree.hpp"
  12 #include "rftreenode.hpp"
  13
  14 //**********************************************************************************************************************
  15 vector<string> ClassifySharedCommand::setParameters(){
  16         try {
  17                 //CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
  18         CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none","summary",false,true,true); parameters.push_back(pshared);
  19         CommandParameter pdesign("design", "InputTypes", "", "", "none", "none", "none","",false,true,true); parameters.push_back(pdesign);
  20         CommandParameter potupersplit("otupersplit", "Multiple", "log2-squareroot", "log2", "", "", "","",false,false); parameters.push_back(potupersplit);
  21         CommandParameter psplitcriteria("splitcriteria", "Multiple", "gainratio-infogain", "gainratio", "", "", "","",false,false); parameters.push_back(psplitcriteria);
  22                 CommandParameter pnumtrees("numtrees", "Number", "", "100", "", "", "","",false,false); parameters.push_back(pnumtrees);
  23
  24         CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
  25                 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
  26                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
  27                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
  28
  29                 vector<string> myArray;
  30                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  31                 return myArray;
  32         }
  33         catch(exception& e) {
  34                 m->errorOut(e, "ClassifySharedCommand", "setParameters");
  35                 exit(1);
  36         }
  37 }
  38 //**********************************************************************************************************************
  39 string ClassifySharedCommand::getHelpString(){
  40         try {
  41                 string helpString = "";
  42                 helpString += "The classify.shared command allows you to ....\n";
  43                 helpString += "The classify.shared command parameters are: shared, design, label, groups, otupersplit.\n";
  44         helpString += "The label parameter is used to analyze specific labels in your input.\n";
  45                 helpString += "The groups parameter allows you to specify which of the groups in your designfile you would like analyzed.\n";
  46                 helpString += "The classify.shared should be in the following format: \n";
  47                 helpString += "classify.shared(shared=yourSharedFile, design=yourDesignFile)\n";
  48                 return helpString;
  49         }
  50         catch(exception& e) {
  51                 m->errorOut(e, "ClassifySharedCommand", "getHelpString");
  52                 exit(1);
  53         }
  54 }
  55 //**********************************************************************************************************************
  56 string ClassifySharedCommand::getOutputPattern(string type) {
  57     try {
  58         string pattern = "";
  59
  60         if (type == "summary") {  pattern = "[filename],[distance],summary"; } //makes file like: amazon.0.03.fasta
  61         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
  62
  63         return pattern;
  64     }
  65     catch(exception& e) {
  66         m->errorOut(e, "ClassifySharedCommand", "getOutputPattern");
  67         exit(1);
  68     }
  69 }
  70 //**********************************************************************************************************************
  71
  72 ClassifySharedCommand::ClassifySharedCommand() {
  73   try {
  74     abort = true; calledHelp = true;
  75     setParameters();
  76     vector<string> tempOutNames;
  77     outputTypes["summary"] = tempOutNames;
  78   }
  79   catch(exception& e) {
  80     m->errorOut(e, "ClassifySharedCommand", "ClassifySharedCommand");
  81     exit(1);
  82   }
  83 }
  84 //**********************************************************************************************************************
  85 ClassifySharedCommand::ClassifySharedCommand(string option) {
  86   try {
  87     abort = false; calledHelp = false;
  88     allLines = 1;
  89
  90       //allow user to run help
  91     if(option == "help") { help(); abort = true; calledHelp = true; }
  92     else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  93
  94     else {
  95         //valid paramters for this command
  96       vector<string> myArray = setParameters();
  97
  98       OptionParser parser(option);
  99       map<string,string> parameters = parser.getParameters();
 100
 101       ValidParameters validParameter;
 102       map<string,string>::iterator it;
 103         //check to make sure all parameters are valid for command
 104       for (it = parameters.begin(); it != parameters.end(); it++) {
 105         if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 106       }
 107
 108         vector<string> tempOutNames;
 109         outputTypes["summary"] = tempOutNames;
 110
 111         //if the user changes the input directory command factory will send this info to us in the output parameter
 112       string inputDir = validParameter.validFile(parameters, "inputdir", false);
 113       if (inputDir == "not found"){     inputDir = "";          }
 114       else {
 115         string path;
 116         it = parameters.find("shared");
 117           //user has given a shared file
 118         if(it != parameters.end()){
 119           path = m->hasPath(it->second);
 120             //if the user has not given a path then, add inputdir. else leave path alone.
 121           if (path == "") {     parameters["shared"] = inputDir + it->second;           }
 122         }
 123
 124         it = parameters.find("design");
 125           //user has given a design file
 126         if(it != parameters.end()){
 127           path = m->hasPath(it->second);
 128             //if the user has not given a path then, add inputdir. else leave path alone.
 129           if (path == "") {     parameters["design"] = inputDir + it->second;           }
 130         }
 131
 132       }
 133
 134         //check for parameters
 135         //get shared file, it is required
 136       sharedfile = validParameter.validFile(parameters, "shared", true);
 137       if (sharedfile == "not open") { sharedfile = ""; abort = true; }
 138       else if (sharedfile == "not found") {
 139           //if there is a current shared file, use it
 140         sharedfile = m->getSharedFile();
 141         if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 142         else {  m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
 143       }else { m->setSharedFile(sharedfile); }
 144
 145         //get design file, it is required
 146       designfile = validParameter.validFile(parameters, "design", true);
 147       if (designfile == "not open") { sharedfile = ""; abort = true; }
 148       else if (designfile == "not found") {
 149           //if there is a current shared file, use it
 150         designfile = m->getDesignFile();
 151         if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
 152         else {  m->mothurOut("You have no current designfile and the design parameter is required."); m->mothurOutEndLine(); abort = true; }
 153       }else { m->setDesignFile(designfile); }
 154
 155
 156         //if the user changes the output directory command factory will send this info to us in the output parameter
 157       outputDir = validParameter.validFile(parameters, "outputdir", false);             if (outputDir == "not found"){
 158         outputDir = m->hasPath(sharedfile); //if user entered a file with a path then preserve it
 159       }
 160
 161
 162         // NEW CODE for OTU per split selection criteria
 163       otupersplit = validParameter.validFile(parameters, "otupersplit", false);
 164       if (otupersplit == "not found") { otupersplit = "log2"; }
 165       if ((otupersplit == "squareroot") || (otupersplit == "log2")) {
 166         optimumFeatureSubsetSelectionCriteria = otupersplit;
 167       }else { m->mothurOut("Not a valid OTU per split selection method. Valid OTU per split selection methods are 'log2' and 'squareroot'."); m->mothurOutEndLine(); abort = true; }
 168
 169         // splitcriteria
 170       splitcriteria = validParameter.validFile(parameters, "splitcriteria", false);
 171       if (splitcriteria == "not found") { splitcriteria = "gainratio"; }
 172       if ((splitcriteria == "gainratio") || (splitcriteria == "infogain")) {
 173         treeSplitCriterion = splitcriteria;
 174       }else { m->mothurOut("Not a valid tree splitting criterio. Valid tree splitting criteria are 'gainratio' and 'infogain'."); m->mothurOutEndLine(); abort = true; }
 175
 176
 177       string temp = validParameter.validFile(parameters, "numtrees", false); if (temp == "not found"){  temp = "100";   }
 178       m->mothurConvert(temp, numDecisionTrees);
 179
 180         //Groups must be checked later to make sure they are valid. SharedUtilities has functions of check the validity, just make to so m->setGroups() after the checks.  If you are using these with a shared file no need to check the SharedRAbundVector class will call SharedUtilites for you, kinda nice, huh?
 181       string groups = validParameter.validFile(parameters, "groups", false);
 182       if (groups == "not found") { groups = ""; }
 183       else { m->splitAtDash(groups, Groups); }
 184       m->setGroups(Groups);
 185
 186         //Commonly used to process list, rabund, sabund, shared and relabund files.  Look at "smart distancing" examples below in the execute function.
 187       string label = validParameter.validFile(parameters, "label", false);
 188       if (label == "not found") { label = ""; }
 189       else {
 190         if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 191         else { allLines = 1;  }
 192       }
 193     }
 194
 195   }
 196   catch(exception& e) {
 197     m->errorOut(e, "ClassifySharedCommand", "ClassifySharedCommand");
 198     exit(1);
 199   }
 200 }
 201 //**********************************************************************************************************************
 202 int ClassifySharedCommand::execute() {
 203   try {
 204
 205     if (abort == true) { if (calledHelp) { return 0; }  return 2;       }
 206
 207     InputData input(sharedfile, "sharedfile");
 208     vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
 209
 210     //read design file
 211     designMap.readDesignMap(designfile);
 212
 213     string lastLabel = lookup[0]->getLabel();
 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];  }  return 0; }
 221
 222       if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 223
 224         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 225
 226         processSharedAndDesignData(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.getSharedRAbundVectors(lastLabel);
 237         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 238
 239         processSharedAndDesignData(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) { return 0; }
 253
 254         //get next line to process
 255       lookup = input.getSharedRAbundVectors();
 256     }
 257
 258     if (m->control_pressed) {  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.getSharedRAbundVectors(lastLabel);
 277
 278       m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 279
 280       processSharedAndDesignData(lookup);
 281
 282       for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 283
 284     }
 285
 286       m->mothurOutEndLine();
 287       m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 288       for (int i = 0; i < outputNames.size(); i++) {    m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 289       m->mothurOutEndLine();
 290
 291     return 0;
 292
 293   }
 294   catch(exception& e) {
 295     m->errorOut(e, "ClassifySharedCommand", "execute");
 296     exit(1);
 297   }
 298 }
 299 //**********************************************************************************************************************
 300
 301 void ClassifySharedCommand::processSharedAndDesignData(vector<SharedRAbundVector*> lookup){
 302     try {
 303 //    for (int i = 0; i < designMap->getNamesOfGroups().size(); i++) {
 304 //      string groupName = designMap->getNamesOfGroups()[i];
 305 //      cout << groupName << endl;
 306 //    }
 307
 308 //    for (int i = 0; i < designMap->getNumSeqs(); i++) {
 309 //      string sharedGroupName = designMap->getNamesSeqs()[i];
 310 //      string treatmentName = designMap->getGroup(sharedGroupName);
 311 //      cout << sharedGroupName << " : " << treatmentName <<  endl;
 312 //    }
 313
 314         map<string, int> treatmentToIntMap;
 315         map<int, string> intToTreatmentMap;
 316         for (int  i = 0; i < designMap.getNumGroups(); i++) {
 317             string treatmentName = designMap.getNamesOfGroups()[i];
 318             treatmentToIntMap[treatmentName] = i;
 319             intToTreatmentMap[i] = treatmentName;
 320         }
 321
 322         int numSamples = lookup.size();
 323         int numFeatures = lookup[0]->getNumBins();
 324
 325         int numRows = numSamples;
 326         int numColumns = numFeatures + 1;           // extra one space needed for the treatment/outcome
 327
 328         vector< vector<int> > dataSet(numRows, vector<int>(numColumns, 0));
 329
 330         for (int i = 0; i < lookup.size(); i++) {
 331             string sharedGroupName = lookup[i]->getGroup();
 332             string treatmentName = designMap.getGroup(sharedGroupName);
 333
 334             int j = 0;
 335             for (; j < lookup[i]->getNumBins(); j++) {
 336                 int otuCount = lookup[i]->getAbundance(j);
 337                 dataSet[i][j] = otuCount;
 338             }
 339             dataSet[i][j] = treatmentToIntMap[treatmentName];
 340         }
 341
 342         RandomForest randomForest(dataSet, numDecisionTrees, treeSplitCriterion);
 343         randomForest.populateDecisionTrees();
 344         randomForest.calcForrestErrorRate();
 345
 346         map<string, string> variables;
 347         variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
 348         variables["[distance]"] = lookup[0]->getLabel();
 349         string filename = getOutputFileName("summary", variables);
 350         outputNames.push_back(filename); outputTypes["summary"].push_back(filename);
 351
 352         randomForest.calcForrestVariableImportance(filename);
 353
 354         m->mothurOutEndLine();
 355     }
 356     catch(exception& e) {
 357         m->errorOut(e, "ClassifySharedCommand", "processSharedAndDesignData");
 358         exit(1);
 359     }
 360 }
 361 //**********************************************************************************************************************
 362