corraxescommand.cpp

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