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