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+1; 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                                 double p = (numCoor - numDisCoor) / (float) count;
 659
 660                                 out << '\t' << p;
 661                                 pValues[j] = p;
 662
 663                         }
 664
 665                         double sum = 0;
 666                         for(int k=0;k<numaxes;k++){
 667                                 sum += pValues[k] * pValues[k];
 668                         }
 669                         out << '\t' << sqrt(sum) << endl;
 670                 }
 671
 672                 return 0;
 673         }
 674         catch(exception& e) {
 675                 m->errorOut(e, "CorrAxesCommand", "calcKendall");
 676                 exit(1);
 677         }
 678 }
 679 //**********************************************************************************************************************
 680 int CorrAxesCommand::getSharedFloat(InputData* input){
 681         try {
 682                 lookupFloat = input->getSharedRAbundFloatVectors();
 683                 string lastLabel = lookupFloat[0]->getLabel();
 684
 685                 if (label == "") { label = lastLabel;  return 0; }
 686
 687                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 688                 set<string> labels; labels.insert(label);
 689                 set<string> processedLabels;
 690                 set<string> userLabels = labels;
 691
 692                 //as long as you are not at the end of the file or done wih the lines you want
 693                 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
 694
 695                         if (m->control_pressed) {  return 0;  }
 696
 697                         if(labels.count(lookupFloat[0]->getLabel()) == 1){
 698                                 processedLabels.insert(lookupFloat[0]->getLabel());
 699                                 userLabels.erase(lookupFloat[0]->getLabel());
 700                                 break;
 701                         }
 702
 703                         if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 704                                 string saveLabel = lookupFloat[0]->getLabel();
 705
 706                                 for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 707                                 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 708
 709                                 processedLabels.insert(lookupFloat[0]->getLabel());
 710                                 userLabels.erase(lookupFloat[0]->getLabel());
 711
 712                                 //restore real lastlabel to save below
 713                                 lookupFloat[0]->setLabel(saveLabel);
 714                                 break;
 715                         }
 716
 717                         lastLabel = lookupFloat[0]->getLabel();
 718
 719                         //get next line to process
 720                         //prevent memory leak
 721                         for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 722                         lookupFloat = input->getSharedRAbundFloatVectors();
 723                 }
 724
 725
 726                 if (m->control_pressed) { return 0;  }
 727
 728                 //output error messages about any remaining user labels
 729                 set<string>::iterator it;
 730                 bool needToRun = false;
 731                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 732                         m->mothurOut("Your file does not include the label " + *it);
 733                         if (processedLabels.count(lastLabel) != 1) {
 734                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 735                                 needToRun = true;
 736                         }else {
 737                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 738                         }
 739                 }
 740
 741                 //run last label if you need to
 742                 if (needToRun == true)  {
 743                         for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } }
 744                         lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 745                 }
 746
 747                 return 0;
 748         }
 749         catch(exception& e) {
 750                 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
 751                 exit(1);
 752         }
 753 }
 754 //**********************************************************************************************************************
 755 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
 756         try {
 757
 758                 vector<SharedRAbundFloatVector*> newLookup;
 759                 for (int i = 0; i < thislookup.size(); i++) {
 760                         SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
 761                         temp->setLabel(thislookup[i]->getLabel());
 762                         temp->setGroup(thislookup[i]->getGroup());
 763                         newLookup.push_back(temp);
 764                 }
 765
 766                 //for each bin
 767                 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
 768                         if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
 769
 770                         //look at each sharedRabund and make sure they are not all zero
 771                         bool allZero = true;
 772                         for (int j = 0; j < thislookup.size(); j++) {
 773                                 if (thislookup[j]->getAbundance(i) != 0) { allZero = false;  break;  }
 774                         }
 775
 776                         //if they are not all zero add this bin
 777                         if (!allZero) {
 778                                 for (int j = 0; j < thislookup.size(); j++) {
 779                                         newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
 780                                 }
 781                         }
 782                 }
 783
 784                 for (int j = 0; j < thislookup.size(); j++) {  delete thislookup[j];  }
 785
 786                 thislookup = newLookup;
 787
 788                 return 0;
 789
 790         }
 791         catch(exception& e) {
 792                 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
 793                 exit(1);
 794         }
 795 }
 796 /*****************************************************************/
 797 map<string, vector<float> > CorrAxesCommand::readAxes(){
 798         try {
 799                 map<string, vector<float> > axes;
 800
 801                 ifstream in;
 802                 m->openInputFile(axesfile, in);
 803
 804                 string headerLine = m->getline(in); m->gobble(in);
 805
 806                 //count the number of axis you are reading
 807                 bool done = false;
 808                 int count = 0;
 809                 while (!done) {
 810                         int pos = headerLine.find("axis");
 811                         if (pos != string::npos) {
 812                                 count++;
 813                                 headerLine = headerLine.substr(pos+4);
 814                         }else { done = true; }
 815                 }
 816
 817                 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
 818
 819                 while (!in.eof()) {
 820
 821                         if (m->control_pressed) { in.close(); return axes; }
 822
 823                         string group = "";
 824                         in >> group; m->gobble(in);
 825
 826                         vector<float> thisGroupsAxes;
 827                         for (int i = 0; i < count; i++) {
 828                                 float temp = 0.0;
 829                                 in >> temp;
 830
 831                                 //only save the axis we want
 832                                 if (i < numaxes) {  thisGroupsAxes.push_back(temp); }
 833                         }
 834
 835                         //save group if its one the user selected
 836                         if (names.count(group) != 0) {
 837                                 map<string, vector<float> >::iterator it = axes.find(group);
 838
 839                                 if (it == axes.end()) {
 840                                         axes[group] = thisGroupsAxes;
 841                                 }else {
 842                                         m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
 843                                 }
 844                         }
 845
 846                         m->gobble(in);
 847                 }
 848                 in.close();
 849
 850                 return axes;
 851         }
 852         catch(exception& e) {
 853                 m->errorOut(e, "CorrAxesCommand", "readAxes");
 854                 exit(1);
 855         }
 856 }
 857 /*****************************************************************/
 858 int CorrAxesCommand::getMetadata(){
 859         try {
 860                 vector<string> groupNames;
 861
 862                 ifstream in;
 863                 m->openInputFile(axesfile, in);
 864
 865                 string headerLine = m->getline(in); m->gobble(in);
 866                 istringstream iss (headerLine,istringstream::in);
 867
 868                 //read the first label, because it refers to the groups
 869                 string columnLabel;
 870                 iss >> columnLabel; m->gobble(iss);
 871
 872                 //save names of columns you are reading
 873                 while (!iss.eof()) {
 874                         iss >> columnLabel; m->gobble(iss);
 875                         metadataLabels.push_back(columnLabel);
 876                 }
 877                 int count = metadataLabels.size();
 878
 879                 //read rest of file
 880                 while (!in.eof()) {
 881
 882                         if (m->control_pressed) { in.close(); return 0; }
 883
 884                         string group = "";
 885                         in >> group; m->gobble(in);
 886                         groupNames.push_back(group);
 887
 888                         SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
 889                         tempLookup->setGroup(group);
 890                         tempLookup->setLabel("1");
 891
 892                         for (int i = 0; i < count; i++) {
 893                                 float temp = 0.0;
 894                                 in >> temp;
 895
 896                                 tempLookup->push_back(temp, group);
 897                         }
 898
 899                         lookupFloat.push_back(tempLookup);
 900
 901                         m->gobble(in);
 902                 }
 903                 in.close();
 904
 905                 //remove any groups the user does not want, and set globaldata->groups with only valid groups
 906                 SharedUtil* util;
 907                 util = new SharedUtil();
 908
 909                 util->setGroups(globaldata->Groups, groupNames);
 910
 911                 for (int i = 0; i < lookupFloat.size(); i++) {
 912                         //if this sharedrabund is not from a group the user wants then delete it.
 913                         if (util->isValidGroup(lookupFloat[i]->getGroup(), globaldata->Groups) == false) {
 914                                 delete lookupFloat[i]; lookupFloat[i] = NULL;
 915                                 lookupFloat.erase(lookupFloat.begin()+i);
 916                                 i--;
 917                         }
 918                 }
 919
 920                 delete util;
 921
 922                 return 0;
 923         }
 924         catch(exception& e) {
 925                 m->errorOut(e, "CorrAxesCommand", "getMetadata");
 926                 exit(1);
 927         }
 928 }
 929 /*****************************************************************/
 930
 931
 932
 933
 934
 935