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