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