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