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 pmetadata("metadata", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pmetadata);
  19         CommandParameter pcutoff("cutoff", "Number", "", "10", "", "", "",false,false); parameters.push_back(pcutoff);
  20                 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
  21                 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
  22                 CommandParameter pmethod("method", "Multiple", "pearson-spearman-kendall", "pearson", "", "", "",false,false); parameters.push_back(pmethod);
  23                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  24                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  25
  26                 vector<string> myArray;
  27                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  28                 return myArray;
  29         }
  30         catch(exception& e) {
  31                 m->errorOut(e, "OTUAssociationCommand", "setParameters");
  32                 exit(1);
  33         }
  34 }
  35 //**********************************************************************************************************************
  36 string OTUAssociationCommand::getHelpString(){
  37         try {
  38                 string helpString = "";
  39                 helpString += "The otu.association command reads a shared or relabund file and calculates the correlation coefficients between otus.\n";
  40         helpString += "If you provide a metadata file, mothur will calculate te correlation bewteen the metadata and the otus.\n";
  41                 helpString += "The otu.association command parameters are shared, relabund, metadata, groups, method, cutoff and label.  The shared or relabund parameter is required.\n";
  42                 helpString += "The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes.\n";
  43                 helpString += "The label parameter allows you to select what distances level you would like used, and are also separated by dashes.\n";
  44         helpString += "The cutoff parameter allows you to set a pvalue at which the otu will be reported.\n";
  45                 helpString += "The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n";
  46                 helpString += "The otu.association command should be in the following format: otu.association(shared=yourSharedFile, method=yourMethod).\n";
  47                 helpString += "Example otu.association(shared=genus.pool.shared, method=kendall).\n";
  48                 helpString += "The otu.association command outputs a .otu.corr file.\n";
  49                 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
  50                 return helpString;
  51         }
  52         catch(exception& e) {
  53                 m->errorOut(e, "OTUAssociationCommand", "getHelpString");
  54                 exit(1);
  55         }
  56 }
  57 //**********************************************************************************************************************
  58 string OTUAssociationCommand::getOutputFileNameTag(string type, string inputName=""){
  59         try {
  60         string outputFileName = "";
  61                 map<string, vector<string> >::iterator it;
  62
  63         //is this a type this command creates
  64         it = outputTypes.find(type);
  65         if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
  66         else {
  67             if (type == "otucorr") {  outputFileName =  "otu.corr"; }
  68             else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
  69         }
  70         return outputFileName;
  71         }
  72         catch(exception& e) {
  73                 m->errorOut(e, "OTUAssociationCommand", "getOutputFileNameTag");
  74                 exit(1);
  75         }
  76 }
  77 //**********************************************************************************************************************
  78 OTUAssociationCommand::OTUAssociationCommand(){
  79         try {
  80                 abort = true; calledHelp = true;
  81                 setParameters();
  82                 vector<string> tempOutNames;
  83                 outputTypes["otucorr"] = tempOutNames;
  84         }
  85         catch(exception& e) {
  86                 m->errorOut(e, "OTUAssociationCommand", "OTUAssociationCommand");
  87                 exit(1);
  88         }
  89 }
  90 //**********************************************************************************************************************
  91 OTUAssociationCommand::OTUAssociationCommand(string option)  {
  92         try {
  93                 abort = false; calledHelp = false;
  94                 allLines = 1;
  95
  96                 //allow user to run help
  97                 if(option == "help") { help(); abort = true; calledHelp = true; }
  98                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  99
 100                 else {
 101                         vector<string> myArray = setParameters();
 102
 103                         OptionParser parser(option);
 104                         map<string, string> parameters = parser.getParameters();
 105
 106                         ValidParameters validParameter;
 107                         map<string, string>::iterator it;
 108
 109                         //check to make sure all parameters are valid for command
 110                         for (it = parameters.begin(); it != parameters.end(); it++) {
 111                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 112                         }
 113
 114                         vector<string> tempOutNames;
 115                         outputTypes["otucorr"] = tempOutNames;
 116
 117                         //if the user changes the input directory command factory will send this info to us in the output parameter
 118                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 119                         if (inputDir == "not found"){   inputDir = "";          }
 120                         else {
 121                                 string path;
 122                                 it = parameters.find("shared");
 123                                 //user has given a template file
 124                                 if(it != parameters.end()){
 125                                         path = m->hasPath(it->second);
 126                                         //if the user has not given a path then, add inputdir. else leave path alone.
 127                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 128                                 }
 129
 130                                 it = parameters.find("relabund");
 131                                 //user has given a template file
 132                                 if(it != parameters.end()){
 133                                         path = m->hasPath(it->second);
 134                                         //if the user has not given a path then, add inputdir. else leave path alone.
 135                                         if (path == "") {       parameters["relabund"] = inputDir + it->second;         }
 136                                 }
 137
 138                 it = parameters.find("metadata");
 139                                 //user has given a template file
 140                                 if(it != parameters.end()){
 141                                         path = m->hasPath(it->second);
 142                                         //if the user has not given a path then, add inputdir. else leave path alone.
 143                                         if (path == "") {       parameters["metadata"] = inputDir + it->second;         }
 144                                 }
 145                         }
 146
 147
 148                         //check for required parameters
 149                         sharedfile = validParameter.validFile(parameters, "shared", true);
 150                         if (sharedfile == "not open") { abort = true; }
 151                         else if (sharedfile == "not found") { sharedfile = ""; }
 152                         else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
 153
 154                         relabundfile = validParameter.validFile(parameters, "relabund", true);
 155                         if (relabundfile == "not open") { abort = true; }
 156                         else if (relabundfile == "not found") { relabundfile = ""; }
 157                         else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
 158
 159             metadatafile = validParameter.validFile(parameters, "metadata", true);
 160                         if (metadatafile == "not open") { abort = true; metadatafile = ""; }
 161                         else if (metadatafile == "not found") { metadatafile = ""; }
 162
 163                         groups = validParameter.validFile(parameters, "groups", false);
 164                         if (groups == "not found") { groups = "";  pickedGroups = false;  }
 165                         else {
 166                                 pickedGroups = true;
 167                                 m->splitAtDash(groups, Groups);
 168                         }
 169                         m->setGroups(Groups);
 170
 171                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(inputFileName);  }
 172
 173                         label = validParameter.validFile(parameters, "label", false);
 174                         if (label == "not found") { label = ""; }
 175                         else {
 176                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 177                                 else { allLines = 1;  }
 178                         }
 179
 180                         if ((relabundfile == "") && (sharedfile == "")) {
 181                                 //is there are current file available for any of these?
 182                                 //give priority to shared, then relabund
 183                                 //if there is a current shared file, use it
 184                                 sharedfile = m->getSharedFile();
 185                                 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 186                                 else {
 187                                         relabundfile = m->getRelAbundFile();
 188                                         if (relabundfile != "") { inputFileName = relabundfile;  m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
 189                                         else {
 190                                                 m->mothurOut("You must provide either a shared or relabund file."); m->mothurOutEndLine(); abort = true;
 191                                         }
 192                                 }
 193                         }
 194
 195
 196                         if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared or relabund file."); m->mothurOutEndLine(); abort = true;  }
 197
 198                         method = validParameter.validFile(parameters, "method", false);         if (method == "not found"){     method = "pearson";             }
 199
 200             string temp = validParameter.validFile(parameters, "cutoff", false);
 201                         if (temp == "not found") { temp = "10"; }
 202                         m->mothurConvert(temp, cutoff);
 203
 204                         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; }
 205
 206                 }
 207         }
 208         catch(exception& e) {
 209                 m->errorOut(e, "OTUAssociationCommand", "OTUAssociationCommand");
 210                 exit(1);
 211         }
 212 }
 213 //**********************************************************************************************************************
 214
 215 int OTUAssociationCommand::execute(){
 216         try {
 217
 218                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 219
 220         if (metadatafile != "") {  readMetadata(); }
 221
 222                 //function are identical just different datatypes
 223                 if (sharedfile != "")                   {  processShared();             }
 224                 else if (relabundfile != "")    {  processRelabund();   }
 225
 226                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       m->mothurRemove(outputNames[i]); } return 0; }
 227
 228                 m->mothurOutEndLine();
 229                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 230                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 231                 m->mothurOutEndLine();
 232
 233                 return 0;
 234         }
 235         catch(exception& e) {
 236                 m->errorOut(e, "OTUAssociationCommand", "execute");
 237                 exit(1);
 238         }
 239 }
 240 //**********************************************************************************************************************
 241 int OTUAssociationCommand::processShared(){
 242         try {
 243                 InputData* input = new InputData(sharedfile, "sharedfile");
 244                 vector<SharedRAbundVector*> lookup = input->getSharedRAbundVectors();
 245                 string lastLabel = lookup[0]->getLabel();
 246
 247         if (metadatafile != "") {
 248             getMetadata();
 249             bool error = false;
 250             if (metadata[0].size() != lookup.size()) { m->mothurOut("[ERROR]: You have selected to use " + toString(metadata[0].size()) + " data rows from the metadata file, but " + toString(lookup.size()) + " from the shared file.\n");  m->control_pressed = true; error=true;}
 251             if (error) {
 252                 //maybe add extra info here?? compare groups in each file??
 253             }
 254         }
 255
 256                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 257                 set<string> processedLabels;
 258                 set<string> userLabels = labels;
 259
 260                 //as long as you are not at the end of the file or done wih the lines you want
 261                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 262
 263                         if (m->control_pressed) {  delete input; return 0;  }
 264
 265                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 266                                 processedLabels.insert(lookup[0]->getLabel());
 267                                 userLabels.erase(lookup[0]->getLabel());
 268
 269                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 270                                 process(lookup);
 271                         }
 272
 273                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 274                                 string saveLabel = lookup[0]->getLabel();
 275
 276                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 277                                 lookup = input->getSharedRAbundVectors(lastLabel);
 278
 279                                 processedLabels.insert(lookup[0]->getLabel());
 280                                 userLabels.erase(lookup[0]->getLabel());
 281
 282                                 //restore real lastlabel to save below
 283                                 lookup[0]->setLabel(saveLabel);
 284
 285                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 286                                 process(lookup);
 287                         }
 288
 289                         lastLabel = lookup[0]->getLabel();
 290
 291                         //get next line to process
 292                         //prevent memory leak
 293                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 294                         lookup = input->getSharedRAbundVectors();
 295                 }
 296
 297
 298                 if (m->control_pressed) { delete input; return 0;  }
 299
 300                 //output error messages about any remaining user labels
 301                 set<string>::iterator it;
 302                 bool needToRun = false;
 303                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 304                         m->mothurOut("Your file does not include the label " + *it);
 305                         if (processedLabels.count(lastLabel) != 1) {
 306                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 307                                 needToRun = true;
 308                         }else {
 309                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 310                         }
 311                 }
 312
 313                 //run last label if you need to
 314                 if (needToRun == true)  {
 315                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } }
 316                         lookup = input->getSharedRAbundVectors(lastLabel);
 317
 318                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 319                         process(lookup);
 320                 }
 321
 322                 delete input;
 323
 324                 return 0;
 325         }
 326         catch(exception& e) {
 327                 m->errorOut(e, "OTUAssociationCommand", "processShared");
 328                 exit(1);
 329         }
 330 }
 331 //**********************************************************************************************************************
 332 int OTUAssociationCommand::process(vector<SharedRAbundVector*>& lookup){
 333         try {
 334
 335                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + "." + method + "." + getOutputFileNameTag("otucorr");
 336                 outputNames.push_back(outputFileName); outputTypes["otucorr"].push_back(outputFileName);
 337
 338                 ofstream out;
 339                 m->openOutputFile(outputFileName, out);
 340                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 341
 342                 //column headings
 343                 if (metadatafile == "") { out << "OTUA\tOTUB\t" << method << "Coef\tSignificance\n"; }
 344         else { out << "OTUA\tMetadata\t" << method << "Coef\tSignificance\n";  }
 345
 346
 347                 vector< vector<double> > xy; xy.resize(lookup[0]->getNumBins());
 348                 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)); } }
 349
 350                 LinearAlgebra linear;
 351         if (metadatafile == "") {//compare otus
 352             for (int i = 0; i < xy.size(); i++) {
 353
 354                 for (int k = 0; k < i; k++) {
 355
 356                     if (m->control_pressed) { out.close(); return 0; }
 357
 358                     double coef = 0.0;
 359                     double sig = 0.0;
 360                     if (method == "spearman")           {   coef = linear.calcSpearman(xy[i], xy[k], sig);      }
 361                     else if (method == "pearson")       {       coef = linear.calcPearson(xy[i], xy[k], sig);   }
 362                     else if (method == "kendall")       {       coef = linear.calcKendall(xy[i], xy[k], sig);   }
 363                     else { m->mothurOut("[ERROR]: invalid method, choices are spearman, pearson or kendall."); m->mothurOutEndLine(); m->control_pressed = true; }
 364
 365                     if (sig < cutoff) { out << m->binLabelsInFile[i] << '\t' << m->binLabelsInFile[k] << '\t' << coef << '\t' << sig << endl; }
 366                 }
 367             }
 368                 }else { //compare otus to metadata
 369             for (int i = 0; i < xy.size(); i++) {
 370
 371                 for (int k = 0; k < metadata.size(); k++) {
 372
 373                     if (m->control_pressed) { out.close(); return 0; }
 374
 375                     double coef = 0.0;
 376                     double sig = 0.0;
 377                     if (method == "spearman")           {   coef = linear.calcSpearman(xy[i], metadata[k], sig);        }
 378                     else if (method == "pearson")       {       coef = linear.calcPearson(xy[i], metadata[k], sig);     }
 379                     else if (method == "kendall")       {       coef = linear.calcKendall(xy[i], metadata[k], sig);     }
 380                     else { m->mothurOut("[ERROR]: invalid method, choices are spearman, pearson or kendall."); m->mothurOutEndLine(); m->control_pressed = true; }
 381
 382                     if (sig < cutoff) { out << m->binLabelsInFile[i] << '\t' << metadataLabels[k] << '\t' << coef << '\t' << sig << endl; }
 383                 }
 384             }
 385
 386         }
 387                 out.close();
 388
 389
 390                 return 0;
 391
 392         }
 393         catch(exception& e) {
 394                 m->errorOut(e, "OTUAssociationCommand", "process");
 395                 exit(1);
 396         }
 397 }
 398 //**********************************************************************************************************************
 399 int OTUAssociationCommand::processRelabund(){
 400         try {
 401                 InputData* input = new InputData(relabundfile, "relabund");
 402                 vector<SharedRAbundFloatVector*> lookup = input->getSharedRAbundFloatVectors();
 403                 string lastLabel = lookup[0]->getLabel();
 404
 405         if (metadatafile != "") {
 406             getMetadata();
 407             bool error = false;
 408             if (metadata[0].size() != lookup.size()) { m->mothurOut("[ERROR]: You have selected to use " + toString(metadata[0].size()) + " data rows from the metadata file, but " + toString(lookup.size()) + " from the relabund file.\n");  m->control_pressed = true; error=true;}
 409             if (error) {
 410                 //maybe add extra info here?? compare groups in each file??
 411             }
 412         }
 413
 414
 415
 416                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 417                 set<string> processedLabels;
 418                 set<string> userLabels = labels;
 419
 420                 //as long as you are not at the end of the file or done wih the lines you want
 421                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 422
 423                         if (m->control_pressed) {  delete input; return 0;  }
 424
 425                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 426                                 processedLabels.insert(lookup[0]->getLabel());
 427                                 userLabels.erase(lookup[0]->getLabel());
 428
 429                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 430                                 process(lookup);
 431                         }
 432
 433                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 434                                 string saveLabel = lookup[0]->getLabel();
 435
 436                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 437                                 lookup = input->getSharedRAbundFloatVectors(lastLabel);
 438
 439                                 processedLabels.insert(lookup[0]->getLabel());
 440                                 userLabels.erase(lookup[0]->getLabel());
 441
 442                                 //restore real lastlabel to save below
 443                                 lookup[0]->setLabel(saveLabel);
 444
 445                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 446                                 process(lookup);
 447                         }
 448
 449                         lastLabel = lookup[0]->getLabel();
 450
 451                         //get next line to process
 452                         //prevent memory leak
 453                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 454                         lookup = input->getSharedRAbundFloatVectors();
 455                 }
 456
 457
 458                 if (m->control_pressed) { delete input; return 0;  }
 459
 460                 //output error messages about any remaining user labels
 461                 set<string>::iterator it;
 462                 bool needToRun = false;
 463                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 464                         m->mothurOut("Your file does not include the label " + *it);
 465                         if (processedLabels.count(lastLabel) != 1) {
 466                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 467                                 needToRun = true;
 468                         }else {
 469                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 470                         }
 471                 }
 472
 473                 //run last label if you need to
 474                 if (needToRun == true)  {
 475                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } }
 476                         lookup = input->getSharedRAbundFloatVectors(lastLabel);
 477
 478                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 479                         process(lookup);
 480                 }
 481
 482                 delete input;
 483
 484                 return 0;
 485         }
 486         catch(exception& e) {
 487                 m->errorOut(e, "OTUAssociationCommand", "processRelabund");
 488                 exit(1);
 489         }
 490 }
 491 //**********************************************************************************************************************
 492 int OTUAssociationCommand::process(vector<SharedRAbundFloatVector*>& lookup){
 493         try {
 494
 495                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + "." + method + "." + getOutputFileNameTag("otucorr");
 496                 outputNames.push_back(outputFileName); outputTypes["otucorr"].push_back(outputFileName);
 497
 498                 ofstream out;
 499                 m->openOutputFile(outputFileName, out);
 500                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 501
 502                 //column headings
 503                 if (metadatafile == "") { out << "OTUA\tOTUB\t" << method << "Coef\tSignificance\n"; }
 504         else { out << "OTUA\tMetadata\t" << method << "Coef\tSignificance\n";  }
 505
 506                 vector< vector<double> > xy; xy.resize(lookup[0]->getNumBins());
 507                 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)); } }
 508
 509                 LinearAlgebra linear;
 510         if (metadatafile == "") {//compare otus
 511             for (int i = 0; i < xy.size(); i++) {
 512
 513                 for (int k = 0; k < i; k++) {
 514
 515                     if (m->control_pressed) { out.close(); return 0; }
 516
 517                     double coef = 0.0;
 518                     double sig = 0.0;
 519                     if (method == "spearman")           {   coef = linear.calcSpearman(xy[i], xy[k], sig);      }
 520                     else if (method == "pearson")       {       coef = linear.calcPearson(xy[i], xy[k], sig);   }
 521                     else if (method == "kendall")       {       coef = linear.calcKendall(xy[i], xy[k], sig);   }
 522                     else { m->mothurOut("[ERROR]: invalid method, choices are spearman, pearson or kendall."); m->mothurOutEndLine(); m->control_pressed = true; }
 523
 524                     if (sig < cutoff) { out << m->binLabelsInFile[i] << '\t' << m->binLabelsInFile[k] << '\t' << coef << '\t' << sig << endl; }
 525                 }
 526             }
 527                 }else { //compare otus to metadata
 528             for (int i = 0; i < xy.size(); i++) {
 529
 530                 for (int k = 0; k < metadata.size(); k++) {
 531
 532                     if (m->control_pressed) { out.close(); return 0; }
 533
 534                     double coef = 0.0;
 535                     double sig = 0.0;
 536                     if (method == "spearman")           {   coef = linear.calcSpearman(xy[i], metadata[k], sig);        }
 537                     else if (method == "pearson")       {       coef = linear.calcPearson(xy[i], metadata[k], sig);     }
 538                     else if (method == "kendall")       {       coef = linear.calcKendall(xy[i], metadata[k], sig);     }
 539                     else { m->mothurOut("[ERROR]: invalid method, choices are spearman, pearson or kendall."); m->mothurOutEndLine(); m->control_pressed = true; }
 540
 541                     if (sig < cutoff) { out << m->binLabelsInFile[i] << '\t' << metadataLabels[k] << '\t' << coef << '\t' << sig << endl; }
 542                 }
 543             }
 544
 545         }
 546
 547                 out.close();
 548
 549                 return 0;
 550
 551         }
 552         catch(exception& e) {
 553                 m->errorOut(e, "OTUAssociationCommand", "process");
 554                 exit(1);
 555         }
 556 }
 557 /*****************************************************************/
 558 int OTUAssociationCommand::readMetadata(){
 559         try {
 560                 ifstream in;
 561                 m->openInputFile(metadatafile, in);
 562
 563                 string headerLine = m->getline(in); m->gobble(in);
 564                 istringstream iss (headerLine,istringstream::in);
 565
 566                 //read the first label, because it refers to the groups
 567                 string columnLabel;
 568                 iss >> columnLabel; m->gobble(iss);
 569
 570                 //save names of columns you are reading
 571                 while (!iss.eof()) {
 572                         iss >> columnLabel; m->gobble(iss);
 573                         metadataLabels.push_back(columnLabel);
 574                 }
 575                 int count = metadataLabels.size();
 576
 577                 //read rest of file
 578                 while (!in.eof()) {
 579
 580                         if (m->control_pressed) { in.close(); return 0; }
 581
 582                         string group = "";
 583                         in >> group; m->gobble(in);
 584
 585                         SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
 586                         tempLookup->setGroup(group);
 587                         tempLookup->setLabel("1");
 588
 589                         for (int i = 0; i < count; i++) {
 590                                 float temp = 0.0;
 591                                 in >> temp;
 592                                 tempLookup->push_back(temp, group);
 593                         }
 594
 595                         metadataLookup.push_back(tempLookup);
 596
 597                         m->gobble(in);
 598                 }
 599                 in.close();
 600
 601                 return 0;
 602         }
 603         catch(exception& e) {
 604                 m->errorOut(e, "OTUAssociationCommand", "readMetadata");
 605                 exit(1);
 606         }
 607 }
 608 /*****************************************************************/
 609 //eliminate groups user did not pick, remove zeroed out otus, fill metadata vector.
 610 int OTUAssociationCommand::getMetadata(){
 611         try {
 612
 613                 vector<string> mGroups = m->getGroups();
 614
 615                 bool remove = false;
 616                 for (int i = 0; i < metadataLookup.size(); i++) {
 617                         //if this sharedrabund is not from a group the user wants then delete it.
 618                         if (!(m->inUsersGroups(metadataLookup[i]->getGroup(), mGroups))) {
 619                                 delete metadataLookup[i]; metadataLookup[i] = NULL;
 620                                 metadataLookup.erase(metadataLookup.begin()+i);
 621                                 i--;
 622                                 remove = true;
 623                         }
 624                 }
 625
 626         vector<SharedRAbundFloatVector*> newLookup;
 627                 for (int i = 0; i < metadataLookup.size(); i++) {
 628                         SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
 629                         temp->setLabel(metadataLookup[i]->getLabel());
 630                         temp->setGroup(metadataLookup[i]->getGroup());
 631                         newLookup.push_back(temp);
 632                 }
 633
 634                 //for each bin
 635         vector<string> newBinLabels;
 636                 for (int i = 0; i < metadataLookup[0]->getNumBins(); i++) {
 637                         if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
 638
 639                         //look at each sharedRabund and make sure they are not all zero
 640                         bool allZero = true;
 641                         for (int j = 0; j < metadataLookup.size(); j++) {
 642                                 if (metadataLookup[j]->getAbundance(i) != 0) { allZero = false;  break;  }
 643                         }
 644
 645                         //if they are not all zero add this bin
 646                         if (!allZero) {
 647                                 for (int j = 0; j < metadataLookup.size(); j++) {
 648                                         newLookup[j]->push_back(metadataLookup[j]->getAbundance(i), metadataLookup[j]->getGroup());
 649                                 }
 650                 newBinLabels.push_back(metadataLabels[i]);
 651                         }
 652                 }
 653
 654         metadataLabels = newBinLabels;
 655
 656                 for (int j = 0; j < metadataLookup.size(); j++) {  delete metadataLookup[j];  }
 657         metadataLookup.clear();
 658
 659         metadata.resize(newLookup[0]->getNumBins());
 660                 for (int i = 0; i < newLookup[0]->getNumBins(); i++) { for (int j = 0; j < newLookup.size(); j++) { metadata[i].push_back(newLookup[j]->getAbundance(i)); } }
 661
 662         for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  }
 663
 664         return 0;
 665     }
 666         catch(exception& e) {
 667                 m->errorOut(e, "OTUAssociationCommand", "getMetadata");
 668                 exit(1);
 669         }
 670 }
 671 /*****************************************************************/
 672
 673
 674
 675
 676
 677
 678
 679