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