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