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(){
188 m->mothurOut("The corr.axes command reads a shared or relabund file as well as a pcoa file and calculates the correlation coefficient.\n");
189 m->mothurOut("The corr.axes command parameters are shared, relabund, axes, metadata, groups, method, numaxes and label. The shared or relabund and axes parameters are required. If shared is given the relative abundance is calculated.\n");
190 m->mothurOut("The metadata parameter.....\n");
191 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");
192 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");
193 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");
194 m->mothurOut("The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n");
195 m->mothurOut("The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n");
196 m->mothurOut("Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n");
197 m->mothurOut("The corr.axes command outputs a .corr.axes file.\n");
198 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
200 catch(exception& e) {
201 m->errorOut(e, "CorrAxesCommand", "help");
206 //**********************************************************************************************************************
208 CorrAxesCommand::~CorrAxesCommand(){}
210 //**********************************************************************************************************************
212 int CorrAxesCommand::execute(){
215 if (abort == true) { return 0; }
217 /*************************************************************************************/
218 // use smart distancing to get right sharedRabund and convert to relabund if needed //
219 /************************************************************************************/
220 if (sharedfile != "") {
222 if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; }
223 if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
225 //fills lookupFloat with relative abundance values from lookup
228 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
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; }
234 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
237 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
239 //this is for a sanity check to make sure the axes file and shared file match
240 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
242 /*************************************************************************************/
243 // read axes file and check for file mismatches with shared or relabund file //
244 /************************************************************************************/
247 map<string, vector<float> > axes = readAxes();
249 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
251 //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
252 if (axes.size() != lookupFloat.size()) {
253 map<string, vector<float> >::iterator it;
254 for (int i = 0; i < lookupFloat.size(); i++) {
255 it = axes.find(lookupFloat[i]->getGroup());
256 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(); }
258 m->control_pressed = true;
261 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
263 /*************************************************************************************/
264 // calc the r values //
265 /************************************************************************************/
267 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
268 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
270 m->openOutputFile(outputFileName, out);
271 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
275 for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; }
278 if (method == "pearson") { calcPearson(axes, out); }
279 else if (method == "spearman") { calcSpearman(axes, out); }
280 else if (method == "kendall") { calcKendall(axes, out); }
281 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
284 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
286 if (m->control_pressed) { return 0; }
288 m->mothurOutEndLine();
289 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
290 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
291 m->mothurOutEndLine();
295 catch(exception& e) {
296 m->errorOut(e, "CorrAxesCommand", "execute");
300 //**********************************************************************************************************************
301 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
304 //find average of each axis - X
305 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
306 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
307 vector<float> temp = it->second;
308 for (int i = 0; i < temp.size(); i++) {
309 averageAxes[i] += temp[i];
313 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
316 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
320 //find the averages this otu - Y
322 for (int j = 0; j < lookupFloat.size(); j++) {
323 sumOtu += lookupFloat[j]->getAbundance(i);
325 float Ybar = sumOtu / (float) lookupFloat.size();
327 //find r value for each axis
328 for (int k = 0; k < averageAxes.size(); k++) {
331 double numerator = 0.0;
332 double denomTerm1 = 0.0;
333 double denomTerm2 = 0.0;
334 for (int j = 0; j < lookupFloat.size(); j++) {
335 float Yi = lookupFloat[j]->getAbundance(i);
336 float Xi = axes[lookupFloat[j]->getGroup()][k];
338 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
339 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
340 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
343 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
345 r = numerator / denom;
355 catch(exception& e) {
356 m->errorOut(e, "CorrAxesCommand", "calcPearson");
360 //**********************************************************************************************************************
361 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
365 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
366 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
367 vector<float> temp = it->second;
368 for (int i = 0; i < temp.size(); i++) {
369 spearmanRank member(it->first, temp[i]);
370 scores[i].push_back(member);
375 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
377 //find ranks of xi in each axis
378 map<string, vector<float> > rankAxes;
379 for (int i = 0; i < numaxes; i++) {
381 vector<spearmanRank> ties;
383 for (int j = 0; j < scores[i].size(); j++) {
385 ties.push_back(scores[i][j]);
387 if (j != scores.size()-1) { // you are not the last so you can look ahead
388 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
389 for (int k = 0; k < ties.size(); k++) {
390 float thisrank = rankTotal / (float) ties.size();
391 rankAxes[ties[k].name].push_back(thisrank);
396 }else { // you are the last one
397 for (int k = 0; k < ties.size(); k++) {
398 float thisrank = rankTotal / (float) ties.size();
399 rankAxes[ties[k].name].push_back(thisrank);
408 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
412 //find the ranks of this otu - Y
413 vector<spearmanRank> otuScores;
414 for (int j = 0; j < lookupFloat.size(); j++) {
415 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
416 otuScores.push_back(member);
419 sort(otuScores.begin(), otuScores.end(), compareSpearman);
421 map<string, float> rankOtus;
422 vector<spearmanRank> ties;
424 for (int j = 0; j < otuScores.size(); j++) {
426 ties.push_back(otuScores[j]);
428 if (j != scores.size()-1) { // you are not the last so you can look ahead
429 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
430 for (int k = 0; k < ties.size(); k++) {
431 float thisrank = rankTotal / (float) ties.size();
432 rankOtus[ties[k].name] = thisrank;
437 }else { // you are the last one
438 for (int k = 0; k < ties.size(); k++) {
439 float thisrank = rankTotal / (float) ties.size();
440 rankOtus[ties[k].name] = thisrank;
445 //calc spearman ranks for each axis for this otu
446 for (int j = 0; j < numaxes; j++) {
449 for (int k = 0; k < lookupFloat.size(); k++) {
451 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
452 float yi = rankOtus[lookupFloat[k]->getGroup()];
454 di += ((xi - yi) * (xi - yi));
457 int n = lookupFloat.size();
458 double p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1)));
469 catch(exception& e) {
470 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
474 //**********************************************************************************************************************
475 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
479 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
480 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
481 vector<float> temp = it->second;
482 for (int i = 0; i < temp.size(); i++) {
483 spearmanRank member(it->first, temp[i]);
484 scores[i].push_back(member);
489 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
491 //convert scores to ranks of xi in each axis
492 for (int i = 0; i < numaxes; i++) {
494 vector<spearmanRank*> ties;
496 for (int j = 0; j < scores[i].size(); j++) {
498 ties.push_back(&(scores[i][j]));
500 if (j != scores.size()-1) { // you are not the last so you can look ahead
501 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
502 for (int k = 0; k < ties.size(); k++) {
503 float thisrank = rankTotal / (float) ties.size();
504 (*ties[k]).score = thisrank;
509 }else { // you are the last one
510 for (int k = 0; k < ties.size(); k++) {
511 float thisrank = rankTotal / (float) ties.size();
512 (*ties[k]).score = thisrank;
519 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
523 //find the ranks of this otu - Y
524 vector<spearmanRank> otuScores;
525 for (int j = 0; j < lookupFloat.size(); j++) {
526 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
527 otuScores.push_back(member);
530 sort(otuScores.begin(), otuScores.end(), compareSpearman);
532 map<string, float> rankOtus;
533 vector<spearmanRank> ties;
535 for (int j = 0; j < otuScores.size(); j++) {
537 ties.push_back(otuScores[j]);
539 if (j != scores.size()-1) { // you are not the last so you can look ahead
540 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
541 for (int k = 0; k < ties.size(); k++) {
542 float thisrank = rankTotal / (float) ties.size();
543 rankOtus[ties[k].name] = thisrank;
548 }else { // you are the last one
549 for (int k = 0; k < ties.size(); k++) {
550 float thisrank = rankTotal / (float) ties.size();
551 rankOtus[ties[k].name] = thisrank;
556 //calc spearman ranks for each axis for this otu
557 for (int j = 0; j < numaxes; j++) {
560 //assemble otus ranks in same order as axis ranks
561 vector<spearmanRank> otus;
562 for (int l = 0; l < scores[j].size(); l++) {
563 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
564 otus.push_back(member);
567 for (int l = 0; l < scores[j].size(); l++) {
569 int numWithHigherRank = 0;
570 float thisrank = otus[l].score;
572 for (int u = l; u < scores[j].size(); u++) {
573 if (otus[u].score > thisrank) { numWithHigherRank++; }
576 P += numWithHigherRank;
579 int n = lookupFloat.size();
581 double p = ( (4 * P) / (float) (n * (n - 1)) ) - 1.0;
591 catch(exception& e) {
592 m->errorOut(e, "CorrAxesCommand", "calcKendall");
596 //**********************************************************************************************************************
597 int CorrAxesCommand::getShared(){
599 InputData* input = new InputData(sharedfile, "sharedfile");
600 lookup = input->getSharedRAbundVectors();
601 string lastLabel = lookup[0]->getLabel();
603 if (label == "") { label = lastLabel; delete input; return 0; }
605 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
606 set<string> labels; labels.insert(label);
607 set<string> processedLabels;
608 set<string> userLabels = labels;
610 //as long as you are not at the end of the file or done wih the lines you want
611 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
612 if (m->control_pressed) { delete input; return 0; }
614 if(labels.count(lookup[0]->getLabel()) == 1){
615 processedLabels.insert(lookup[0]->getLabel());
616 userLabels.erase(lookup[0]->getLabel());
620 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
621 string saveLabel = lookup[0]->getLabel();
623 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
624 lookup = input->getSharedRAbundVectors(lastLabel);
626 processedLabels.insert(lookup[0]->getLabel());
627 userLabels.erase(lookup[0]->getLabel());
629 //restore real lastlabel to save below
630 lookup[0]->setLabel(saveLabel);
634 lastLabel = lookup[0]->getLabel();
636 //get next line to process
637 //prevent memory leak
638 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
639 lookup = input->getSharedRAbundVectors();
643 if (m->control_pressed) { delete input; return 0; }
645 //output error messages about any remaining user labels
646 set<string>::iterator it;
647 bool needToRun = false;
648 for (it = userLabels.begin(); it != userLabels.end(); it++) {
649 m->mothurOut("Your file does not include the label " + *it);
650 if (processedLabels.count(lastLabel) != 1) {
651 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
654 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
658 //run last label if you need to
659 if (needToRun == true) {
660 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
661 lookup = input->getSharedRAbundVectors(lastLabel);
667 catch(exception& e) {
668 m->errorOut(e, "CorrAxesCommand", "getShared");
672 //**********************************************************************************************************************
673 int CorrAxesCommand::getSharedFloat(){
675 InputData* input = new InputData(relabundfile, "relabund");
676 lookupFloat = input->getSharedRAbundFloatVectors();
677 string lastLabel = lookupFloat[0]->getLabel();
679 if (label == "") { label = lastLabel; delete input; return 0; }
681 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
682 set<string> labels; labels.insert(label);
683 set<string> processedLabels;
684 set<string> userLabels = labels;
686 //as long as you are not at the end of the file or done wih the lines you want
687 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
689 if (m->control_pressed) { delete input; return 0; }
691 if(labels.count(lookupFloat[0]->getLabel()) == 1){
692 processedLabels.insert(lookupFloat[0]->getLabel());
693 userLabels.erase(lookupFloat[0]->getLabel());
697 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
698 string saveLabel = lookupFloat[0]->getLabel();
700 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
701 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
703 processedLabels.insert(lookupFloat[0]->getLabel());
704 userLabels.erase(lookupFloat[0]->getLabel());
706 //restore real lastlabel to save below
707 lookupFloat[0]->setLabel(saveLabel);
711 lastLabel = lookupFloat[0]->getLabel();
713 //get next line to process
714 //prevent memory leak
715 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
716 lookupFloat = input->getSharedRAbundFloatVectors();
720 if (m->control_pressed) { delete input; return 0; }
722 //output error messages about any remaining user labels
723 set<string>::iterator it;
724 bool needToRun = false;
725 for (it = userLabels.begin(); it != userLabels.end(); it++) {
726 m->mothurOut("Your file does not include the label " + *it);
727 if (processedLabels.count(lastLabel) != 1) {
728 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
731 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
735 //run last label if you need to
736 if (needToRun == true) {
737 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
738 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
744 catch(exception& e) {
745 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
749 /*****************************************************************/
750 int CorrAxesCommand::convertToRelabund(){
753 vector<SharedRAbundFloatVector*> newLookup;
754 for (int i = 0; i < lookup.size(); i++) {
755 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
756 temp->setLabel(lookup[i]->getLabel());
757 temp->setGroup(lookup[i]->getGroup());
758 newLookup.push_back(temp);
761 for (int i = 0; i < lookup.size(); i++) {
763 for (int j = 0; j < lookup[i]->getNumBins(); j++) {
765 if (m->control_pressed) { return 0; }
767 int abund = lookup[i]->getAbundance(j);
769 float relabund = abund / (float) lookup[i]->getNumSeqs();
771 newLookup[i]->push_back(relabund, lookup[i]->getGroup());
775 if (pickedGroups) { eliminateZeroOTUS(newLookup); }
777 lookupFloat = newLookup;
781 catch(exception& e) {
782 m->errorOut(e, "CorrAxesCommand", "convertToRelabund");
786 //**********************************************************************************************************************
787 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
790 vector<SharedRAbundFloatVector*> newLookup;
791 for (int i = 0; i < thislookup.size(); i++) {
792 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
793 temp->setLabel(thislookup[i]->getLabel());
794 temp->setGroup(thislookup[i]->getGroup());
795 newLookup.push_back(temp);
799 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
800 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
802 //look at each sharedRabund and make sure they are not all zero
804 for (int j = 0; j < thislookup.size(); j++) {
805 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
808 //if they are not all zero add this bin
810 for (int j = 0; j < thislookup.size(); j++) {
811 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
816 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
818 thislookup = newLookup;
823 catch(exception& e) {
824 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
828 /*****************************************************************/
829 map<string, vector<float> > CorrAxesCommand::readAxes(){
831 map<string, vector<float> > axes;
834 m->openInputFile(axesfile, in);
836 string headerLine = m->getline(in); m->gobble(in);
838 //count the number of axis you are reading
842 int pos = headerLine.find("axis");
843 if (pos != string::npos) {
845 headerLine = headerLine.substr(pos+4);
846 }else { done = true; }
849 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
853 if (m->control_pressed) { in.close(); return axes; }
856 in >> group; m->gobble(in);
858 vector<float> thisGroupsAxes;
859 for (int i = 0; i < count; i++) {
863 //only save the axis we want
864 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
867 //save group if its one the user selected
868 if (names.count(group) != 0) {
869 map<string, vector<float> >::iterator it = axes.find(group);
871 if (it == axes.end()) {
872 axes[group] = thisGroupsAxes;
874 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
884 catch(exception& e) {
885 m->errorOut(e, "CorrAxesCommand", "convertToRelabund");
889 /*****************************************************************/