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                 return 0;
 318
 319         }
 320         catch(exception& e) {
 321                 m->errorOut(e, "OTUAssociationCommand", "process");
 322                 exit(1);
 323         }
 324 }
 325 //**********************************************************************************************************************
 326 int OTUAssociationCommand::processRelabund(){
 327         try {
 328                 InputData* input = new InputData(relabundfile, "relabund");
 329                 vector<SharedRAbundFloatVector*> lookup = input->getSharedRAbundFloatVectors();
 330                 string lastLabel = lookup[0]->getLabel();
 331
 332                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 333                 set<string> processedLabels;
 334                 set<string> userLabels = labels;
 335
 336                 //as long as you are not at the end of the file or done wih the lines you want
 337                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 338
 339                         if (m->control_pressed) {  delete input; return 0;  }
 340
 341                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 342                                 processedLabels.insert(lookup[0]->getLabel());
 343                                 userLabels.erase(lookup[0]->getLabel());
 344
 345                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 346                                 process(lookup);
 347                         }
 348
 349                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 350                                 string saveLabel = lookup[0]->getLabel();
 351
 352                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 353                                 lookup = input->getSharedRAbundFloatVectors(lastLabel);
 354
 355                                 processedLabels.insert(lookup[0]->getLabel());
 356                                 userLabels.erase(lookup[0]->getLabel());
 357
 358                                 //restore real lastlabel to save below
 359                                 lookup[0]->setLabel(saveLabel);
 360
 361                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 362                                 process(lookup);
 363                         }
 364
 365                         lastLabel = lookup[0]->getLabel();
 366
 367                         //get next line to process
 368                         //prevent memory leak
 369                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 370                         lookup = input->getSharedRAbundFloatVectors();
 371                 }
 372
 373
 374                 if (m->control_pressed) { delete input; return 0;  }
 375
 376                 //output error messages about any remaining user labels
 377                 set<string>::iterator it;
 378                 bool needToRun = false;
 379                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 380                         m->mothurOut("Your file does not include the label " + *it);
 381                         if (processedLabels.count(lastLabel) != 1) {
 382                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 383                                 needToRun = true;
 384                         }else {
 385                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 386                         }
 387                 }
 388
 389                 //run last label if you need to
 390                 if (needToRun == true)  {
 391                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } }
 392                         lookup = input->getSharedRAbundFloatVectors(lastLabel);
 393
 394                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 395                         process(lookup);
 396                 }
 397
 398                 delete input;
 399
 400                 return 0;
 401         }
 402         catch(exception& e) {
 403                 m->errorOut(e, "OTUAssociationCommand", "processRelabund");
 404                 exit(1);
 405         }
 406 }
 407 //**********************************************************************************************************************
 408 int OTUAssociationCommand::process(vector<SharedRAbundFloatVector*>& lookup){
 409         try {
 410
 411                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + "." + method + ".otu.corr";
 412                 outputNames.push_back(outputFileName); outputTypes["shared"].push_back(outputFileName);
 413
 414                 ofstream out;
 415                 m->openOutputFile(outputFileName, out);
 416                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 417
 418                 //column headings
 419                 out << "OTUA\tOTUB\t" << method << "Coef\tSignificance\n";
 420
 421                 vector< vector<double> > xy; xy.resize(lookup[0]->getNumBins());
 422                 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)); } }
 423
 424                 LinearAlgebra linear;
 425                 for (int i = 0; i < xy.size(); i++) {
 426
 427                         for (int k = 0; k < i; k++) {
 428
 429                                 if (m->control_pressed) { out.close(); return 0; }
 430
 431                                 double coef = 0.0;
 432                                 double sig = 0.0;
 433                                 if (method == "spearman")               {   coef = linear.calcSpearman(xy[i], xy[k], sig);      }
 434                                 else if (method == "pearson")   {       coef = linear.calcPearson(xy[i], xy[k], sig);   }
 435                                 else if (method == "kendall")   {       coef = linear.calcKendall(xy[i], xy[k], sig);   }
 436                                 else { m->mothurOut("[ERROR]: invalid method, choices are spearman, pearson or kendall."); m->mothurOutEndLine(); m->control_pressed = true; }
 437
 438                                 out << m->binLabelsInFile[i] << '\t' << m->binLabelsInFile[k] << coef << '\t' << sig << endl;
 439                         }
 440                 }
 441
 442                 out.close();
 443
 444                 return 0;
 445
 446         }
 447         catch(exception& e) {
 448                 m->errorOut(e, "OTUAssociationCommand", "process");
 449                 exit(1);
 450         }
 451 }
 452 /*****************************************************************/
 453
 454
 455
 456
 457
 458
 459
 460