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