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