corraxescommand.cpp

   1 /*
   2  *  corraxescommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 12/22/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "corraxescommand.h"
  11 #include "sharedutilities.h"
  12
  13 //********************************************************************************************************************
  14 //sorts highes to lowest
  15 inline bool compareSpearman(spearmanRank left, spearmanRank right){
  16         return (left.score > right.score);
  17 }
  18 //**********************************************************************************************************************
  19 vector<string> CorrAxesCommand::getValidParameters(){
  20         try {
  21                 string Array[] =  {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
  22                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  23                 return myArray;
  24         }
  25         catch(exception& e) {
  26                 m->errorOut(e, "CorrAxesCommand", "getValidParameters");
  27                 exit(1);
  28         }
  29 }
  30 //**********************************************************************************************************************
  31 vector<string> CorrAxesCommand::getRequiredParameters(){
  32         try {
  33                 string Array[] =  {"axes"};
  34                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  35                 return myArray;
  36         }
  37         catch(exception& e) {
  38                 m->errorOut(e, "CorrAxesCommand", "getRequiredParameters");
  39                 exit(1);
  40         }
  41 }
  42 //**********************************************************************************************************************
  43 CorrAxesCommand::CorrAxesCommand(){
  44         try {
  45                 abort = true;
  46                 //initialize outputTypes
  47                 vector<string> tempOutNames;
  48                 outputTypes["corr.axes"] = tempOutNames;
  49         }
  50         catch(exception& e) {
  51                 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
  52                 exit(1);
  53         }
  54 }
  55
  56 //**********************************************************************************************************************
  57 vector<string> CorrAxesCommand::getRequiredFiles(){
  58         try {
  59                 vector<string> myArray;
  60                 return myArray;
  61         }
  62         catch(exception& e) {
  63                 m->errorOut(e, "CorrAxesCommand", "getRequiredFiles");
  64                 exit(1);
  65         }
  66 }
  67 //**********************************************************************************************************************
  68 CorrAxesCommand::CorrAxesCommand(string option)  {
  69         try {
  70                 abort = false;
  71                 globaldata = GlobalData::getInstance();
  72
  73                 //allow user to run help
  74                 if(option == "help") { help(); abort = true; }
  75
  76                 else {
  77                         //valid paramters for this command
  78                         string Array[] =  {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
  79                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  80
  81                         OptionParser parser(option);
  82                         map<string, string> parameters = parser.getParameters();
  83
  84                         ValidParameters validParameter;
  85                         map<string, string>::iterator it;
  86
  87                         //check to make sure all parameters are valid for command
  88                         for (it = parameters.begin(); it != parameters.end(); it++) {
  89                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  90                         }
  91
  92                         vector<string> tempOutNames;
  93                         outputTypes["corr.axes"] = tempOutNames;
  94
  95                         //if the user changes the input directory command factory will send this info to us in the output parameter
  96                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  97                         if (inputDir == "not found"){   inputDir = "";          }
  98                         else {
  99                                 string path;
 100                                 it = parameters.find("axes");
 101                                 //user has given a template file
 102                                 if(it != parameters.end()){
 103                                         path = m->hasPath(it->second);
 104                                         //if the user has not given a path then, add inputdir. else leave path alone.
 105                                         if (path == "") {       parameters["axes"] = inputDir + it->second;             }
 106                                 }
 107
 108                                 it = parameters.find("shared");
 109                                 //user has given a template file
 110                                 if(it != parameters.end()){
 111                                         path = m->hasPath(it->second);
 112                                         //if the user has not given a path then, add inputdir. else leave path alone.
 113                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 114                                 }
 115
 116                                 it = parameters.find("relabund");
 117                                 //user has given a template file
 118                                 if(it != parameters.end()){
 119                                         path = m->hasPath(it->second);
 120                                         //if the user has not given a path then, add inputdir. else leave path alone.
 121                                         if (path == "") {       parameters["relabund"] = inputDir + it->second;         }
 122                                 }
 123
 124                                 it = parameters.find("metadata");
 125                                 //user has given a template file
 126                                 if(it != parameters.end()){
 127                                         path = m->hasPath(it->second);
 128                                         //if the user has not given a path then, add inputdir. else leave path alone.
 129                                         if (path == "") {       parameters["metadata"] = inputDir + it->second;         }
 130                                 }
 131                         }
 132
 133
 134                         //check for required parameters
 135                         axesfile = validParameter.validFile(parameters, "axes", true);
 136                         if (axesfile == "not open") { abort = true; }
 137                         else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true;  }
 138
 139                         sharedfile = validParameter.validFile(parameters, "shared", true);
 140                         if (sharedfile == "not open") { abort = true; }
 141                         else if (sharedfile == "not found") { sharedfile = ""; }
 142                         else { inputFileName = sharedfile; }
 143
 144                         relabundfile = validParameter.validFile(parameters, "relabund", true);
 145                         if (relabundfile == "not open") { abort = true; }
 146                         else if (relabundfile == "not found") { relabundfile = ""; }
 147                         else { inputFileName = relabundfile; }
 148
 149                         metadatafile = validParameter.validFile(parameters, "metadata", true);
 150                         if (metadatafile == "not open") { abort = true; }
 151                         else if (metadatafile == "not found") { metadatafile = ""; }
 152                         else { inputFileName = metadatafile; }
 153
 154                         groups = validParameter.validFile(parameters, "groups", false);
 155                         if (groups == "not found") { groups = "";  pickedGroups = false;  }
 156                         else {
 157                                 pickedGroups = true;
 158                                 m->splitAtDash(groups, Groups);
 159                         }
 160                         globaldata->Groups = Groups;
 161
 162                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(inputFileName);  }
 163
 164                         label = validParameter.validFile(parameters, "label", false);
 165                         if (label == "not found") { label = ""; m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); label=""; }
 166
 167                         if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) { m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true;  }
 168
 169                         if (metadatafile != "") {
 170                                 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true;  }
 171                         }else {
 172                                 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true;  }
 173                         }
 174                         string temp;
 175                         temp = validParameter.validFile(parameters, "numaxes", false);          if (temp == "not found"){       temp = "3";                             }
 176                         convert(temp, numaxes);
 177
 178                         method = validParameter.validFile(parameters, "method", false);         if (method == "not found"){     method = "pearson";             }
 179
 180                         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; }
 181
 182                 }
 183         }
 184         catch(exception& e) {
 185                 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
 186                 exit(1);
 187         }
 188 }
 189 //**********************************************************************************************************************
 190
 191 void CorrAxesCommand::help(){
 192         try {
 193                 m->mothurOut("The corr.axes command reads a shared or relabund file as well as a pcoa file and calculates the correlation coefficient.\n");
 194                 m->mothurOut("The corr.axes command parameters are shared, relabund, axes, metadata, groups, method, numaxes and label.  The shared, relabund or metadata and axes parameters are required.  If shared is given the relative abundance is calculated.\n");
 195                 m->mothurOut("The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes.\n");
 196                 m->mothurOut("The label parameter allows you to select what distance level you would like used, if none is given the first distance is used.\n");
 197                 m->mothurOut("The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n");
 198                 m->mothurOut("The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n");
 199                 m->mothurOut("The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n");
 200                 m->mothurOut("Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n");
 201                 m->mothurOut("The corr.axes command outputs a .corr.axes file.\n");
 202                 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
 203         }
 204         catch(exception& e) {
 205                 m->errorOut(e, "CorrAxesCommand", "help");
 206                 exit(1);
 207         }
 208 }
 209
 210 //**********************************************************************************************************************
 211
 212 CorrAxesCommand::~CorrAxesCommand(){}
 213
 214 //**********************************************************************************************************************
 215
 216 int CorrAxesCommand::execute(){
 217         try {
 218
 219                 if (abort == true) { return 0; }
 220
 221                 /*************************************************************************************/
 222                 // use smart distancing to get right sharedRabund and convert to relabund if needed  //
 223                 /************************************************************************************/
 224                 if (sharedfile != "") {
 225                         InputData* input = new InputData(sharedfile, "sharedfile");
 226                         getSharedFloat(input);
 227                         delete input;
 228
 229                         if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
 230                         if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
 231
 232                 }else if (relabundfile != "") {
 233                         InputData* input = new InputData(relabundfile, "relabund");
 234                         getSharedFloat(input);
 235                         delete input;
 236
 237                         if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
 238                         if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
 239
 240                 }else if (metadatafile != "") {
 241                         getMetadata();  //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
 242                         if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
 243                         if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
 244
 245                         if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
 246
 247                 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
 248
 249                 if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
 250
 251                 //this is for a sanity check to make sure the axes file and shared file match
 252                 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
 253
 254                 /*************************************************************************************/
 255                 // read axes file and check for file mismatches with shared or relabund file         //
 256                 /************************************************************************************/
 257
 258                 //read axes file
 259                 map<string, vector<float> > axes = readAxes();
 260
 261                 if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
 262
 263                 //sanity check, the read only adds groups that are in the shared or relabund file, but we want to make sure the axes file isn't missing anyone
 264                 if (axes.size() != lookupFloat.size()) {
 265                         map<string, vector<float> >::iterator it;
 266                         for (int i = 0; i < lookupFloat.size(); i++) {
 267                                 it = axes.find(lookupFloat[i]->getGroup());
 268                                 if (it == axes.end()) { m->mothurOut(lookupFloat[i]->getGroup() + " is in your shared of relabund file but not in your axes file, please correct."); m->mothurOutEndLine(); }
 269                         }
 270                         m->control_pressed = true;
 271                 }
 272
 273                 if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
 274
 275                 /*************************************************************************************/
 276                 // calc the r values                                                                                                                            //
 277                 /************************************************************************************/
 278
 279                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
 280                 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
 281                 ofstream out;
 282                 m->openOutputFile(outputFileName, out);
 283                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 284
 285                 //output headings
 286                 if (metadatafile == "") {  out << "OTU\t";      }
 287                 else {  out << "Feature\t";                                             }
 288
 289                 for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; }
 290                 out << endl;
 291
 292                 if (method == "pearson")                {  calcPearson(axes, out);      }
 293                 else if (method == "spearman")  {  calcSpearman(axes, out); }
 294                 else if (method == "kendall")   {  calcKendall(axes, out);      }
 295                 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
 296
 297                 out.close();
 298                 for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 299
 300                 if (m->control_pressed) {  return 0; }
 301
 302                 m->mothurOutEndLine();
 303                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 304                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 305                 m->mothurOutEndLine();
 306
 307                 return 0;
 308         }
 309         catch(exception& e) {
 310                 m->errorOut(e, "CorrAxesCommand", "execute");
 311                 exit(1);
 312         }
 313 }
 314 //**********************************************************************************************************************
 315 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
 316    try {
 317
 318            //find average of each axis - X
 319            vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
 320            for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
 321                    vector<float> temp = it->second;
 322                    for (int i = 0; i < temp.size(); i++) {
 323                            averageAxes[i] += temp[i];
 324                    }
 325            }
 326
 327            for (int i = 0; i < averageAxes.size(); i++) {  averageAxes[i] = averageAxes[i] / (float) axes.size(); }
 328
 329            //for each otu
 330            for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
 331
 332                    if (metadatafile == "") {  out << i+1 << '\t';       }
 333                    else {  out << metadataLabels[i] << '\t';            }
 334
 335                    //find the averages this otu - Y
 336                    float sumOtu = 0.0;
 337                    for (int j = 0; j < lookupFloat.size(); j++) {
 338                            sumOtu += lookupFloat[j]->getAbundance(i);
 339                    }
 340                    float Ybar = sumOtu / (float) lookupFloat.size();
 341
 342                    //find r value for each axis
 343                    for (int k = 0; k < averageAxes.size(); k++) {
 344
 345                            double r = 0.0;
 346                            double numerator = 0.0;
 347                            double denomTerm1 = 0.0;
 348                            double denomTerm2 = 0.0;
 349                            for (int j = 0; j < lookupFloat.size(); j++) {
 350                                    float Yi = lookupFloat[j]->getAbundance(i);
 351                                    float Xi = axes[lookupFloat[j]->getGroup()][k];
 352
 353                                    numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
 354                                    denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
 355                                    denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
 356                            }
 357
 358                            double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
 359
 360                            r = numerator / denom;
 361
 362                            out << r << '\t';
 363                    }
 364
 365                    out << endl;
 366            }
 367
 368            return 0;
 369    }
 370    catch(exception& e) {
 371            m->errorOut(e, "CorrAxesCommand", "calcPearson");
 372            exit(1);
 373    }
 374 }
 375 //**********************************************************************************************************************
 376 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
 377         try {
 378
 379                 //format data
 380                 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
 381                 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
 382                         vector<float> temp = it->second;
 383                         for (int i = 0; i < temp.size(); i++) {
 384                                 spearmanRank member(it->first, temp[i]);
 385                                 scores[i].push_back(member);
 386                         }
 387                 }
 388
 389                 //sort each axis
 390                 for (int i = 0; i < numaxes; i++) {  sort(scores[i].begin(), scores[i].end(), compareSpearman); }
 391
 392                 //find ranks of xi in each axis
 393                 map<string, vector<float> > rankAxes;
 394                 for (int i = 0; i < numaxes; i++) {
 395
 396                         vector<spearmanRank> ties;
 397                         int rankTotal = 0;
 398                         for (int j = 0; j < scores[i].size(); j++) {
 399                                 rankTotal += j;
 400                                 ties.push_back(scores[i][j]);
 401
 402                                 if (j != scores.size()-1) { // you are not the last so you can look ahead
 403                                         if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
 404                                                 for (int k = 0; k < ties.size(); k++) {
 405                                                         float thisrank = rankTotal / (float) ties.size();
 406                                                         rankAxes[ties[k].name].push_back(thisrank);
 407                                                 }
 408                                                 ties.clear();
 409                                                 rankTotal = 0;
 410                                         }
 411                                 }else { // you are the last one
 412                                         for (int k = 0; k < ties.size(); k++) {
 413                                                 float thisrank = rankTotal / (float) ties.size();
 414                                                 rankAxes[ties[k].name].push_back(thisrank);
 415                                         }
 416                                 }
 417                         }
 418                 }
 419
 420
 421
 422                 //for each otu
 423                 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
 424
 425                         if (metadatafile == "") {  out << i+1 << '\t';  }
 426                         else {  out << metadataLabels[i] << '\t';               }
 427
 428                         //find the ranks of this otu - Y
 429                         vector<spearmanRank> otuScores;
 430                         for (int j = 0; j < lookupFloat.size(); j++) {
 431                                 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
 432                                 otuScores.push_back(member);
 433                         }
 434
 435                         sort(otuScores.begin(), otuScores.end(), compareSpearman);
 436
 437                         map<string, float> rankOtus;
 438                         vector<spearmanRank> ties;
 439                         int rankTotal = 0;
 440                         for (int j = 0; j < otuScores.size(); j++) {
 441                                 rankTotal += j;
 442                                 ties.push_back(otuScores[j]);
 443
 444                                 if (j != scores.size()-1) { // you are not the last so you can look ahead
 445                                         if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
 446                                                 for (int k = 0; k < ties.size(); k++) {
 447                                                         float thisrank = rankTotal / (float) ties.size();
 448                                                         rankOtus[ties[k].name] = thisrank;
 449                                                 }
 450                                                 ties.clear();
 451                                                 rankTotal = 0;
 452                                         }
 453                                 }else { // you are the last one
 454                                         for (int k = 0; k < ties.size(); k++) {
 455                                                 float thisrank = rankTotal / (float) ties.size();
 456                                                 rankOtus[ties[k].name] = thisrank;
 457                                         }
 458                                 }
 459                         }
 460
 461                         //calc spearman ranks for each axis for this otu
 462                         for (int j = 0; j < numaxes; j++) {
 463
 464                                 double di = 0.0;
 465                                 for (int k = 0; k < lookupFloat.size(); k++) {
 466
 467                                         float xi = rankAxes[lookupFloat[k]->getGroup()][j];
 468                                         float yi = rankOtus[lookupFloat[k]->getGroup()];
 469
 470                                         di += ((xi - yi) * (xi - yi));
 471                                 }
 472
 473                                 int n = lookupFloat.size();
 474                                 double p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1)));
 475
 476                                 out << p << '\t';
 477                         }
 478
 479
 480                         out << endl;
 481                 }
 482
 483                 return 0;
 484         }
 485         catch(exception& e) {
 486                 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
 487                 exit(1);
 488         }
 489 }
 490 //**********************************************************************************************************************
 491 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
 492         try {
 493
 494                 //format data
 495                 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
 496                 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
 497                         vector<float> temp = it->second;
 498                         for (int i = 0; i < temp.size(); i++) {
 499                                 spearmanRank member(it->first, temp[i]);
 500                                 scores[i].push_back(member);
 501                         }
 502                 }
 503
 504                 //sort each axis
 505                 for (int i = 0; i < numaxes; i++) {  sort(scores[i].begin(), scores[i].end(), compareSpearman); }
 506
 507                 //convert scores to ranks of xi in each axis
 508                 for (int i = 0; i < numaxes; i++) {
 509
 510                         vector<spearmanRank*> ties;
 511                         int rankTotal = 0;
 512                         for (int j = 0; j < scores[i].size(); j++) {
 513                                 rankTotal += j;
 514                                 ties.push_back(&(scores[i][j]));
 515
 516                                 if (j != scores.size()-1) { // you are not the last so you can look ahead
 517                                         if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
 518                                                 for (int k = 0; k < ties.size(); k++) {
 519                                                         float thisrank = rankTotal / (float) ties.size();
 520                                                         (*ties[k]).score = thisrank;
 521                                                 }
 522                                                 ties.clear();
 523                                                 rankTotal = 0;
 524                                         }
 525                                 }else { // you are the last one
 526                                         for (int k = 0; k < ties.size(); k++) {
 527                                                 float thisrank = rankTotal / (float) ties.size();
 528                                                 (*ties[k]).score = thisrank;
 529                                         }
 530                                 }
 531                         }
 532                 }
 533
 534                 //for each otu
 535                 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
 536
 537                         if (metadatafile == "") {  out << i+1 << '\t';  }
 538                         else {  out << metadataLabels[i] << '\t';               }
 539
 540                         //find the ranks of this otu - Y
 541                         vector<spearmanRank> otuScores;
 542                         for (int j = 0; j < lookupFloat.size(); j++) {
 543                                 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
 544                                 otuScores.push_back(member);
 545                         }
 546
 547                         sort(otuScores.begin(), otuScores.end(), compareSpearman);
 548
 549                         map<string, float> rankOtus;
 550                         vector<spearmanRank> ties;
 551                         int rankTotal = 0;
 552                         for (int j = 0; j < otuScores.size(); j++) {
 553                                 rankTotal += j;
 554                                 ties.push_back(otuScores[j]);
 555
 556                                 if (j != scores.size()-1) { // you are not the last so you can look ahead
 557                                         if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
 558                                                 for (int k = 0; k < ties.size(); k++) {
 559                                                         float thisrank = rankTotal / (float) ties.size();
 560                                                         rankOtus[ties[k].name] = thisrank;
 561                                                 }
 562                                                 ties.clear();
 563                                                 rankTotal = 0;
 564                                         }
 565                                 }else { // you are the last one
 566                                         for (int k = 0; k < ties.size(); k++) {
 567                                                 float thisrank = rankTotal / (float) ties.size();
 568                                                 rankOtus[ties[k].name] = thisrank;
 569                                         }
 570                                 }
 571                         }
 572
 573                         //calc spearman ranks for each axis for this otu
 574                         for (int j = 0; j < numaxes; j++) {
 575
 576                                 int P = 0;
 577                                 //assemble otus ranks in same order as axis ranks
 578                                 vector<spearmanRank> otus;
 579                                 for (int l = 0; l < scores[j].size(); l++) {
 580                                         spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
 581                                         otus.push_back(member);
 582                                 }
 583
 584                                 for (int l = 0; l < scores[j].size(); l++) {
 585
 586                                         int numWithHigherRank = 0;
 587                                         float thisrank = otus[l].score;
 588
 589                                         for (int u = l; u < scores[j].size(); u++) {
 590                                                 if (otus[u].score > thisrank) { numWithHigherRank++; }
 591                                         }
 592
 593                                         P += numWithHigherRank;
 594                                 }
 595
 596                                 int n = lookupFloat.size();
 597
 598                                 double p = ( (4 * P) / (float) (n * (n - 1)) ) - 1.0;
 599
 600                                 out << p << '\t';
 601                         }
 602
 603                         out << endl;
 604                 }
 605
 606                 return 0;
 607         }
 608         catch(exception& e) {
 609                 m->errorOut(e, "CorrAxesCommand", "calcKendall");
 610                 exit(1);
 611         }
 612 }
 613 //**********************************************************************************************************************
 614 int CorrAxesCommand::getSharedFloat(InputData* input){
 615         try {
 616                 lookupFloat = input->getSharedRAbundFloatVectors();
 617                 string lastLabel = lookupFloat[0]->getLabel();
 618
 619                 if (label == "") { label = lastLabel;  return 0; }
 620
 621                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 622                 set<string> labels; labels.insert(label);
 623                 set<string> processedLabels;
 624                 set<string> userLabels = labels;
 625
 626                 //as long as you are not at the end of the file or done wih the lines you want
 627                 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
 628
 629                         if (m->control_pressed) {  return 0;  }
 630
 631                         if(labels.count(lookupFloat[0]->getLabel()) == 1){
 632                                 processedLabels.insert(lookupFloat[0]->getLabel());
 633                                 userLabels.erase(lookupFloat[0]->getLabel());
 634                                 break;
 635                         }
 636
 637                         if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 638                                 string saveLabel = lookupFloat[0]->getLabel();
 639
 640                                 for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 641                                 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 642
 643                                 processedLabels.insert(lookupFloat[0]->getLabel());
 644                                 userLabels.erase(lookupFloat[0]->getLabel());
 645
 646                                 //restore real lastlabel to save below
 647                                 lookupFloat[0]->setLabel(saveLabel);
 648                                 break;
 649                         }
 650
 651                         lastLabel = lookupFloat[0]->getLabel();
 652
 653                         //get next line to process
 654                         //prevent memory leak
 655                         for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 656                         lookupFloat = input->getSharedRAbundFloatVectors();
 657                 }
 658
 659
 660                 if (m->control_pressed) { return 0;  }
 661
 662                 //output error messages about any remaining user labels
 663                 set<string>::iterator it;
 664                 bool needToRun = false;
 665                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 666                         m->mothurOut("Your file does not include the label " + *it);
 667                         if (processedLabels.count(lastLabel) != 1) {
 668                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 669                                 needToRun = true;
 670                         }else {
 671                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 672                         }
 673                 }
 674
 675                 //run last label if you need to
 676                 if (needToRun == true)  {
 677                         for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } }
 678                         lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 679                 }
 680
 681                 return 0;
 682         }
 683         catch(exception& e) {
 684                 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
 685                 exit(1);
 686         }
 687 }
 688 //**********************************************************************************************************************
 689 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
 690         try {
 691
 692                 vector<SharedRAbundFloatVector*> newLookup;
 693                 for (int i = 0; i < thislookup.size(); i++) {
 694                         SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
 695                         temp->setLabel(thislookup[i]->getLabel());
 696                         temp->setGroup(thislookup[i]->getGroup());
 697                         newLookup.push_back(temp);
 698                 }
 699
 700                 //for each bin
 701                 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
 702                         if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
 703
 704                         //look at each sharedRabund and make sure they are not all zero
 705                         bool allZero = true;
 706                         for (int j = 0; j < thislookup.size(); j++) {
 707                                 if (thislookup[j]->getAbundance(i) != 0) { allZero = false;  break;  }
 708                         }
 709
 710                         //if they are not all zero add this bin
 711                         if (!allZero) {
 712                                 for (int j = 0; j < thislookup.size(); j++) {
 713                                         newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
 714                                 }
 715                         }
 716                 }
 717
 718                 for (int j = 0; j < thislookup.size(); j++) {  delete thislookup[j];  }
 719
 720                 thislookup = newLookup;
 721
 722                 return 0;
 723
 724         }
 725         catch(exception& e) {
 726                 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
 727                 exit(1);
 728         }
 729 }
 730 /*****************************************************************/
 731 map<string, vector<float> > CorrAxesCommand::readAxes(){
 732         try {
 733                 map<string, vector<float> > axes;
 734
 735                 ifstream in;
 736                 m->openInputFile(axesfile, in);
 737
 738                 string headerLine = m->getline(in); m->gobble(in);
 739
 740                 //count the number of axis you are reading
 741                 bool done = false;
 742                 int count = 0;
 743                 while (!done) {
 744                         int pos = headerLine.find("axis");
 745                         if (pos != string::npos) {
 746                                 count++;
 747                                 headerLine = headerLine.substr(pos+4);
 748                         }else { done = true; }
 749                 }
 750
 751                 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
 752
 753                 while (!in.eof()) {
 754
 755                         if (m->control_pressed) { in.close(); return axes; }
 756
 757                         string group = "";
 758                         in >> group; m->gobble(in);
 759
 760                         vector<float> thisGroupsAxes;
 761                         for (int i = 0; i < count; i++) {
 762                                 float temp = 0.0;
 763                                 in >> temp;
 764
 765                                 //only save the axis we want
 766                                 if (i < numaxes) {  thisGroupsAxes.push_back(temp); }
 767                         }
 768
 769                         //save group if its one the user selected
 770                         if (names.count(group) != 0) {
 771                                 map<string, vector<float> >::iterator it = axes.find(group);
 772
 773                                 if (it == axes.end()) {
 774                                         axes[group] = thisGroupsAxes;
 775                                 }else {
 776                                         m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
 777                                 }
 778                         }
 779
 780                         m->gobble(in);
 781                 }
 782                 in.close();
 783
 784                 return axes;
 785         }
 786         catch(exception& e) {
 787                 m->errorOut(e, "CorrAxesCommand", "readAxes");
 788                 exit(1);
 789         }
 790 }
 791 /*****************************************************************/
 792 int CorrAxesCommand::getMetadata(){
 793         try {
 794                 vector<string> groupNames;
 795
 796                 ifstream in;
 797                 m->openInputFile(axesfile, in);
 798
 799                 string headerLine = m->getline(in); m->gobble(in);
 800                 istringstream iss (headerLine,istringstream::in);
 801
 802                 //read the first label, because it refers to the groups
 803                 string columnLabel;
 804                 iss >> columnLabel; m->gobble(iss);
 805
 806                 //save names of columns you are reading
 807                 while (!iss.eof()) {
 808                         iss >> columnLabel; m->gobble(iss);
 809                         metadataLabels.push_back(columnLabel);
 810                 }
 811                 int count = metadataLabels.size();
 812
 813                 //read rest of file
 814                 while (!in.eof()) {
 815
 816                         if (m->control_pressed) { in.close(); return 0; }
 817
 818                         string group = "";
 819                         in >> group; m->gobble(in);
 820                         groupNames.push_back(group);
 821
 822                         SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
 823                         tempLookup->setGroup(group);
 824                         tempLookup->setLabel("1");
 825
 826                         for (int i = 0; i < count; i++) {
 827                                 float temp = 0.0;
 828                                 in >> temp;
 829
 830                                 tempLookup->push_back(temp, group);
 831                         }
 832
 833                         lookupFloat.push_back(tempLookup);
 834
 835                         m->gobble(in);
 836                 }
 837                 in.close();
 838
 839                 //remove any groups the user does not want, and set globaldata->groups with only valid groups
 840                 SharedUtil* util;
 841                 util = new SharedUtil();
 842
 843                 util->setGroups(globaldata->Groups, groupNames);
 844
 845                 for (int i = 0; i < lookupFloat.size(); i++) {
 846                         //if this sharedrabund is not from a group the user wants then delete it.
 847                         if (util->isValidGroup(lookupFloat[i]->getGroup(), globaldata->Groups) == false) {
 848                                 delete lookupFloat[i]; lookupFloat[i] = NULL;
 849                                 lookupFloat.erase(lookupFloat.begin()+i);
 850                                 i--;
 851                         }
 852                 }
 853
 854                 delete util;
 855
 856                 return 0;
 857         }
 858         catch(exception& e) {
 859                 m->errorOut(e, "CorrAxesCommand", "getMetadata");
 860                 exit(1);
 861         }
 862 }
 863 /*****************************************************************/
 864
 865
 866
 867
 868
 869