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 highest to lowest
15 inline bool compareSpearman(spearmanRank left, spearmanRank right){
16 return (left.score > right.score);
18 //********************************************************************************************************************
19 //sorts lowest to highest
20 inline bool compareSpearmanReverse(spearmanRank left, spearmanRank right){
21 return (left.score < right.score);
23 //**********************************************************************************************************************
24 vector<string> CorrAxesCommand::getValidParameters(){
26 string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
27 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
31 m->errorOut(e, "CorrAxesCommand", "getValidParameters");
35 //**********************************************************************************************************************
36 vector<string> CorrAxesCommand::getRequiredParameters(){
38 string Array[] = {"axes"};
39 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
43 m->errorOut(e, "CorrAxesCommand", "getRequiredParameters");
47 //**********************************************************************************************************************
48 CorrAxesCommand::CorrAxesCommand(){
51 //initialize outputTypes
52 vector<string> tempOutNames;
53 outputTypes["corr.axes"] = tempOutNames;
56 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
61 //**********************************************************************************************************************
62 vector<string> CorrAxesCommand::getRequiredFiles(){
64 vector<string> myArray;
68 m->errorOut(e, "CorrAxesCommand", "getRequiredFiles");
72 //**********************************************************************************************************************
73 CorrAxesCommand::CorrAxesCommand(string option) {
76 globaldata = GlobalData::getInstance();
78 //allow user to run help
79 if(option == "help") { help(); abort = true; }
82 //valid paramters for this command
83 string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
84 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
86 OptionParser parser(option);
87 map<string, string> parameters = parser.getParameters();
89 ValidParameters validParameter;
90 map<string, string>::iterator it;
92 //check to make sure all parameters are valid for command
93 for (it = parameters.begin(); it != parameters.end(); it++) {
94 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
97 vector<string> tempOutNames;
98 outputTypes["corr.axes"] = tempOutNames;
100 //if the user changes the input directory command factory will send this info to us in the output parameter
101 string inputDir = validParameter.validFile(parameters, "inputdir", false);
102 if (inputDir == "not found"){ inputDir = ""; }
105 it = parameters.find("axes");
106 //user has given a template file
107 if(it != parameters.end()){
108 path = m->hasPath(it->second);
109 //if the user has not given a path then, add inputdir. else leave path alone.
110 if (path == "") { parameters["axes"] = inputDir + it->second; }
113 it = parameters.find("shared");
114 //user has given a template file
115 if(it != parameters.end()){
116 path = m->hasPath(it->second);
117 //if the user has not given a path then, add inputdir. else leave path alone.
118 if (path == "") { parameters["shared"] = inputDir + it->second; }
121 it = parameters.find("relabund");
122 //user has given a template file
123 if(it != parameters.end()){
124 path = m->hasPath(it->second);
125 //if the user has not given a path then, add inputdir. else leave path alone.
126 if (path == "") { parameters["relabund"] = inputDir + it->second; }
129 it = parameters.find("metadata");
130 //user has given a template file
131 if(it != parameters.end()){
132 path = m->hasPath(it->second);
133 //if the user has not given a path then, add inputdir. else leave path alone.
134 if (path == "") { parameters["metadata"] = inputDir + it->second; }
139 //check for required parameters
140 axesfile = validParameter.validFile(parameters, "axes", true);
141 if (axesfile == "not open") { abort = true; }
142 else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true; }
144 sharedfile = validParameter.validFile(parameters, "shared", true);
145 if (sharedfile == "not open") { abort = true; }
146 else if (sharedfile == "not found") { sharedfile = ""; }
147 else { inputFileName = sharedfile; }
149 relabundfile = validParameter.validFile(parameters, "relabund", true);
150 if (relabundfile == "not open") { abort = true; }
151 else if (relabundfile == "not found") { relabundfile = ""; }
152 else { inputFileName = relabundfile; }
154 metadatafile = validParameter.validFile(parameters, "metadata", true);
155 if (metadatafile == "not open") { abort = true; }
156 else if (metadatafile == "not found") { metadatafile = ""; }
157 else { inputFileName = metadatafile; }
159 groups = validParameter.validFile(parameters, "groups", false);
160 if (groups == "not found") { groups = ""; pickedGroups = false; }
163 m->splitAtDash(groups, Groups);
165 globaldata->Groups = Groups;
167 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
169 label = validParameter.validFile(parameters, "label", false);
170 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=""; }
172 if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) { m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true; }
174 if (metadatafile != "") {
175 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
177 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
180 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
181 convert(temp, numaxes);
183 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
185 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; }
189 catch(exception& e) {
190 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
194 //**********************************************************************************************************************
196 void CorrAxesCommand::help(){
198 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");
199 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");
200 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");
201 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");
202 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");
203 m->mothurOut("The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n");
204 m->mothurOut("The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n");
205 m->mothurOut("Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n");
206 m->mothurOut("The corr.axes command outputs a .corr.axes file.\n");
207 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
209 catch(exception& e) {
210 m->errorOut(e, "CorrAxesCommand", "help");
215 //**********************************************************************************************************************
217 CorrAxesCommand::~CorrAxesCommand(){}
219 //**********************************************************************************************************************
221 int CorrAxesCommand::execute(){
224 if (abort == true) { return 0; }
226 /*************************************************************************************/
227 // use smart distancing to get right sharedRabund and convert to relabund if needed //
228 /************************************************************************************/
229 if (sharedfile != "") {
230 InputData* input = new InputData(sharedfile, "sharedfile");
231 getSharedFloat(input);
234 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
235 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
237 }else if (relabundfile != "") {
238 InputData* input = new InputData(relabundfile, "relabund");
239 getSharedFloat(input);
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 relabund file."); m->mothurOutEndLine(); return 0; }
245 }else if (metadatafile != "") {
246 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
247 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
248 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
250 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
252 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
254 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
256 //this is for a sanity check to make sure the axes file and shared file match
257 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
259 /*************************************************************************************/
260 // read axes file and check for file mismatches with shared or relabund file //
261 /************************************************************************************/
264 map<string, vector<float> > axes = readAxes();
266 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
268 //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
269 if (axes.size() != lookupFloat.size()) {
270 map<string, vector<float> >::iterator it;
271 for (int i = 0; i < lookupFloat.size(); i++) {
272 it = axes.find(lookupFloat[i]->getGroup());
273 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(); }
275 m->control_pressed = true;
278 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
280 /*************************************************************************************/
281 // calc the r values //
282 /************************************************************************************/
284 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
285 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
287 m->openOutputFile(outputFileName, out);
288 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
291 if (metadatafile == "") { out << "OTU"; }
292 else { out << "Feature"; }
294 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1); }
295 out << "\tlength" << endl;
297 if (method == "pearson") { calcPearson(axes, out); }
298 else if (method == "spearman") { calcSpearman(axes, out); }
299 else if (method == "kendall") { calcKendall(axes, out); }
300 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
303 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
305 if (m->control_pressed) { return 0; }
307 m->mothurOutEndLine();
308 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
309 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
310 m->mothurOutEndLine();
314 catch(exception& e) {
315 m->errorOut(e, "CorrAxesCommand", "execute");
319 //**********************************************************************************************************************
320 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
323 //find average of each axis - X
324 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
325 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
326 vector<float> temp = it->second;
327 for (int i = 0; i < temp.size(); i++) {
328 averageAxes[i] += temp[i];
332 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
335 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
337 if (metadatafile == "") { out << i+1; }
338 else { out << metadataLabels[i]; }
340 //find the averages this otu - Y
342 for (int j = 0; j < lookupFloat.size(); j++) {
343 sumOtu += lookupFloat[j]->getAbundance(i);
345 float Ybar = sumOtu / (float) lookupFloat.size();
347 vector<float> rValues(averageAxes.size());
349 //find r value for each axis
350 for (int k = 0; k < averageAxes.size(); k++) {
353 double numerator = 0.0;
354 double denomTerm1 = 0.0;
355 double denomTerm2 = 0.0;
356 for (int j = 0; j < lookupFloat.size(); j++) {
357 float Yi = lookupFloat[j]->getAbundance(i);
358 float Xi = axes[lookupFloat[j]->getGroup()][k];
360 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
361 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
362 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
365 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
367 r = numerator / denom;
373 for(int k=0;k<rValues.size();k++){
374 sum += rValues[k] * rValues[k];
376 out << '\t' << sqrt(sum) << endl;
381 catch(exception& e) {
382 m->errorOut(e, "CorrAxesCommand", "calcPearson");
386 //**********************************************************************************************************************
387 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
391 vector< map<float, int> > tableX; tableX.resize(numaxes);
392 map<float, int>::iterator itTable;
393 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
394 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
395 vector<float> temp = it->second;
396 for (int i = 0; i < temp.size(); i++) {
397 spearmanRank member(it->first, temp[i]);
398 scores[i].push_back(member);
400 //count number of repeats
401 itTable = tableX[i].find(temp[i]);
402 if (itTable == tableX[i].end()) {
403 tableX[i][temp[i]] = 1;
405 tableX[i][temp[i]]++;
412 vector<double> Lx; Lx.resize(numaxes, 0.0);
413 for (int i = 0; i < numaxes; i++) {
414 for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) {
415 double tx = (double) itTable->second;
416 Lx[i] += ((pow(tx, 3.0) - tx) / 12.0);
421 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
423 //find ranks of xi in each axis
424 map<string, vector<float> > rankAxes;
425 for (int i = 0; i < numaxes; i++) {
427 vector<spearmanRank> ties;
429 for (int j = 0; j < scores[i].size(); j++) {
431 ties.push_back(scores[i][j]);
433 if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
434 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
436 for (int k = 0; k < ties.size(); k++) {
437 float thisrank = rankTotal / (float) ties.size();
438 rankAxes[ties[k].name].push_back(thisrank);
443 }else { // you are the last one
445 for (int k = 0; k < ties.size(); k++) {
446 float thisrank = rankTotal / (float) ties.size();
447 rankAxes[ties[k].name].push_back(thisrank);
456 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
458 if (metadatafile == "") { out << i+1; }
459 else { out << metadataLabels[i]; }
461 //find the ranks of this otu - Y
462 vector<spearmanRank> otuScores;
463 map<float, int> tableY;
464 for (int j = 0; j < lookupFloat.size(); j++) {
465 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
466 otuScores.push_back(member);
468 itTable = tableY.find(member.score);
469 if (itTable == tableY.end()) {
470 tableY[member.score] = 1;
472 tableY[member.score]++;
479 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
480 double ty = (double) itTable->second;
481 Ly += ((pow(ty, 3.0) - ty) / 12.0);
484 sort(otuScores.begin(), otuScores.end(), compareSpearman);
486 map<string, float> rankOtus;
487 vector<spearmanRank> ties;
489 for (int j = 0; j < otuScores.size(); j++) {
491 ties.push_back(otuScores[j]);
493 if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
494 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
496 for (int k = 0; k < ties.size(); k++) {
497 float thisrank = rankTotal / (float) ties.size();
498 rankOtus[ties[k].name] = thisrank;
503 }else { // you are the last one
505 for (int k = 0; k < ties.size(); k++) {
506 float thisrank = rankTotal / (float) ties.size();
507 rankOtus[ties[k].name] = thisrank;
511 vector<double> pValues(numaxes);
513 //calc spearman ranks for each axis for this otu
514 for (int j = 0; j < numaxes; j++) {
517 for (int k = 0; k < lookupFloat.size(); k++) {
519 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
520 float yi = rankOtus[lookupFloat[k]->getGroup()];
522 di += ((xi - yi) * (xi - yi));
527 double n = (double) lookupFloat.size();
529 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
530 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
532 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
540 for(int k=0;k<numaxes;k++){
541 sum += pValues[k] * pValues[k];
543 out << '\t' << sqrt(sum) << endl;
548 catch(exception& e) {
549 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
553 //**********************************************************************************************************************
554 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
558 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
559 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
560 vector<float> temp = it->second;
561 for (int i = 0; i < temp.size(); i++) {
562 spearmanRank member(it->first, temp[i]);
563 scores[i].push_back(member);
568 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
570 //convert scores to ranks of xi in each axis
571 for (int i = 0; i < numaxes; i++) {
573 vector<spearmanRank*> ties;
575 for (int j = 0; j < scores[i].size(); j++) {
577 ties.push_back(&(scores[i][j]));
579 if (j != scores[i].size()-1) { // you are not the last so you can look ahead
580 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
581 for (int k = 0; k < ties.size(); k++) {
582 float thisrank = rankTotal / (float) ties.size();
583 (*ties[k]).score = thisrank;
588 }else { // you are the last one
589 for (int k = 0; k < ties.size(); k++) {
590 float thisrank = rankTotal / (float) ties.size();
591 (*ties[k]).score = thisrank;
598 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
600 if (metadatafile == "") { out << i+1; }
601 else { out << metadataLabels[i]; }
603 //find the ranks of this otu - Y
604 vector<spearmanRank> otuScores;
605 for (int j = 0; j < lookupFloat.size(); j++) {
606 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
607 otuScores.push_back(member);
610 sort(otuScores.begin(), otuScores.end(), compareSpearman);
612 map<string, float> rankOtus;
613 vector<spearmanRank> ties;
615 for (int j = 0; j < otuScores.size(); j++) {
617 ties.push_back(otuScores[j]);
619 if (j != otuScores.size()-1) { // you are not the last so you can look ahead
620 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
621 for (int k = 0; k < ties.size(); k++) {
622 float thisrank = rankTotal / (float) ties.size();
623 rankOtus[ties[k].name] = thisrank;
628 }else { // you are the last one
629 for (int k = 0; k < ties.size(); k++) {
630 float thisrank = rankTotal / (float) ties.size();
631 rankOtus[ties[k].name] = thisrank;
637 vector<double> pValues(numaxes);
639 //calc spearman ranks for each axis for this otu
640 for (int j = 0; j < numaxes; j++) {
645 vector<spearmanRank> otus;
646 vector<spearmanRank> otusTemp;
647 for (int l = 0; l < scores[j].size(); l++) {
648 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
649 otus.push_back(member);
653 for (int l = 0; l < scores[j].size(); l++) {
655 int numWithHigherRank = 0;
656 int numWithLowerRank = 0;
657 float thisrank = otus[l].score;
659 for (int u = l; u < scores[j].size(); u++) {
660 if (otus[u].score > thisrank) { numWithHigherRank++; }
661 else if (otus[u].score < thisrank) { numWithLowerRank++; }
665 numCoor += numWithHigherRank;
666 numDisCoor += numWithLowerRank;
669 //comparing to yourself
670 count -= lookupFloat.size();
672 double p = (numCoor - numDisCoor) / (float) count;
680 for(int k=0;k<numaxes;k++){
681 sum += pValues[k] * pValues[k];
683 out << '\t' << sqrt(sum) << endl;
688 catch(exception& e) {
689 m->errorOut(e, "CorrAxesCommand", "calcKendall");
693 //**********************************************************************************************************************
694 int CorrAxesCommand::getSharedFloat(InputData* input){
696 lookupFloat = input->getSharedRAbundFloatVectors();
697 string lastLabel = lookupFloat[0]->getLabel();
699 if (label == "") { label = lastLabel; return 0; }
701 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
702 set<string> labels; labels.insert(label);
703 set<string> processedLabels;
704 set<string> userLabels = labels;
706 //as long as you are not at the end of the file or done wih the lines you want
707 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
709 if (m->control_pressed) { return 0; }
711 if(labels.count(lookupFloat[0]->getLabel()) == 1){
712 processedLabels.insert(lookupFloat[0]->getLabel());
713 userLabels.erase(lookupFloat[0]->getLabel());
717 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
718 string saveLabel = lookupFloat[0]->getLabel();
720 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
721 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
723 processedLabels.insert(lookupFloat[0]->getLabel());
724 userLabels.erase(lookupFloat[0]->getLabel());
726 //restore real lastlabel to save below
727 lookupFloat[0]->setLabel(saveLabel);
731 lastLabel = lookupFloat[0]->getLabel();
733 //get next line to process
734 //prevent memory leak
735 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
736 lookupFloat = input->getSharedRAbundFloatVectors();
740 if (m->control_pressed) { return 0; }
742 //output error messages about any remaining user labels
743 set<string>::iterator it;
744 bool needToRun = false;
745 for (it = userLabels.begin(); it != userLabels.end(); it++) {
746 m->mothurOut("Your file does not include the label " + *it);
747 if (processedLabels.count(lastLabel) != 1) {
748 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
751 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
755 //run last label if you need to
756 if (needToRun == true) {
757 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
758 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
763 catch(exception& e) {
764 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
768 //**********************************************************************************************************************
769 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
772 vector<SharedRAbundFloatVector*> newLookup;
773 for (int i = 0; i < thislookup.size(); i++) {
774 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
775 temp->setLabel(thislookup[i]->getLabel());
776 temp->setGroup(thislookup[i]->getGroup());
777 newLookup.push_back(temp);
781 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
782 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
784 //look at each sharedRabund and make sure they are not all zero
786 for (int j = 0; j < thislookup.size(); j++) {
787 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
790 //if they are not all zero add this bin
792 for (int j = 0; j < thislookup.size(); j++) {
793 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
798 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
800 thislookup = newLookup;
805 catch(exception& e) {
806 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
810 /*****************************************************************/
811 map<string, vector<float> > CorrAxesCommand::readAxes(){
813 map<string, vector<float> > axes;
816 m->openInputFile(axesfile, in);
818 string headerLine = m->getline(in); m->gobble(in);
820 //count the number of axis you are reading
824 int pos = headerLine.find("axis");
825 if (pos != string::npos) {
827 headerLine = headerLine.substr(pos+4);
828 }else { done = true; }
831 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
835 if (m->control_pressed) { in.close(); return axes; }
838 in >> group; m->gobble(in);
840 vector<float> thisGroupsAxes;
841 for (int i = 0; i < count; i++) {
845 //only save the axis we want
846 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
849 //save group if its one the user selected
850 if (names.count(group) != 0) {
851 map<string, vector<float> >::iterator it = axes.find(group);
853 if (it == axes.end()) {
854 axes[group] = thisGroupsAxes;
856 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
866 catch(exception& e) {
867 m->errorOut(e, "CorrAxesCommand", "readAxes");
871 /*****************************************************************/
872 int CorrAxesCommand::getMetadata(){
874 vector<string> groupNames;
877 m->openInputFile(axesfile, in);
879 string headerLine = m->getline(in); m->gobble(in);
880 istringstream iss (headerLine,istringstream::in);
882 //read the first label, because it refers to the groups
884 iss >> columnLabel; m->gobble(iss);
886 //save names of columns you are reading
888 iss >> columnLabel; m->gobble(iss);
889 metadataLabels.push_back(columnLabel);
891 int count = metadataLabels.size();
896 if (m->control_pressed) { in.close(); return 0; }
899 in >> group; m->gobble(in);
900 groupNames.push_back(group);
902 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
903 tempLookup->setGroup(group);
904 tempLookup->setLabel("1");
906 for (int i = 0; i < count; i++) {
910 tempLookup->push_back(temp, group);
913 lookupFloat.push_back(tempLookup);
919 //remove any groups the user does not want, and set globaldata->groups with only valid groups
921 util = new SharedUtil();
923 util->setGroups(globaldata->Groups, groupNames);
925 for (int i = 0; i < lookupFloat.size(); i++) {
926 //if this sharedrabund is not from a group the user wants then delete it.
927 if (util->isValidGroup(lookupFloat[i]->getGroup(), globaldata->Groups) == false) {
928 delete lookupFloat[i]; lookupFloat[i] = NULL;
929 lookupFloat.erase(lookupFloat.begin()+i);
938 catch(exception& e) {
939 m->errorOut(e, "CorrAxesCommand", "getMetadata");
943 /*****************************************************************/