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, relabund or metadata file as well as an axes 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";        }
 287                 else {  out << "Feature";                                               }
 288
 289                 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1); }
 290                 out << "\tlength" << 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;       }
 333                    else {  out << metadataLabels[i];            }
 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                    vector<float> rValues(averageAxes.size());
 343
 344                    //find r value for each axis
 345                    for (int k = 0; k < averageAxes.size(); k++) {
 346
 347                            double r = 0.0;
 348                            double numerator = 0.0;
 349                            double denomTerm1 = 0.0;
 350                            double denomTerm2 = 0.0;
 351                            for (int j = 0; j < lookupFloat.size(); j++) {
 352                                    float Yi = lookupFloat[j]->getAbundance(i);
 353                                    float Xi = axes[lookupFloat[j]->getGroup()][k];
 354
 355                                    numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
 356                                    denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
 357                                    denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
 358                            }
 359
 360                            double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
 361
 362                            r = numerator / denom;
 363                            rValues[k] = r;
 364                            out << '\t' << r;
 365                    }
 366
 367                    double sum = 0;
 368                    for(int k=0;k<rValues.size();k++){
 369                            sum += rValues[k] * rValues[k];
 370                    }
 371                    out << '\t' << sqrt(sum) << endl;
 372            }
 373
 374            return 0;
 375    }
 376    catch(exception& e) {
 377            m->errorOut(e, "CorrAxesCommand", "calcPearson");
 378            exit(1);
 379    }
 380 }
 381 //**********************************************************************************************************************
 382 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
 383         try {
 384
 385                 //format data
 386                 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
 387                 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
 388                         vector<float> temp = it->second;
 389                         for (int i = 0; i < temp.size(); i++) {
 390                                 spearmanRank member(it->first, temp[i]);
 391                                 scores[i].push_back(member);
 392                         }
 393                 }
 394
 395                 //sort each axis
 396                 for (int i = 0; i < numaxes; i++) {  sort(scores[i].begin(), scores[i].end(), compareSpearman); }
 397
 398                 //find ranks of xi in each axis
 399                 map<string, vector<float> > rankAxes;
 400                 for (int i = 0; i < numaxes; i++) {
 401
 402                         vector<spearmanRank> ties;
 403                         int rankTotal = 0;
 404                         for (int j = 0; j < scores[i].size(); j++) {
 405                                 rankTotal += j;
 406                                 ties.push_back(scores[i][j]);
 407
 408                                 if (j != scores.size()-1) { // you are not the last so you can look ahead
 409                                         if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
 410                                                 for (int k = 0; k < ties.size(); k++) {
 411                                                         float thisrank = rankTotal / (float) ties.size();
 412                                                         rankAxes[ties[k].name].push_back(thisrank);
 413                                                 }
 414                                                 ties.clear();
 415                                                 rankTotal = 0;
 416                                         }
 417                                 }else { // you are the last one
 418                                         for (int k = 0; k < ties.size(); k++) {
 419                                                 float thisrank = rankTotal / (float) ties.size();
 420                                                 rankAxes[ties[k].name].push_back(thisrank);
 421                                         }
 422                                 }
 423                         }
 424                 }
 425
 426
 427
 428                 //for each otu
 429                 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
 430
 431                         if (metadatafile == "") {  out << i+1;  }
 432                         else {  out << metadataLabels[i];               }
 433
 434                         //find the ranks of this otu - Y
 435                         vector<spearmanRank> otuScores;
 436                         for (int j = 0; j < lookupFloat.size(); j++) {
 437                                 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
 438                                 otuScores.push_back(member);
 439                         }
 440
 441                         sort(otuScores.begin(), otuScores.end(), compareSpearman);
 442
 443                         map<string, float> rankOtus;
 444                         vector<spearmanRank> ties;
 445                         int rankTotal = 0;
 446                         for (int j = 0; j < otuScores.size(); j++) {
 447                                 rankTotal += j;
 448                                 ties.push_back(otuScores[j]);
 449
 450                                 if (j != scores.size()-1) { // you are not the last so you can look ahead
 451                                         if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
 452                                                 for (int k = 0; k < ties.size(); k++) {
 453                                                         float thisrank = rankTotal / (float) ties.size();
 454                                                         rankOtus[ties[k].name] = thisrank;
 455                                                 }
 456                                                 ties.clear();
 457                                                 rankTotal = 0;
 458                                         }
 459                                 }else { // you are the last one
 460                                         for (int k = 0; k < ties.size(); k++) {
 461                                                 float thisrank = rankTotal / (float) ties.size();
 462                                                 rankOtus[ties[k].name] = thisrank;
 463                                         }
 464                                 }
 465                         }
 466
 467                         vector<double> pValues(numaxes);
 468                         //calc spearman ranks for each axis for this otu
 469                         for (int j = 0; j < numaxes; j++) {
 470
 471                                 double di = 0.0;
 472                                 for (int k = 0; k < lookupFloat.size(); k++) {
 473
 474                                         float xi = rankAxes[lookupFloat[k]->getGroup()][j];
 475                                         float yi = rankOtus[lookupFloat[k]->getGroup()];
 476
 477                                         di += ((xi - yi) * (xi - yi));
 478                                 }
 479
 480                                 int n = lookupFloat.size();
 481                                 double p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1)));
 482
 483                                 out  << '\t' << p;
 484                                 pValues[j] = p;
 485
 486                         }
 487
 488                         double sum = 0;
 489                         for(int k=0;k<numaxes;k++){
 490                                 sum += pValues[k] * pValues[k];
 491                         }
 492                         out << '\t' << sqrt(sum) << endl;
 493                 }
 494
 495                 return 0;
 496         }
 497         catch(exception& e) {
 498                 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
 499                 exit(1);
 500         }
 501 }
 502 //**********************************************************************************************************************
 503 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
 504         try {
 505
 506                 //format data
 507                 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
 508                 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
 509                         vector<float> temp = it->second;
 510                         for (int i = 0; i < temp.size(); i++) {
 511                                 spearmanRank member(it->first, temp[i]);
 512                                 scores[i].push_back(member);
 513                         }
 514                 }
 515
 516                 //sort each axis
 517                 for (int i = 0; i < numaxes; i++) {  sort(scores[i].begin(), scores[i].end(), compareSpearman); }
 518
 519                 //convert scores to ranks of xi in each axis
 520                 for (int i = 0; i < numaxes; i++) {
 521
 522                         vector<spearmanRank*> ties;
 523                         int rankTotal = 0;
 524                         for (int j = 0; j < scores[i].size(); j++) {
 525                                 rankTotal += j;
 526                                 ties.push_back(&(scores[i][j]));
 527
 528                                 if (j != scores.size()-1) { // you are not the last so you can look ahead
 529                                         if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
 530                                                 for (int k = 0; k < ties.size(); k++) {
 531                                                         float thisrank = rankTotal / (float) ties.size();
 532                                                         (*ties[k]).score = thisrank;
 533                                                 }
 534                                                 ties.clear();
 535                                                 rankTotal = 0;
 536                                         }
 537                                 }else { // you are the last one
 538                                         for (int k = 0; k < ties.size(); k++) {
 539                                                 float thisrank = rankTotal / (float) ties.size();
 540                                                 (*ties[k]).score = thisrank;
 541                                         }
 542                                 }
 543                         }
 544                 }
 545
 546                 //for each otu
 547                 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
 548
 549                         if (metadatafile == "") {  out << i+1;  }
 550                         else {  out << metadataLabels[i];               }
 551
 552                         //find the ranks of this otu - Y
 553                         vector<spearmanRank> otuScores;
 554                         for (int j = 0; j < lookupFloat.size(); j++) {
 555                                 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
 556                                 otuScores.push_back(member);
 557                         }
 558
 559                         sort(otuScores.begin(), otuScores.end(), compareSpearman);
 560
 561                         map<string, float> rankOtus;
 562                         vector<spearmanRank> ties;
 563                         int rankTotal = 0;
 564                         for (int j = 0; j < otuScores.size(); j++) {
 565                                 rankTotal += j;
 566                                 ties.push_back(otuScores[j]);
 567
 568                                 if (j != scores.size()-1) { // you are not the last so you can look ahead
 569                                         if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
 570                                                 for (int k = 0; k < ties.size(); k++) {
 571                                                         float thisrank = rankTotal / (float) ties.size();
 572                                                         rankOtus[ties[k].name] = thisrank;
 573                                                 }
 574                                                 ties.clear();
 575                                                 rankTotal = 0;
 576                                         }
 577                                 }else { // you are the last one
 578                                         for (int k = 0; k < ties.size(); k++) {
 579                                                 float thisrank = rankTotal / (float) ties.size();
 580                                                 rankOtus[ties[k].name] = thisrank;
 581                                         }
 582                                 }
 583                         }
 584                         vector<double> pValues(numaxes);
 585
 586                         //calc spearman ranks for each axis for this otu
 587                         for (int j = 0; j < numaxes; j++) {
 588
 589                                 int P = 0;
 590                                 //assemble otus ranks in same order as axis ranks
 591                                 vector<spearmanRank> otus;
 592                                 for (int l = 0; l < scores[j].size(); l++) {
 593                                         spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
 594                                         otus.push_back(member);
 595                                 }
 596
 597                                 for (int l = 0; l < scores[j].size(); l++) {
 598
 599                                         int numWithHigherRank = 0;
 600                                         float thisrank = otus[l].score;
 601
 602                                         for (int u = l; u < scores[j].size(); u++) {
 603                                                 if (otus[u].score > thisrank) { numWithHigherRank++; }
 604                                         }
 605
 606                                         P += numWithHigherRank;
 607                                 }
 608
 609                                 int n = lookupFloat.size();
 610
 611                                 double p = ( (4 * P) / (float) (n * (n - 1)) ) - 1.0;
 612
 613                                 out << '\t' << p;
 614                                 pValues[j] = p;
 615
 616                         }
 617
 618                         double sum = 0;
 619                         for(int k=0;k<numaxes;k++){
 620                                 sum += pValues[k] * pValues[k];
 621                         }
 622                         out << '\t' << sqrt(sum) << endl;
 623                 }
 624
 625                 return 0;
 626         }
 627         catch(exception& e) {
 628                 m->errorOut(e, "CorrAxesCommand", "calcKendall");
 629                 exit(1);
 630         }
 631 }
 632 //**********************************************************************************************************************
 633 int CorrAxesCommand::getSharedFloat(InputData* input){
 634         try {
 635                 lookupFloat = input->getSharedRAbundFloatVectors();
 636                 string lastLabel = lookupFloat[0]->getLabel();
 637
 638                 if (label == "") { label = lastLabel;  return 0; }
 639
 640                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 641                 set<string> labels; labels.insert(label);
 642                 set<string> processedLabels;
 643                 set<string> userLabels = labels;
 644
 645                 //as long as you are not at the end of the file or done wih the lines you want
 646                 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
 647
 648                         if (m->control_pressed) {  return 0;  }
 649
 650                         if(labels.count(lookupFloat[0]->getLabel()) == 1){
 651                                 processedLabels.insert(lookupFloat[0]->getLabel());
 652                                 userLabels.erase(lookupFloat[0]->getLabel());
 653                                 break;
 654                         }
 655
 656                         if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 657                                 string saveLabel = lookupFloat[0]->getLabel();
 658
 659                                 for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 660                                 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 661
 662                                 processedLabels.insert(lookupFloat[0]->getLabel());
 663                                 userLabels.erase(lookupFloat[0]->getLabel());
 664
 665                                 //restore real lastlabel to save below
 666                                 lookupFloat[0]->setLabel(saveLabel);
 667                                 break;
 668                         }
 669
 670                         lastLabel = lookupFloat[0]->getLabel();
 671
 672                         //get next line to process
 673                         //prevent memory leak
 674                         for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 675                         lookupFloat = input->getSharedRAbundFloatVectors();
 676                 }
 677
 678
 679                 if (m->control_pressed) { return 0;  }
 680
 681                 //output error messages about any remaining user labels
 682                 set<string>::iterator it;
 683                 bool needToRun = false;
 684                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 685                         m->mothurOut("Your file does not include the label " + *it);
 686                         if (processedLabels.count(lastLabel) != 1) {
 687                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 688                                 needToRun = true;
 689                         }else {
 690                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 691                         }
 692                 }
 693
 694                 //run last label if you need to
 695                 if (needToRun == true)  {
 696                         for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } }
 697                         lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 698                 }
 699
 700                 return 0;
 701         }
 702         catch(exception& e) {
 703                 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
 704                 exit(1);
 705         }
 706 }
 707 //**********************************************************************************************************************
 708 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
 709         try {
 710
 711                 vector<SharedRAbundFloatVector*> newLookup;
 712                 for (int i = 0; i < thislookup.size(); i++) {
 713                         SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
 714                         temp->setLabel(thislookup[i]->getLabel());
 715                         temp->setGroup(thislookup[i]->getGroup());
 716                         newLookup.push_back(temp);
 717                 }
 718
 719                 //for each bin
 720                 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
 721                         if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
 722
 723                         //look at each sharedRabund and make sure they are not all zero
 724                         bool allZero = true;
 725                         for (int j = 0; j < thislookup.size(); j++) {
 726                                 if (thislookup[j]->getAbundance(i) != 0) { allZero = false;  break;  }
 727                         }
 728
 729                         //if they are not all zero add this bin
 730                         if (!allZero) {
 731                                 for (int j = 0; j < thislookup.size(); j++) {
 732                                         newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
 733                                 }
 734                         }
 735                 }
 736
 737                 for (int j = 0; j < thislookup.size(); j++) {  delete thislookup[j];  }
 738
 739                 thislookup = newLookup;
 740
 741                 return 0;
 742
 743         }
 744         catch(exception& e) {
 745                 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
 746                 exit(1);
 747         }
 748 }
 749 /*****************************************************************/
 750 map<string, vector<float> > CorrAxesCommand::readAxes(){
 751         try {
 752                 map<string, vector<float> > axes;
 753
 754                 ifstream in;
 755                 m->openInputFile(axesfile, in);
 756
 757                 string headerLine = m->getline(in); m->gobble(in);
 758
 759                 //count the number of axis you are reading
 760                 bool done = false;
 761                 int count = 0;
 762                 while (!done) {
 763                         int pos = headerLine.find("axis");
 764                         if (pos != string::npos) {
 765                                 count++;
 766                                 headerLine = headerLine.substr(pos+4);
 767                         }else { done = true; }
 768                 }
 769
 770                 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
 771
 772                 while (!in.eof()) {
 773
 774                         if (m->control_pressed) { in.close(); return axes; }
 775
 776                         string group = "";
 777                         in >> group; m->gobble(in);
 778
 779                         vector<float> thisGroupsAxes;
 780                         for (int i = 0; i < count; i++) {
 781                                 float temp = 0.0;
 782                                 in >> temp;
 783
 784                                 //only save the axis we want
 785                                 if (i < numaxes) {  thisGroupsAxes.push_back(temp); }
 786                         }
 787
 788                         //save group if its one the user selected
 789                         if (names.count(group) != 0) {
 790                                 map<string, vector<float> >::iterator it = axes.find(group);
 791
 792                                 if (it == axes.end()) {
 793                                         axes[group] = thisGroupsAxes;
 794                                 }else {
 795                                         m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
 796                                 }
 797                         }
 798
 799                         m->gobble(in);
 800                 }
 801                 in.close();
 802
 803                 return axes;
 804         }
 805         catch(exception& e) {
 806                 m->errorOut(e, "CorrAxesCommand", "readAxes");
 807                 exit(1);
 808         }
 809 }
 810 /*****************************************************************/
 811 int CorrAxesCommand::getMetadata(){
 812         try {
 813                 vector<string> groupNames;
 814
 815                 ifstream in;
 816                 m->openInputFile(axesfile, in);
 817
 818                 string headerLine = m->getline(in); m->gobble(in);
 819                 istringstream iss (headerLine,istringstream::in);
 820
 821                 //read the first label, because it refers to the groups
 822                 string columnLabel;
 823                 iss >> columnLabel; m->gobble(iss);
 824
 825                 //save names of columns you are reading
 826                 while (!iss.eof()) {
 827                         iss >> columnLabel; m->gobble(iss);
 828                         metadataLabels.push_back(columnLabel);
 829                 }
 830                 int count = metadataLabels.size();
 831
 832                 //read rest of file
 833                 while (!in.eof()) {
 834
 835                         if (m->control_pressed) { in.close(); return 0; }
 836
 837                         string group = "";
 838                         in >> group; m->gobble(in);
 839                         groupNames.push_back(group);
 840
 841                         SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
 842                         tempLookup->setGroup(group);
 843                         tempLookup->setLabel("1");
 844
 845                         for (int i = 0; i < count; i++) {
 846                                 float temp = 0.0;
 847                                 in >> temp;
 848
 849                                 tempLookup->push_back(temp, group);
 850                         }
 851
 852                         lookupFloat.push_back(tempLookup);
 853
 854                         m->gobble(in);
 855                 }
 856                 in.close();
 857
 858                 //remove any groups the user does not want, and set globaldata->groups with only valid groups
 859                 SharedUtil* util;
 860                 util = new SharedUtil();
 861
 862                 util->setGroups(globaldata->Groups, groupNames);
 863
 864                 for (int i = 0; i < lookupFloat.size(); i++) {
 865                         //if this sharedrabund is not from a group the user wants then delete it.
 866                         if (util->isValidGroup(lookupFloat[i]->getGroup(), globaldata->Groups) == false) {
 867                                 delete lookupFloat[i]; lookupFloat[i] = NULL;
 868                                 lookupFloat.erase(lookupFloat.begin()+i);
 869                                 i--;
 870                         }
 871                 }
 872
 873                 delete util;
 874
 875                 return 0;
 876         }
 877         catch(exception& e) {
 878                 m->errorOut(e, "CorrAxesCommand", "getMetadata");
 879                 exit(1);
 880         }
 881 }
 882 /*****************************************************************/
 883
 884
 885
 886
 887
 888