metastatscommand.cpp

   1 /*
   2  *  metastatscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/16/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "metastatscommand.h"
  11 #include "metastats.h"
  12
  13 //**********************************************************************************************************************
  14
  15 MetaStatsCommand::MetaStatsCommand(string option) {
  16         try {
  17                 globaldata = GlobalData::getInstance();
  18                 abort = false;
  19                 allLines = 1;
  20                 labels.clear();
  21
  22                 //allow user to run help
  23                 if(option == "help") { help(); abort = true; }
  24
  25                 else {
  26                         //valid paramters for this command
  27                         string AlignArray[] =  {"groups","label","outputdir","iters","threshold","g","design","inputdir"};
  28                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  29
  30                         OptionParser parser(option);
  31                         map<string,string> parameters = parser.getParameters();
  32
  33                         ValidParameters validParameter;
  34
  35                         //check to make sure all parameters are valid for command
  36                         map<string,string>::iterator it;
  37                         for (it = parameters.begin(); it != parameters.end(); it++) {
  38                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  39                         }
  40
  41                         //if the user changes the output directory command factory will send this info to us in the output parameter
  42                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
  43                                 outputDir = "";
  44                                 outputDir += m->hasPath(globaldata->inputFileName); //if user entered a file with a path then preserve it
  45                         }
  46
  47                         //if the user changes the input directory command factory will send this info to us in the output parameter
  48                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  49                         if (inputDir == "not found"){   inputDir = "";          }
  50                         else {
  51                                 string path;
  52                                 it = parameters.find("design");
  53                                 //user has given a template file
  54                                 if(it != parameters.end()){
  55                                         path = m->hasPath(it->second);
  56                                         //if the user has not given a path then, add inputdir. else leave path alone.
  57                                         if (path == "") {       parameters["design"] = inputDir + it->second;           }
  58                                 }
  59                         }
  60
  61                         //check for required parameters
  62                         designfile = validParameter.validFile(parameters, "design", true);
  63                         if (designfile == "not open") { abort = true; }
  64                         else if (designfile == "not found") {  designfile = "";  m->mothurOut("You must provide an design file."); m->mothurOutEndLine(); abort = true; }
  65
  66                         //make sure the user has already run the read.otu command
  67                         if ((globaldata->getSharedFile() == "")) {
  68                                  m->mothurOut("You must read a list and a group, or a shared file before you can use the normalize.shared command."); m->mothurOutEndLine(); abort = true;
  69                         }
  70
  71                         //check for optional parameter and set defaults
  72                         // ...at some point should added some additional type checking...
  73                         label = validParameter.validFile(parameters, "label", false);
  74                         if (label == "not found") { label = ""; }
  75                         else {
  76                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
  77                                 else { allLines = 1;  }
  78                         }
  79
  80                         //if the user has not specified any labels use the ones from read.otu
  81                         if (label == "") {
  82                                 allLines = globaldata->allLines;
  83                                 labels = globaldata->labels;
  84                         }
  85
  86                         groups = validParameter.validFile(parameters, "groups", false);
  87                         if (groups == "not found") { groups = ""; pickedGroups = false; }
  88                         else {
  89                                 pickedGroups = true;
  90                                 m->splitAtDash(groups, Groups);
  91                                 globaldata->Groups = Groups;
  92                         }
  93
  94                         string temp = validParameter.validFile(parameters, "iters", false);                     if (temp == "not found") { temp = "1000"; }
  95                         convert(temp, iters);
  96
  97                         temp = validParameter.validFile(parameters, "threshold", false);                        if (temp == "not found") { temp = "0.05"; }
  98                         convert(temp, threshold);
  99
 100                         //temp = validParameter.validFile(parameters, "g", false);                                      if (temp == "not found") { temp = "0"; }
 101                         //convert(temp, g);
 102                 }
 103
 104         }
 105         catch(exception& e) {
 106                 m->errorOut(e, "MetaStatsCommand", "MetaStatsCommand");
 107                 exit(1);
 108         }
 109 }
 110
 111 //**********************************************************************************************************************
 112
 113 void MetaStatsCommand::help(){
 114         try {
 115                 m->mothurOut("Paper reference goes here....\n");
 116                 m->mothurOut("The metastats command can only be executed after a successful read.otu command of a list and group or shared file.\n");
 117                 m->mothurOut("The metastats command outputs a .... file. \n");
 118                 m->mothurOut("The metastats command parameters are iters, threshold, groups, design and label.  No parameters are required.\n");
 119                 m->mothurOut("The design parameter allows you to assign your groups to sets when you are running metastat. mothur will run all pairwise comparisons of the sets. It is required. \n");
 120                 m->mothurOut("The design file looks like the group file.  It is a 2 column tab delimited file, where the first column is the group name and the second column is the set the group belongs to.\n");
 121                 m->mothurOut("The iters parameter ...\n");
 122                 m->mothurOut("The threshold parameter ...\n");
 123                 m->mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like included. The group names are separated by dashes.\n");
 124                 m->mothurOut("The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n");
 125                 m->mothurOut("The metastats command should be in the following format: metastats(groups=yourGroups, label=yourLabels).\n");
 126                 m->mothurOut("Example metastats(groups=A-B-C).\n");
 127                 m->mothurOut("The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n");
 128                 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
 129
 130         }
 131         catch(exception& e) {
 132                 m->errorOut(e, "MetaStatsCommand", "help");
 133                 exit(1);
 134         }
 135 }
 136
 137 //**********************************************************************************************************************
 138
 139 MetaStatsCommand::~MetaStatsCommand(){}
 140
 141 //**********************************************************************************************************************
 142
 143 int MetaStatsCommand::execute(){
 144         try {
 145
 146                 if (abort == true) { return 0; }
 147
 148                 designMap = new GroupMap(designfile);
 149                 designMap->readDesignMap();
 150
 151                 read = new ReadOTUFile(globaldata->inputFileName);
 152                 read->read(&*globaldata);
 153                 input = globaldata->ginput;
 154                 lookup = input->getSharedRAbundVectors();
 155                 string lastLabel = lookup[0]->getLabel();
 156
 157                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 158                 set<string> processedLabels;
 159                 set<string> userLabels = labels;
 160
 161                 //as long as you are not at the end of the file or done wih the lines you want
 162                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 163
 164                         if (m->control_pressed) {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } globaldata->Groups.clear(); delete read;  for (int i = 0; i < outputNames.size(); i++) {     remove(outputNames[i].c_str()); } return 0; }
 165
 166                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 167
 168                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 169                                 process(lookup);
 170
 171                                 processedLabels.insert(lookup[0]->getLabel());
 172                                 userLabels.erase(lookup[0]->getLabel());
 173                         }
 174
 175                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 176                                 string saveLabel = lookup[0]->getLabel();
 177
 178                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 179                                 lookup = input->getSharedRAbundVectors(lastLabel);
 180                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 181
 182                                 process(lookup);
 183
 184                                 processedLabels.insert(lookup[0]->getLabel());
 185                                 userLabels.erase(lookup[0]->getLabel());
 186
 187                                 //restore real lastlabel to save below
 188                                 lookup[0]->setLabel(saveLabel);
 189                         }
 190
 191                         lastLabel = lookup[0]->getLabel();
 192                         //prevent memory leak
 193                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 194
 195                         if (m->control_pressed) {  globaldata->Groups.clear(); delete read;   for (int i = 0; i < outputNames.size(); i++) {    remove(outputNames[i].c_str()); } return 0; }
 196
 197                         //get next line to process
 198                         lookup = input->getSharedRAbundVectors();
 199                 }
 200
 201                 if (m->control_pressed) { globaldata->Groups.clear(); delete read;  for (int i = 0; i < outputNames.size(); i++) {      remove(outputNames[i].c_str()); }  return 0; }
 202
 203                 //output error messages about any remaining user labels
 204                 set<string>::iterator it;
 205                 bool needToRun = false;
 206                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 207                         m->mothurOut("Your file does not include the label " + *it);
 208                         if (processedLabels.count(lastLabel) != 1) {
 209                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 210                                 needToRun = true;
 211                         }else {
 212                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 213                         }
 214                 }
 215
 216                 //run last label if you need to
 217                 if (needToRun == true)  {
 218                         for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
 219                         lookup = input->getSharedRAbundVectors(lastLabel);
 220
 221                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 222
 223                         process(lookup);
 224
 225                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 226                 }
 227
 228                 //reset groups parameter
 229                 globaldata->Groups.clear();
 230                 delete input; globaldata->ginput = NULL;
 231                 delete read;
 232
 233                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str()); } return 0;}
 234
 235                 m->mothurOutEndLine();
 236                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 237                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 238                 m->mothurOutEndLine();
 239
 240                 return 0;
 241         }
 242         catch(exception& e) {
 243                 m->errorOut(e, "MetaStatsCommand", "execute");
 244                 exit(1);
 245         }
 246 }
 247 //**********************************************************************************************************************
 248
 249 int MetaStatsCommand::process(vector<SharedRAbundVector*>& thisLookUp){
 250         try {
 251                 if (pickedGroups) { eliminateZeroOTUS(thisLookUp); }
 252
 253                 string outputFileName = outputDir +  m->getRootName(m->getSimpleName(globaldata->inputFileName)) + thisLookUp[0]->getLabel() + ".metastats";
 254                 outputNames.push_back(outputFileName);
 255                 int nameLength = outputFileName.length();
 256                 char * output = new char[nameLength];
 257                 strcpy(output, outputFileName.c_str());
 258
 259                 //build matrix from shared rabunds
 260                 double** data;
 261                 data = new double*[thisLookUp[0]->getNumBins()];
 262
 263                 for (int j = 0; j < thisLookUp[0]->getNumBins(); j++) {
 264                         data[j] = new double[thisLookUp.size()];
 265                         for (int i = 0; i < thisLookUp.size(); i++) {
 266                                 data[j][i] = (thisLookUp[i]->getAbundance(j));
 267                         }
 268                 }
 269
 270                 if (g == 0) { g = thisLookUp.size() / 2; }
 271
 272                 metastat_main(output, thisLookUp[0]->getNumBins(), thisLookUp.size(), threshold, iters, data, g);
 273
 274                 return 0;
 275         }
 276         catch(exception& e) {
 277                 m->errorOut(e, "MetaStatsCommand", "normalize");
 278                 exit(1);
 279         }
 280 }
 281 //**********************************************************************************************************************
 282 int MetaStatsCommand::eliminateZeroOTUS(vector<SharedRAbundVector*>& thislookup) {
 283         try {
 284
 285                 vector<SharedRAbundVector*> newLookup;
 286                 for (int i = 0; i < thislookup.size(); i++) {
 287                         SharedRAbundVector* temp = new SharedRAbundVector();
 288                         temp->setLabel(thislookup[i]->getLabel());
 289                         temp->setGroup(thislookup[i]->getGroup());
 290                         newLookup.push_back(temp);
 291                 }
 292
 293                 //for each bin
 294                 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
 295                         if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
 296
 297                         //look at each sharedRabund and make sure they are not all zero
 298                         bool allZero = true;
 299                         for (int j = 0; j < thislookup.size(); j++) {
 300                                 if (thislookup[j]->getAbundance(i) != 0) { allZero = false;  break;  }
 301                         }
 302
 303                         //if they are not all zero add this bin
 304                         if (!allZero) {
 305                                 for (int j = 0; j < thislookup.size(); j++) {
 306                                         newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
 307                                 }
 308                         }
 309                 }
 310
 311                 for (int j = 0; j < thislookup.size(); j++) {  delete thislookup[j];  }
 312
 313                 thislookup = newLookup;
 314
 315                 return 0;
 316
 317         }
 318         catch(exception& e) {
 319                 m->errorOut(e, "MetaStatsCommand", "eliminateZeroOTUS");
 320                 exit(1);
 321         }
 322 }
 323 //**********************************************************************************************************************