getrelabundcommand.cpp

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