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