otuassociationcommand.cpp

   1 /*
   2  *  otuassociationcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 1/19/12.
   6  *  Copyright 2012 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "otuassociationcommand.h"
  11 #include "linearalgebra.h"
  12
  13 //**********************************************************************************************************************
  14 vector<string> OTUAssociationCommand::setParameters(){
  15         try {
  16                 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pshared);
  17                 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(prelabund);
  18                 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
  19                 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
  20                 CommandParameter pmethod("method", "Multiple", "pearson-spearman-kendall", "pearson", "", "", "",false,false); parameters.push_back(pmethod);
  21                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  22                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  23
  24                 vector<string> myArray;
  25                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  26                 return myArray;
  27         }
  28         catch(exception& e) {
  29                 m->errorOut(e, "OTUAssociationCommand", "setParameters");
  30                 exit(1);
  31         }
  32 }
  33 //**********************************************************************************************************************
  34 string OTUAssociationCommand::getHelpString(){
  35         try {
  36                 string helpString = "";
  37                 helpString += "The otu.association command reads a shared or relabund file and calculates the correlation coefficients between otus.\n";
  38                 helpString += "The otu.association command parameters are shared, relabund, groups, method and label.  The shared or relabund parameter is required.\n";
  39                 helpString += "The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes.\n";
  40                 helpString += "The label parameter allows you to select what distances level you would like used, and are also separated by dashes.\n";
  41                 helpString += "The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n";
  42                 helpString += "The otu.association command should be in the following format: otu.association(shared=yourSharedFile, method=yourMethod).\n";
  43                 helpString += "Example otu.association(shared=genus.pool.shared, method=kendall).\n";
  44                 helpString += "The otu.association command outputs a .otu.corr file.\n";
  45                 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
  46                 return helpString;
  47         }
  48         catch(exception& e) {
  49                 m->errorOut(e, "OTUAssociationCommand", "getHelpString");
  50                 exit(1);
  51         }
  52 }
  53 //**********************************************************************************************************************
  54 OTUAssociationCommand::OTUAssociationCommand(){
  55         try {
  56                 abort = true; calledHelp = true;
  57                 setParameters();
  58                 vector<string> tempOutNames;
  59                 outputTypes["otu.corr"] = tempOutNames;
  60         }
  61         catch(exception& e) {
  62                 m->errorOut(e, "OTUAssociationCommand", "OTUAssociationCommand");
  63                 exit(1);
  64         }
  65 }
  66 //**********************************************************************************************************************
  67 OTUAssociationCommand::OTUAssociationCommand(string option)  {
  68         try {
  69                 abort = false; calledHelp = false;
  70                 allLines = 1;
  71
  72                 //allow user to run help
  73                 if(option == "help") { help(); abort = true; calledHelp = true; }
  74                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  75
  76                 else {
  77                         vector<string> myArray = setParameters();
  78
  79                         OptionParser parser(option);
  80                         map<string, string> parameters = parser.getParameters();
  81
  82                         ValidParameters validParameter;
  83                         map<string, string>::iterator it;
  84
  85                         //check to make sure all parameters are valid for command
  86                         for (it = parameters.begin(); it != parameters.end(); it++) {
  87                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  88                         }
  89
  90                         vector<string> tempOutNames;
  91                         outputTypes["otu.corr"] = tempOutNames;
  92
  93                         //if the user changes the input directory command factory will send this info to us in the output parameter
  94                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  95                         if (inputDir == "not found"){   inputDir = "";          }
  96                         else {
  97                                 string path;
  98                                 it = parameters.find("shared");
  99                                 //user has given a template file
 100                                 if(it != parameters.end()){
 101                                         path = m->hasPath(it->second);
 102                                         //if the user has not given a path then, add inputdir. else leave path alone.
 103                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 104                                 }
 105
 106                                 it = parameters.find("relabund");
 107                                 //user has given a template file
 108                                 if(it != parameters.end()){
 109                                         path = m->hasPath(it->second);
 110                                         //if the user has not given a path then, add inputdir. else leave path alone.
 111                                         if (path == "") {       parameters["relabund"] = inputDir + it->second;         }
 112                                 }
 113                         }
 114
 115
 116                         //check for required parameters
 117                         sharedfile = validParameter.validFile(parameters, "shared", true);
 118                         if (sharedfile == "not open") { abort = true; }
 119                         else if (sharedfile == "not found") { sharedfile = ""; }
 120                         else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
 121
 122                         relabundfile = validParameter.validFile(parameters, "relabund", true);
 123                         if (relabundfile == "not open") { abort = true; }
 124                         else if (relabundfile == "not found") { relabundfile = ""; }
 125                         else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
 126
 127                         groups = validParameter.validFile(parameters, "groups", false);
 128                         if (groups == "not found") { groups = "";  pickedGroups = false;  }
 129                         else {
 130                                 pickedGroups = true;
 131                                 m->splitAtDash(groups, Groups);
 132                         }
 133                         m->setGroups(Groups);
 134
 135                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(inputFileName);  }
 136
 137                         label = validParameter.validFile(parameters, "label", false);
 138                         if (label == "not found") { label = ""; }
 139                         else {
 140                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 141                                 else { allLines = 1;  }
 142                         }
 143
 144                         if ((relabundfile == "") && (sharedfile == "")) {
 145                                 //is there are current file available for any of these?
 146                                 //give priority to shared, then relabund
 147                                 //if there is a current shared file, use it
 148                                 sharedfile = m->getSharedFile();
 149                                 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 150                                 else {
 151                                         relabundfile = m->getRelAbundFile();
 152                                         if (relabundfile != "") { inputFileName = relabundfile;  m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
 153                                         else {
 154                                                 m->mothurOut("You must provide either a shared or relabund file."); m->mothurOutEndLine(); abort = true;
 155                                         }
 156                                 }
 157                         }
 158
 159
 160                         if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared or relabund file."); m->mothurOutEndLine(); abort = true;  }
 161
 162                         method = validParameter.validFile(parameters, "method", false);         if (method == "not found"){     method = "pearson";             }
 163
 164                         if ((method != "pearson") && (method != "spearman") && (method != "kendall")) { m->mothurOut(method + " is not a valid method. Valid methods are pearson, spearman, and kendall."); m->mothurOutEndLine(); abort = true; }
 165
 166                 }
 167         }
 168         catch(exception& e) {
 169                 m->errorOut(e, "OTUAssociationCommand", "OTUAssociationCommand");
 170                 exit(1);
 171         }
 172 }
 173 //**********************************************************************************************************************
 174
 175 int OTUAssociationCommand::execute(){
 176         try {
 177
 178                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 179
 180                 //function are identical just different datatypes
 181                 if (sharedfile != "")                   {  processShared();             }
 182                 else if (relabundfile != "")    {  processRelabund();   }
 183
 184                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       m->mothurRemove(outputNames[i]); } return 0; }
 185
 186                 m->mothurOutEndLine();
 187                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 188                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 189                 m->mothurOutEndLine();
 190
 191                 return 0;
 192         }
 193         catch(exception& e) {
 194                 m->errorOut(e, "OTUAssociationCommand", "execute");
 195                 exit(1);
 196         }
 197 }
 198 //**********************************************************************************************************************
 199 int OTUAssociationCommand::processShared(){
 200         try {
 201                 InputData* input = new InputData(sharedfile, "sharedfile");
 202                 vector<SharedRAbundVector*> lookup = input->getSharedRAbundVectors();
 203                 string lastLabel = lookup[0]->getLabel();
 204
 205                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 206                 set<string> processedLabels;
 207                 set<string> userLabels = labels;
 208
 209                 //as long as you are not at the end of the file or done wih the lines you want
 210                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 211
 212                         if (m->control_pressed) {  delete input; return 0;  }
 213
 214                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 215                                 processedLabels.insert(lookup[0]->getLabel());
 216                                 userLabels.erase(lookup[0]->getLabel());
 217
 218                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 219                                 process(lookup);
 220                         }
 221
 222                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 223                                 string saveLabel = lookup[0]->getLabel();
 224
 225                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 226                                 lookup = input->getSharedRAbundVectors(lastLabel);
 227
 228                                 processedLabels.insert(lookup[0]->getLabel());
 229                                 userLabels.erase(lookup[0]->getLabel());
 230
 231                                 //restore real lastlabel to save below
 232                                 lookup[0]->setLabel(saveLabel);
 233
 234                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 235                                 process(lookup);
 236                         }
 237
 238                         lastLabel = lookup[0]->getLabel();
 239
 240                         //get next line to process
 241                         //prevent memory leak
 242                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 243                         lookup = input->getSharedRAbundVectors();
 244                 }
 245
 246
 247                 if (m->control_pressed) { delete input; return 0;  }
 248
 249                 //output error messages about any remaining user labels
 250                 set<string>::iterator it;
 251                 bool needToRun = false;
 252                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 253                         m->mothurOut("Your file does not include the label " + *it);
 254                         if (processedLabels.count(lastLabel) != 1) {
 255                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 256                                 needToRun = true;
 257                         }else {
 258                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 259                         }
 260                 }
 261
 262                 //run last label if you need to
 263                 if (needToRun == true)  {
 264                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } }
 265                         lookup = input->getSharedRAbundVectors(lastLabel);
 266
 267                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 268                         process(lookup);
 269                 }
 270
 271                 delete input;
 272
 273                 return 0;
 274         }
 275         catch(exception& e) {
 276                 m->errorOut(e, "OTUAssociationCommand", "processShared");
 277                 exit(1);
 278         }
 279 }
 280 //**********************************************************************************************************************
 281 int OTUAssociationCommand::process(vector<SharedRAbundVector*>& lookup){
 282         try {
 283
 284                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + "." + method + ".otu.corr";
 285                 outputNames.push_back(outputFileName); outputTypes["shared"].push_back(outputFileName);
 286
 287                 ofstream out;
 288                 m->openOutputFile(outputFileName, out);
 289                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 290
 291                 //column headings
 292                 out << "OTUA\tOTUB\t" << method << "Coef\tSignificance\n";
 293
 294                 vector< vector<double> > xy; xy.resize(lookup[0]->getNumBins());
 295                 for (int i = 0; i < lookup[0]->getNumBins(); i++) { for (int j = 0; j < lookup.size(); j++) { xy[i].push_back(lookup[j]->getAbundance(i)); } }
 296
 297                 LinearAlgebra linear;
 298                 for (int i = 0; i < xy.size(); i++) {
 299
 300                         for (int k = 0; k < i; k++) {
 301
 302                                 if (m->control_pressed) { out.close(); return 0; }
 303
 304                                 double coef = 0.0;
 305                                 double sig = 0.0;
 306                                 if (method == "spearman")               {   coef = linear.calcSpearman(xy[i], xy[k], sig);      }
 307                                 else if (method == "pearson")   {       coef = linear.calcPearson(xy[i], xy[k], sig);   }
 308                                 else if (method == "kendall")   {       coef = linear.calcKendall(xy[i], xy[k], sig);   }
 309                                 else { m->mothurOut("[ERROR]: invalid method, choices are spearman, pearson or kendall."); m->mothurOutEndLine(); m->control_pressed = true; }
 310
 311                 out << m->binLabelsInFile[i] << '\t' << m->binLabelsInFile[k] << '\t' << coef << '\t' << sig << endl;
 312                         }
 313                 }
 314
 315                 out.close();
 316
 317
 318                 return 0;
 319
 320         }
 321         catch(exception& e) {
 322                 m->errorOut(e, "OTUAssociationCommand", "process");
 323                 exit(1);
 324         }
 325 }
 326 //**********************************************************************************************************************
 327 int OTUAssociationCommand::processRelabund(){
 328         try {
 329                 InputData* input = new InputData(relabundfile, "relabund");
 330                 vector<SharedRAbundFloatVector*> lookup = input->getSharedRAbundFloatVectors();
 331                 string lastLabel = lookup[0]->getLabel();
 332
 333                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 334                 set<string> processedLabels;
 335                 set<string> userLabels = labels;
 336
 337                 //as long as you are not at the end of the file or done wih the lines you want
 338                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 339
 340                         if (m->control_pressed) {  delete input; return 0;  }
 341
 342                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 343                                 processedLabels.insert(lookup[0]->getLabel());
 344                                 userLabels.erase(lookup[0]->getLabel());
 345
 346                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 347                                 process(lookup);
 348                         }
 349
 350                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 351                                 string saveLabel = lookup[0]->getLabel();
 352
 353                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 354                                 lookup = input->getSharedRAbundFloatVectors(lastLabel);
 355
 356                                 processedLabels.insert(lookup[0]->getLabel());
 357                                 userLabels.erase(lookup[0]->getLabel());
 358
 359                                 //restore real lastlabel to save below
 360                                 lookup[0]->setLabel(saveLabel);
 361
 362                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 363                                 process(lookup);
 364                         }
 365
 366                         lastLabel = lookup[0]->getLabel();
 367
 368                         //get next line to process
 369                         //prevent memory leak
 370                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 371                         lookup = input->getSharedRAbundFloatVectors();
 372                 }
 373
 374
 375                 if (m->control_pressed) { delete input; return 0;  }
 376
 377                 //output error messages about any remaining user labels
 378                 set<string>::iterator it;
 379                 bool needToRun = false;
 380                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 381                         m->mothurOut("Your file does not include the label " + *it);
 382                         if (processedLabels.count(lastLabel) != 1) {
 383                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 384                                 needToRun = true;
 385                         }else {
 386                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 387                         }
 388                 }
 389
 390                 //run last label if you need to
 391                 if (needToRun == true)  {
 392                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } }
 393                         lookup = input->getSharedRAbundFloatVectors(lastLabel);
 394
 395                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 396                         process(lookup);
 397                 }
 398
 399                 delete input;
 400
 401                 return 0;
 402         }
 403         catch(exception& e) {
 404                 m->errorOut(e, "OTUAssociationCommand", "processRelabund");
 405                 exit(1);
 406         }
 407 }
 408 //**********************************************************************************************************************
 409 int OTUAssociationCommand::process(vector<SharedRAbundFloatVector*>& lookup){
 410         try {
 411
 412                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + "." + method + ".otu.corr";
 413                 outputNames.push_back(outputFileName); outputTypes["shared"].push_back(outputFileName);
 414
 415                 ofstream out;
 416                 m->openOutputFile(outputFileName, out);
 417                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 418
 419                 //column headings
 420                 out << "OTUA\tOTUB\t" << method << "Coef\tSignificance\n";
 421
 422                 vector< vector<double> > xy; xy.resize(lookup[0]->getNumBins());
 423                 for (int i = 0; i < lookup[0]->getNumBins(); i++) { for (int j = 0; j < lookup.size(); j++) { xy[i].push_back(lookup[j]->getAbundance(i)); } }
 424
 425                 LinearAlgebra linear;
 426                 for (int i = 0; i < xy.size(); i++) {
 427
 428                         for (int k = 0; k < i; k++) {
 429
 430                                 if (m->control_pressed) { out.close(); return 0; }
 431
 432                                 double coef = 0.0;
 433                                 double sig = 0.0;
 434                                 if (method == "spearman")               {   coef = linear.calcSpearman(xy[i], xy[k], sig);      }
 435                                 else if (method == "pearson")   {       coef = linear.calcPearson(xy[i], xy[k], sig);   }
 436                                 else if (method == "kendall")   {       coef = linear.calcKendall(xy[i], xy[k], sig);   }
 437                                 else { m->mothurOut("[ERROR]: invalid method, choices are spearman, pearson or kendall."); m->mothurOutEndLine(); m->control_pressed = true; }
 438
 439                 out << m->binLabelsInFile[i] << '\t' << m->binLabelsInFile[k] << '\t' << coef << '\t' << sig << endl;
 440                         }
 441                 }
 442
 443                 out.close();
 444
 445                 return 0;
 446
 447         }
 448         catch(exception& e) {
 449                 m->errorOut(e, "OTUAssociationCommand", "process");
 450                 exit(1);
 451         }
 452 }
 453 /*****************************************************************/
 454
 455
 456
 457
 458
 459
 460
 461