corraxescommand.cpp

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