corraxescommand.cpp

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