5 * Created by westcott on 12/22/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "corraxescommand.h"
11 #include "sharedutilities.h"
13 //********************************************************************************************************************
14 //sorts highes to lowest
15 inline bool compareSpearman(spearmanRank left, spearmanRank right){
16 return (left.score > right.score);
18 //**********************************************************************************************************************
19 vector<string> CorrAxesCommand::getValidParameters(){
21 string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
22 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
26 m->errorOut(e, "CorrAxesCommand", "getValidParameters");
30 //**********************************************************************************************************************
31 vector<string> CorrAxesCommand::getRequiredParameters(){
33 string Array[] = {"axes"};
34 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
38 m->errorOut(e, "CorrAxesCommand", "getRequiredParameters");
42 //**********************************************************************************************************************
43 CorrAxesCommand::CorrAxesCommand(){
46 //initialize outputTypes
47 vector<string> tempOutNames;
48 outputTypes["corr.axes"] = tempOutNames;
51 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
56 //**********************************************************************************************************************
57 vector<string> CorrAxesCommand::getRequiredFiles(){
59 vector<string> myArray;
63 m->errorOut(e, "CorrAxesCommand", "getRequiredFiles");
67 //**********************************************************************************************************************
68 CorrAxesCommand::CorrAxesCommand(string option) {
71 globaldata = GlobalData::getInstance();
73 //allow user to run help
74 if(option == "help") { help(); abort = true; }
77 //valid paramters for this command
78 string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
79 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
81 OptionParser parser(option);
82 map<string, string> parameters = parser.getParameters();
84 ValidParameters validParameter;
85 map<string, string>::iterator it;
87 //check to make sure all parameters are valid for command
88 for (it = parameters.begin(); it != parameters.end(); it++) {
89 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
92 vector<string> tempOutNames;
93 outputTypes["corr.axes"] = tempOutNames;
95 //if the user changes the input directory command factory will send this info to us in the output parameter
96 string inputDir = validParameter.validFile(parameters, "inputdir", false);
97 if (inputDir == "not found"){ inputDir = ""; }
100 it = parameters.find("axes");
101 //user has given a template file
102 if(it != parameters.end()){
103 path = m->hasPath(it->second);
104 //if the user has not given a path then, add inputdir. else leave path alone.
105 if (path == "") { parameters["axes"] = inputDir + it->second; }
108 it = parameters.find("shared");
109 //user has given a template file
110 if(it != parameters.end()){
111 path = m->hasPath(it->second);
112 //if the user has not given a path then, add inputdir. else leave path alone.
113 if (path == "") { parameters["shared"] = inputDir + it->second; }
116 it = parameters.find("relabund");
117 //user has given a template file
118 if(it != parameters.end()){
119 path = m->hasPath(it->second);
120 //if the user has not given a path then, add inputdir. else leave path alone.
121 if (path == "") { parameters["relabund"] = inputDir + it->second; }
124 it = parameters.find("metadata");
125 //user has given a template file
126 if(it != parameters.end()){
127 path = m->hasPath(it->second);
128 //if the user has not given a path then, add inputdir. else leave path alone.
129 if (path == "") { parameters["metadata"] = inputDir + it->second; }
134 //check for required parameters
135 axesfile = validParameter.validFile(parameters, "axes", true);
136 if (axesfile == "not open") { abort = true; }
137 else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true; }
139 sharedfile = validParameter.validFile(parameters, "shared", true);
140 if (sharedfile == "not open") { abort = true; }
141 else if (sharedfile == "not found") { sharedfile = ""; }
142 else { inputFileName = sharedfile; }
144 relabundfile = validParameter.validFile(parameters, "relabund", true);
145 if (relabundfile == "not open") { abort = true; }
146 else if (relabundfile == "not found") { relabundfile = ""; }
147 else { inputFileName = relabundfile; }
149 metadatafile = validParameter.validFile(parameters, "metadata", true);
150 if (metadatafile == "not open") { abort = true; }
151 else if (metadatafile == "not found") { metadatafile = ""; }
152 else { inputFileName = metadatafile; }
154 groups = validParameter.validFile(parameters, "groups", false);
155 if (groups == "not found") { groups = ""; pickedGroups = false; }
158 m->splitAtDash(groups, Groups);
160 globaldata->Groups = Groups;
162 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
164 label = validParameter.validFile(parameters, "label", false);
165 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=""; }
167 if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) { m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true; }
169 if (metadatafile != "") {
170 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
172 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
175 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
176 convert(temp, numaxes);
178 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
180 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; }
184 catch(exception& e) {
185 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
189 //**********************************************************************************************************************
191 void CorrAxesCommand::help(){
193 m->mothurOut("The corr.axes command reads a shared, relabund or metadata file as well as an axes file and calculates the correlation coefficient.\n");
194 m->mothurOut("The corr.axes command parameters are shared, relabund, axes, metadata, groups, method, numaxes and label. The shared, relabund or metadata and axes parameters are required. If shared is given the relative abundance is calculated.\n");
195 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");
196 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");
197 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");
198 m->mothurOut("The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n");
199 m->mothurOut("The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n");
200 m->mothurOut("Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n");
201 m->mothurOut("The corr.axes command outputs a .corr.axes file.\n");
202 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
204 catch(exception& e) {
205 m->errorOut(e, "CorrAxesCommand", "help");
210 //**********************************************************************************************************************
212 CorrAxesCommand::~CorrAxesCommand(){}
214 //**********************************************************************************************************************
216 int CorrAxesCommand::execute(){
219 if (abort == true) { return 0; }
221 /*************************************************************************************/
222 // use smart distancing to get right sharedRabund and convert to relabund if needed //
223 /************************************************************************************/
224 if (sharedfile != "") {
225 InputData* input = new InputData(sharedfile, "sharedfile");
226 getSharedFloat(input);
229 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
230 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
232 }else if (relabundfile != "") {
233 InputData* input = new InputData(relabundfile, "relabund");
234 getSharedFloat(input);
237 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
238 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
240 }else if (metadatafile != "") {
241 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
242 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
243 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
245 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
247 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
249 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
251 //this is for a sanity check to make sure the axes file and shared file match
252 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
254 /*************************************************************************************/
255 // read axes file and check for file mismatches with shared or relabund file //
256 /************************************************************************************/
259 map<string, vector<float> > axes = readAxes();
261 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
263 //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
264 if (axes.size() != lookupFloat.size()) {
265 map<string, vector<float> >::iterator it;
266 for (int i = 0; i < lookupFloat.size(); i++) {
267 it = axes.find(lookupFloat[i]->getGroup());
268 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(); }
270 m->control_pressed = true;
273 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
275 /*************************************************************************************/
276 // calc the r values //
277 /************************************************************************************/
279 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
280 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
282 m->openOutputFile(outputFileName, out);
283 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
286 if (metadatafile == "") { out << "OTU\t"; }
287 else { out << "Feature\t"; }
289 for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; }
292 if (method == "pearson") { calcPearson(axes, out); }
293 else if (method == "spearman") { calcSpearman(axes, out); }
294 else if (method == "kendall") { calcKendall(axes, out); }
295 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
298 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
300 if (m->control_pressed) { return 0; }
302 m->mothurOutEndLine();
303 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
304 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
305 m->mothurOutEndLine();
309 catch(exception& e) {
310 m->errorOut(e, "CorrAxesCommand", "execute");
314 //**********************************************************************************************************************
315 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
318 //find average of each axis - X
319 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
320 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
321 vector<float> temp = it->second;
322 for (int i = 0; i < temp.size(); i++) {
323 averageAxes[i] += temp[i];
327 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
330 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
332 if (metadatafile == "") { out << i+1 << '\t'; }
333 else { out << metadataLabels[i] << '\t'; }
335 //find the averages this otu - Y
337 for (int j = 0; j < lookupFloat.size(); j++) {
338 sumOtu += lookupFloat[j]->getAbundance(i);
340 float Ybar = sumOtu / (float) lookupFloat.size();
342 //find r value for each axis
343 for (int k = 0; k < averageAxes.size(); k++) {
346 double numerator = 0.0;
347 double denomTerm1 = 0.0;
348 double denomTerm2 = 0.0;
349 for (int j = 0; j < lookupFloat.size(); j++) {
350 float Yi = lookupFloat[j]->getAbundance(i);
351 float Xi = axes[lookupFloat[j]->getGroup()][k];
353 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
354 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
355 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
358 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
360 r = numerator / denom;
370 catch(exception& e) {
371 m->errorOut(e, "CorrAxesCommand", "calcPearson");
375 //**********************************************************************************************************************
376 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
380 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
381 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
382 vector<float> temp = it->second;
383 for (int i = 0; i < temp.size(); i++) {
384 spearmanRank member(it->first, temp[i]);
385 scores[i].push_back(member);
390 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
392 //find ranks of xi in each axis
393 map<string, vector<float> > rankAxes;
394 for (int i = 0; i < numaxes; i++) {
396 vector<spearmanRank> ties;
398 for (int j = 0; j < scores[i].size(); j++) {
400 ties.push_back(scores[i][j]);
402 if (j != scores.size()-1) { // you are not the last so you can look ahead
403 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
404 for (int k = 0; k < ties.size(); k++) {
405 float thisrank = rankTotal / (float) ties.size();
406 rankAxes[ties[k].name].push_back(thisrank);
411 }else { // you are the last one
412 for (int k = 0; k < ties.size(); k++) {
413 float thisrank = rankTotal / (float) ties.size();
414 rankAxes[ties[k].name].push_back(thisrank);
423 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
425 if (metadatafile == "") { out << i+1 << '\t'; }
426 else { out << metadataLabels[i] << '\t'; }
428 //find the ranks of this otu - Y
429 vector<spearmanRank> otuScores;
430 for (int j = 0; j < lookupFloat.size(); j++) {
431 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
432 otuScores.push_back(member);
435 sort(otuScores.begin(), otuScores.end(), compareSpearman);
437 map<string, float> rankOtus;
438 vector<spearmanRank> ties;
440 for (int j = 0; j < otuScores.size(); j++) {
442 ties.push_back(otuScores[j]);
444 if (j != scores.size()-1) { // you are not the last so you can look ahead
445 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
446 for (int k = 0; k < ties.size(); k++) {
447 float thisrank = rankTotal / (float) ties.size();
448 rankOtus[ties[k].name] = thisrank;
453 }else { // you are the last one
454 for (int k = 0; k < ties.size(); k++) {
455 float thisrank = rankTotal / (float) ties.size();
456 rankOtus[ties[k].name] = thisrank;
461 //calc spearman ranks for each axis for this otu
462 for (int j = 0; j < numaxes; j++) {
465 for (int k = 0; k < lookupFloat.size(); k++) {
467 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
468 float yi = rankOtus[lookupFloat[k]->getGroup()];
470 di += ((xi - yi) * (xi - yi));
473 int n = lookupFloat.size();
474 double p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1)));
485 catch(exception& e) {
486 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
490 //**********************************************************************************************************************
491 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
495 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
496 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
497 vector<float> temp = it->second;
498 for (int i = 0; i < temp.size(); i++) {
499 spearmanRank member(it->first, temp[i]);
500 scores[i].push_back(member);
505 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
507 //convert scores to ranks of xi in each axis
508 for (int i = 0; i < numaxes; i++) {
510 vector<spearmanRank*> ties;
512 for (int j = 0; j < scores[i].size(); j++) {
514 ties.push_back(&(scores[i][j]));
516 if (j != scores.size()-1) { // you are not the last so you can look ahead
517 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
518 for (int k = 0; k < ties.size(); k++) {
519 float thisrank = rankTotal / (float) ties.size();
520 (*ties[k]).score = thisrank;
525 }else { // you are the last one
526 for (int k = 0; k < ties.size(); k++) {
527 float thisrank = rankTotal / (float) ties.size();
528 (*ties[k]).score = thisrank;
535 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
537 if (metadatafile == "") { out << i+1 << '\t'; }
538 else { out << metadataLabels[i] << '\t'; }
540 //find the ranks of this otu - Y
541 vector<spearmanRank> otuScores;
542 for (int j = 0; j < lookupFloat.size(); j++) {
543 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
544 otuScores.push_back(member);
547 sort(otuScores.begin(), otuScores.end(), compareSpearman);
549 map<string, float> rankOtus;
550 vector<spearmanRank> ties;
552 for (int j = 0; j < otuScores.size(); j++) {
554 ties.push_back(otuScores[j]);
556 if (j != scores.size()-1) { // you are not the last so you can look ahead
557 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
558 for (int k = 0; k < ties.size(); k++) {
559 float thisrank = rankTotal / (float) ties.size();
560 rankOtus[ties[k].name] = thisrank;
565 }else { // you are the last one
566 for (int k = 0; k < ties.size(); k++) {
567 float thisrank = rankTotal / (float) ties.size();
568 rankOtus[ties[k].name] = thisrank;
573 //calc spearman ranks for each axis for this otu
574 for (int j = 0; j < numaxes; j++) {
577 //assemble otus ranks in same order as axis ranks
578 vector<spearmanRank> otus;
579 for (int l = 0; l < scores[j].size(); l++) {
580 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
581 otus.push_back(member);
584 for (int l = 0; l < scores[j].size(); l++) {
586 int numWithHigherRank = 0;
587 float thisrank = otus[l].score;
589 for (int u = l; u < scores[j].size(); u++) {
590 if (otus[u].score > thisrank) { numWithHigherRank++; }
593 P += numWithHigherRank;
596 int n = lookupFloat.size();
598 double p = ( (4 * P) / (float) (n * (n - 1)) ) - 1.0;
608 catch(exception& e) {
609 m->errorOut(e, "CorrAxesCommand", "calcKendall");
613 //**********************************************************************************************************************
614 int CorrAxesCommand::getSharedFloat(InputData* input){
616 lookupFloat = input->getSharedRAbundFloatVectors();
617 string lastLabel = lookupFloat[0]->getLabel();
619 if (label == "") { label = lastLabel; return 0; }
621 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
622 set<string> labels; labels.insert(label);
623 set<string> processedLabels;
624 set<string> userLabels = labels;
626 //as long as you are not at the end of the file or done wih the lines you want
627 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
629 if (m->control_pressed) { return 0; }
631 if(labels.count(lookupFloat[0]->getLabel()) == 1){
632 processedLabels.insert(lookupFloat[0]->getLabel());
633 userLabels.erase(lookupFloat[0]->getLabel());
637 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
638 string saveLabel = lookupFloat[0]->getLabel();
640 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
641 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
643 processedLabels.insert(lookupFloat[0]->getLabel());
644 userLabels.erase(lookupFloat[0]->getLabel());
646 //restore real lastlabel to save below
647 lookupFloat[0]->setLabel(saveLabel);
651 lastLabel = lookupFloat[0]->getLabel();
653 //get next line to process
654 //prevent memory leak
655 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
656 lookupFloat = input->getSharedRAbundFloatVectors();
660 if (m->control_pressed) { return 0; }
662 //output error messages about any remaining user labels
663 set<string>::iterator it;
664 bool needToRun = false;
665 for (it = userLabels.begin(); it != userLabels.end(); it++) {
666 m->mothurOut("Your file does not include the label " + *it);
667 if (processedLabels.count(lastLabel) != 1) {
668 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
671 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
675 //run last label if you need to
676 if (needToRun == true) {
677 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
678 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
683 catch(exception& e) {
684 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
688 //**********************************************************************************************************************
689 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
692 vector<SharedRAbundFloatVector*> newLookup;
693 for (int i = 0; i < thislookup.size(); i++) {
694 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
695 temp->setLabel(thislookup[i]->getLabel());
696 temp->setGroup(thislookup[i]->getGroup());
697 newLookup.push_back(temp);
701 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
702 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
704 //look at each sharedRabund and make sure they are not all zero
706 for (int j = 0; j < thislookup.size(); j++) {
707 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
710 //if they are not all zero add this bin
712 for (int j = 0; j < thislookup.size(); j++) {
713 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
718 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
720 thislookup = newLookup;
725 catch(exception& e) {
726 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
730 /*****************************************************************/
731 map<string, vector<float> > CorrAxesCommand::readAxes(){
733 map<string, vector<float> > axes;
736 m->openInputFile(axesfile, in);
738 string headerLine = m->getline(in); m->gobble(in);
740 //count the number of axis you are reading
744 int pos = headerLine.find("axis");
745 if (pos != string::npos) {
747 headerLine = headerLine.substr(pos+4);
748 }else { done = true; }
751 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
755 if (m->control_pressed) { in.close(); return axes; }
758 in >> group; m->gobble(in);
760 vector<float> thisGroupsAxes;
761 for (int i = 0; i < count; i++) {
765 //only save the axis we want
766 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
769 //save group if its one the user selected
770 if (names.count(group) != 0) {
771 map<string, vector<float> >::iterator it = axes.find(group);
773 if (it == axes.end()) {
774 axes[group] = thisGroupsAxes;
776 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
786 catch(exception& e) {
787 m->errorOut(e, "CorrAxesCommand", "readAxes");
791 /*****************************************************************/
792 int CorrAxesCommand::getMetadata(){
794 vector<string> groupNames;
797 m->openInputFile(axesfile, in);
799 string headerLine = m->getline(in); m->gobble(in);
800 istringstream iss (headerLine,istringstream::in);
802 //read the first label, because it refers to the groups
804 iss >> columnLabel; m->gobble(iss);
806 //save names of columns you are reading
808 iss >> columnLabel; m->gobble(iss);
809 metadataLabels.push_back(columnLabel);
811 int count = metadataLabels.size();
816 if (m->control_pressed) { in.close(); return 0; }
819 in >> group; m->gobble(in);
820 groupNames.push_back(group);
822 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
823 tempLookup->setGroup(group);
824 tempLookup->setLabel("1");
826 for (int i = 0; i < count; i++) {
830 tempLookup->push_back(temp, group);
833 lookupFloat.push_back(tempLookup);
839 //remove any groups the user does not want, and set globaldata->groups with only valid groups
841 util = new SharedUtil();
843 util->setGroups(globaldata->Groups, groupNames);
845 for (int i = 0; i < lookupFloat.size(); i++) {
846 //if this sharedrabund is not from a group the user wants then delete it.
847 if (util->isValidGroup(lookupFloat[i]->getGroup(), globaldata->Groups) == false) {
848 delete lookupFloat[i]; lookupFloat[i] = NULL;
849 lookupFloat.erase(lookupFloat.begin()+i);
858 catch(exception& e) {
859 m->errorOut(e, "CorrAxesCommand", "getMetadata");
863 /*****************************************************************/