5 * Created by westcott on 12/22/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "corraxescommand.h"
12 //********************************************************************************************************************
13 //sorts highes to lowest
14 inline bool compareSpearman(spearmanRank left, spearmanRank right){
15 return (left.score > right.score);
17 //**********************************************************************************************************************
18 vector<string> CorrAxesCommand::getValidParameters(){
20 string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
21 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
25 m->errorOut(e, "CorrAxesCommand", "getValidParameters");
29 //**********************************************************************************************************************
30 vector<string> CorrAxesCommand::getRequiredParameters(){
32 string Array[] = {"axes"};
33 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
37 m->errorOut(e, "CorrAxesCommand", "getRequiredParameters");
41 //**********************************************************************************************************************
42 CorrAxesCommand::CorrAxesCommand(){
45 //initialize outputTypes
46 vector<string> tempOutNames;
47 outputTypes["corr.axes"] = tempOutNames;
50 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
55 //**********************************************************************************************************************
56 vector<string> CorrAxesCommand::getRequiredFiles(){
58 vector<string> myArray;
62 m->errorOut(e, "CorrAxesCommand", "getRequiredFiles");
66 //**********************************************************************************************************************
67 CorrAxesCommand::CorrAxesCommand(string option) {
70 globaldata = GlobalData::getInstance();
72 //allow user to run help
73 if(option == "help") { help(); abort = true; }
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)));
80 OptionParser parser(option);
81 map<string, string> parameters = parser.getParameters();
83 ValidParameters validParameter;
84 map<string, string>::iterator it;
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; }
91 vector<string> tempOutNames;
92 outputTypes["corr.axes"] = tempOutNames;
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 = ""; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
148 metadatafile = validParameter.validFile(parameters, "metadata", true);
149 if (metadatafile == "not open") { abort = true; }
150 else if (metadatafile == "not found") { metadatafile = ""; }
152 groups = validParameter.validFile(parameters, "groups", false);
153 if (groups == "not found") { groups = ""; pickedGroups = false; }
156 m->splitAtDash(groups, Groups);
158 globaldata->Groups = Groups;
160 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
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=""; }
165 if ((relabundfile == "") && (sharedfile == "")) { m->mothurOut("You must provide either a shared or relabund file."); m->mothurOutEndLine(); abort = true; }
167 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true; }
170 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
171 convert(temp, numaxes);
173 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
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; }
179 catch(exception& e) {
180 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
184 //**********************************************************************************************************************
186 void CorrAxesCommand::help(){
190 catch(exception& e) {
191 m->errorOut(e, "CorrAxesCommand", "help");
196 //**********************************************************************************************************************
198 CorrAxesCommand::~CorrAxesCommand(){}
200 //**********************************************************************************************************************
202 int CorrAxesCommand::execute(){
205 if (abort == true) { return 0; }
207 /*************************************************************************************/
208 // use smart distancing to get right sharedRabund and convert to relabund if needed //
209 /************************************************************************************/
210 if (sharedfile != "") {
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; }
215 //fills lookupFloat with relative abundance values from lookup
218 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
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; }
224 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
227 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
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()); }
232 /*************************************************************************************/
233 // read axes file and check for file mismatches with shared or relabund file //
234 /************************************************************************************/
237 map<string, vector<float> > axes = readAxes();
239 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
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(); }
248 m->control_pressed = true;
251 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
253 /*************************************************************************************/
254 // calc the r values //
255 /************************************************************************************/
257 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + "corr.axes";
258 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
260 m->openOutputFile(outputFileName, out);
261 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
265 for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; }
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(); }
274 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
276 if (m->control_pressed) { return 0; }
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();
285 catch(exception& e) {
286 m->errorOut(e, "CorrAxesCommand", "execute");
290 //**********************************************************************************************************************
291 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
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];
303 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
306 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
310 //find the averages this otu - Y
312 for (int j = 0; j < lookupFloat.size(); j++) {
313 sumOtu += lookupFloat[j]->getAbundance(i);
315 float Ybar = sumOtu / (float) lookupFloat.size();
317 //find r value for each axis
318 for (int k = 0; k < averageAxes.size(); k++) {
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];
328 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
329 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
330 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
333 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
335 r = numerator / denom;
345 catch(exception& e) {
346 m->errorOut(e, "CorrAxesCommand", "calcPearson");
350 //**********************************************************************************************************************
351 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
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);
365 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
367 //find ranks of xi in each axis
368 map<string, vector<float> > rankAxes;
369 for (int i = 0; i < numaxes; i++) {
371 vector<spearmanRank> ties;
373 for (int j = 0; j < scores[i].size(); j++) {
375 ties.push_back(scores[i][j]);
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);
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);
398 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
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);
409 sort(otuScores.begin(), otuScores.end(), compareSpearman);
411 map<string, float> rankOtus;
412 vector<spearmanRank> ties;
414 for (int j = 0; j < otuScores.size(); j++) {
416 ties.push_back(otuScores[j]);
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;
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;
435 //calc spearman ranks for each axis for this otu
436 for (int j = 0; j < numaxes; j++) {
439 for (int k = 0; k < lookupFloat.size(); k++) {
441 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
442 float yi = rankOtus[lookupFloat[k]->getGroup()];
444 di += ((xi - yi) * (xi - yi));
447 int n = lookupFloat.size();
448 double p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1)));
459 catch(exception& e) {
460 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
464 //**********************************************************************************************************************
465 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
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);
479 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
481 //find ranks of xi in each axis
482 map<string, vector<float> > rankAxes;
483 for (int i = 0; i < numaxes; i++) {
485 vector<spearmanRank> ties;
487 for (int j = 0; j < scores[i].size(); j++) {
489 ties.push_back(scores[i][j]);
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);
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);
512 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
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);
523 sort(otuScores.begin(), otuScores.end(), compareSpearman);
525 map<string, float> rankOtus;
526 vector<spearmanRank> ties;
528 for (int j = 0; j < otuScores.size(); j++) {
530 ties.push_back(otuScores[j]);
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;
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;
549 //calc kendall coeffient for each axis for this otu
550 for (int j = 0; j < numaxes; j++) {
552 int numConcordant = 0;
553 int numDiscordant = 0;
555 for (int f = 0; f < j; f++) {
557 for (int k = 0; k < lookupFloat.size(); k++) {
559 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
560 float yi = rankOtus[lookupFloat[k]->getGroup()];
562 for (int h = 0; h < k; h++) {
564 float xj = rankAxes[lookupFloat[k]->getGroup()][f];
565 float yj = rankOtus[lookupFloat[h]->getGroup()];
567 if ( ((xi > xj) && (yi < yj)) || ((xi < xj) && (yi > yj)) ){ numDiscordant++; }
568 if ( ((xi > xj) && (yi > yj)) || ((xi < xj) && (yi < yj)) ){ numConcordant++; }
573 int n = lookupFloat.size();
574 double p = (numConcordant - numDiscordant) / (float) (0.5 * n * (n - 1));
585 catch(exception& e) {
586 m->errorOut(e, "CorrAxesCommand", "calcKendall");
590 //**********************************************************************************************************************
591 int CorrAxesCommand::getShared(){
593 InputData* input = new InputData(sharedfile, "sharedfile");
594 lookup = input->getSharedRAbundVectors();
595 string lastLabel = lookup[0]->getLabel();
597 if (label == "") { label = lastLabel; delete input; return 0; }
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;
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; }
608 if(labels.count(lookup[0]->getLabel()) == 1){
609 processedLabels.insert(lookup[0]->getLabel());
610 userLabels.erase(lookup[0]->getLabel());
614 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
615 string saveLabel = lookup[0]->getLabel();
617 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
618 lookup = input->getSharedRAbundVectors(lastLabel);
620 processedLabels.insert(lookup[0]->getLabel());
621 userLabels.erase(lookup[0]->getLabel());
623 //restore real lastlabel to save below
624 lookup[0]->setLabel(saveLabel);
628 lastLabel = lookup[0]->getLabel();
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();
637 if (m->control_pressed) { delete input; return 0; }
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();
648 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
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);
661 catch(exception& e) {
662 m->errorOut(e, "CorrAxesCommand", "getShared");
666 //**********************************************************************************************************************
667 int CorrAxesCommand::getSharedFloat(){
669 InputData* input = new InputData(relabundfile, "relabund");
670 lookupFloat = input->getSharedRAbundFloatVectors();
671 string lastLabel = lookupFloat[0]->getLabel();
673 if (label == "") { label = lastLabel; delete input; return 0; }
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;
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)) {
683 if (m->control_pressed) { delete input; return 0; }
685 if(labels.count(lookupFloat[0]->getLabel()) == 1){
686 processedLabels.insert(lookupFloat[0]->getLabel());
687 userLabels.erase(lookupFloat[0]->getLabel());
691 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
692 string saveLabel = lookupFloat[0]->getLabel();
694 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
695 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
697 processedLabels.insert(lookupFloat[0]->getLabel());
698 userLabels.erase(lookupFloat[0]->getLabel());
700 //restore real lastlabel to save below
701 lookupFloat[0]->setLabel(saveLabel);
705 lastLabel = lookupFloat[0]->getLabel();
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();
714 if (m->control_pressed) { delete input; return 0; }
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();
725 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
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);
738 catch(exception& e) {
739 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
743 /*****************************************************************/
744 int CorrAxesCommand::convertToRelabund(){
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);
755 for (int i = 0; i < lookup.size(); i++) {
757 for (int j = 0; j < lookup[i]->getNumBins(); j++) {
759 if (m->control_pressed) { return 0; }
761 int abund = lookup[i]->getAbundance(j);
763 float relabund = abund / (float) lookup[i]->getNumSeqs();
765 newLookup[i]->push_back(relabund, lookup[i]->getGroup());
769 if (pickedGroups) { eliminateZeroOTUS(newLookup); }
771 lookupFloat = newLookup;
775 catch(exception& e) {
776 m->errorOut(e, "CorrAxesCommand", "convertToRelabund");
780 //**********************************************************************************************************************
781 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
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);
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; }
796 //look at each sharedRabund and make sure they are not all zero
798 for (int j = 0; j < thislookup.size(); j++) {
799 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
802 //if they are not all zero add this bin
804 for (int j = 0; j < thislookup.size(); j++) {
805 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
810 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
812 thislookup = newLookup;
817 catch(exception& e) {
818 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
822 /*****************************************************************/
823 map<string, vector<float> > CorrAxesCommand::readAxes(){
825 map<string, vector<float> > axes;
828 m->openInputFile(axesfile, in);
830 string headerLine = m->getline(in); m->gobble(in);
832 //count the number of axis you are reading
836 int pos = headerLine.find("axis");
837 if (pos != string::npos) {
839 headerLine = headerLine.substr(pos+4);
840 }else { done = true; }
845 if (m->control_pressed) { in.close(); return axes; }
848 in >> group; m->gobble(in);
850 vector<float> thisGroupsAxes;
851 for (int i = 0; i < count; i++) {
855 //only save the axis we want
856 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
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);
863 if (it == axes.end()) {
864 axes[group] = thisGroupsAxes;
866 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
876 catch(exception& e) {
877 m->errorOut(e, "CorrAxesCommand", "convertToRelabund");
881 /*****************************************************************/