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 or relabund file as well as a pcoa 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 != "") {
226 if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; }
227 if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
229 //fills lookupFloat with relative abundance values from lookup
232 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
234 }else if (relabundfile != "") {
236 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
237 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
239 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
241 }else if (metadatafile != "") {
242 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
243 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
244 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
246 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
248 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
250 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
252 //this is for a sanity check to make sure the axes file and shared file match
253 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
255 /*************************************************************************************/
256 // read axes file and check for file mismatches with shared or relabund file //
257 /************************************************************************************/
260 map<string, vector<float> > axes = readAxes();
262 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
264 //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
265 if (axes.size() != lookupFloat.size()) {
266 map<string, vector<float> >::iterator it;
267 for (int i = 0; i < lookupFloat.size(); i++) {
268 it = axes.find(lookupFloat[i]->getGroup());
269 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(); }
271 m->control_pressed = true;
274 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
276 /*************************************************************************************/
277 // calc the r values //
278 /************************************************************************************/
280 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
281 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
283 m->openOutputFile(outputFileName, out);
284 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
287 if (metadatafile == "") { out << "OTU\t"; }
288 else { out << "Feature\t"; }
290 for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; }
293 if (method == "pearson") { calcPearson(axes, out); }
294 else if (method == "spearman") { calcSpearman(axes, out); }
295 else if (method == "kendall") { calcKendall(axes, out); }
296 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
299 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
301 if (m->control_pressed) { return 0; }
303 m->mothurOutEndLine();
304 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
305 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
306 m->mothurOutEndLine();
310 catch(exception& e) {
311 m->errorOut(e, "CorrAxesCommand", "execute");
315 //**********************************************************************************************************************
316 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
319 //find average of each axis - X
320 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
321 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
322 vector<float> temp = it->second;
323 for (int i = 0; i < temp.size(); i++) {
324 averageAxes[i] += temp[i];
328 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
331 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
333 if (metadatafile == "") { out << i+1 << '\t'; }
334 else { out << metadataLabels[i] << '\t'; }
336 //find the averages this otu - Y
338 for (int j = 0; j < lookupFloat.size(); j++) {
339 sumOtu += lookupFloat[j]->getAbundance(i);
341 float Ybar = sumOtu / (float) lookupFloat.size();
343 //find r value for each axis
344 for (int k = 0; k < averageAxes.size(); k++) {
347 double numerator = 0.0;
348 double denomTerm1 = 0.0;
349 double denomTerm2 = 0.0;
350 for (int j = 0; j < lookupFloat.size(); j++) {
351 float Yi = lookupFloat[j]->getAbundance(i);
352 float Xi = axes[lookupFloat[j]->getGroup()][k];
354 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
355 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
356 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
359 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
361 r = numerator / denom;
371 catch(exception& e) {
372 m->errorOut(e, "CorrAxesCommand", "calcPearson");
376 //**********************************************************************************************************************
377 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
381 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
382 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
383 vector<float> temp = it->second;
384 for (int i = 0; i < temp.size(); i++) {
385 spearmanRank member(it->first, temp[i]);
386 scores[i].push_back(member);
391 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
393 //find ranks of xi in each axis
394 map<string, vector<float> > rankAxes;
395 for (int i = 0; i < numaxes; i++) {
397 vector<spearmanRank> ties;
399 for (int j = 0; j < scores[i].size(); j++) {
401 ties.push_back(scores[i][j]);
403 if (j != scores.size()-1) { // you are not the last so you can look ahead
404 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
405 for (int k = 0; k < ties.size(); k++) {
406 float thisrank = rankTotal / (float) ties.size();
407 rankAxes[ties[k].name].push_back(thisrank);
412 }else { // you are the last one
413 for (int k = 0; k < ties.size(); k++) {
414 float thisrank = rankTotal / (float) ties.size();
415 rankAxes[ties[k].name].push_back(thisrank);
424 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
426 if (metadatafile == "") { out << i+1 << '\t'; }
427 else { out << metadataLabels[i] << '\t'; }
429 //find the ranks of this otu - Y
430 vector<spearmanRank> otuScores;
431 for (int j = 0; j < lookupFloat.size(); j++) {
432 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
433 otuScores.push_back(member);
436 sort(otuScores.begin(), otuScores.end(), compareSpearman);
438 map<string, float> rankOtus;
439 vector<spearmanRank> ties;
441 for (int j = 0; j < otuScores.size(); j++) {
443 ties.push_back(otuScores[j]);
445 if (j != scores.size()-1) { // you are not the last so you can look ahead
446 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
447 for (int k = 0; k < ties.size(); k++) {
448 float thisrank = rankTotal / (float) ties.size();
449 rankOtus[ties[k].name] = thisrank;
454 }else { // you are the last one
455 for (int k = 0; k < ties.size(); k++) {
456 float thisrank = rankTotal / (float) ties.size();
457 rankOtus[ties[k].name] = thisrank;
462 //calc spearman ranks for each axis for this otu
463 for (int j = 0; j < numaxes; j++) {
466 for (int k = 0; k < lookupFloat.size(); k++) {
468 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
469 float yi = rankOtus[lookupFloat[k]->getGroup()];
471 di += ((xi - yi) * (xi - yi));
474 int n = lookupFloat.size();
475 double p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1)));
486 catch(exception& e) {
487 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
491 //**********************************************************************************************************************
492 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
496 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
497 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
498 vector<float> temp = it->second;
499 for (int i = 0; i < temp.size(); i++) {
500 spearmanRank member(it->first, temp[i]);
501 scores[i].push_back(member);
506 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
508 //convert scores to ranks of xi in each axis
509 for (int i = 0; i < numaxes; i++) {
511 vector<spearmanRank*> ties;
513 for (int j = 0; j < scores[i].size(); j++) {
515 ties.push_back(&(scores[i][j]));
517 if (j != scores.size()-1) { // you are not the last so you can look ahead
518 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
519 for (int k = 0; k < ties.size(); k++) {
520 float thisrank = rankTotal / (float) ties.size();
521 (*ties[k]).score = thisrank;
526 }else { // you are the last one
527 for (int k = 0; k < ties.size(); k++) {
528 float thisrank = rankTotal / (float) ties.size();
529 (*ties[k]).score = thisrank;
536 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
538 if (metadatafile == "") { out << i+1 << '\t'; }
539 else { out << metadataLabels[i] << '\t'; }
541 //find the ranks of this otu - Y
542 vector<spearmanRank> otuScores;
543 for (int j = 0; j < lookupFloat.size(); j++) {
544 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
545 otuScores.push_back(member);
548 sort(otuScores.begin(), otuScores.end(), compareSpearman);
550 map<string, float> rankOtus;
551 vector<spearmanRank> ties;
553 for (int j = 0; j < otuScores.size(); j++) {
555 ties.push_back(otuScores[j]);
557 if (j != scores.size()-1) { // you are not the last so you can look ahead
558 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
559 for (int k = 0; k < ties.size(); k++) {
560 float thisrank = rankTotal / (float) ties.size();
561 rankOtus[ties[k].name] = thisrank;
566 }else { // you are the last one
567 for (int k = 0; k < ties.size(); k++) {
568 float thisrank = rankTotal / (float) ties.size();
569 rankOtus[ties[k].name] = thisrank;
574 //calc spearman ranks for each axis for this otu
575 for (int j = 0; j < numaxes; j++) {
578 //assemble otus ranks in same order as axis ranks
579 vector<spearmanRank> otus;
580 for (int l = 0; l < scores[j].size(); l++) {
581 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
582 otus.push_back(member);
585 for (int l = 0; l < scores[j].size(); l++) {
587 int numWithHigherRank = 0;
588 float thisrank = otus[l].score;
590 for (int u = l; u < scores[j].size(); u++) {
591 if (otus[u].score > thisrank) { numWithHigherRank++; }
594 P += numWithHigherRank;
597 int n = lookupFloat.size();
599 double p = ( (4 * P) / (float) (n * (n - 1)) ) - 1.0;
609 catch(exception& e) {
610 m->errorOut(e, "CorrAxesCommand", "calcKendall");
614 //**********************************************************************************************************************
615 int CorrAxesCommand::getShared(){
617 InputData* input = new InputData(sharedfile, "sharedfile");
618 lookup = input->getSharedRAbundVectors();
619 string lastLabel = lookup[0]->getLabel();
621 if (label == "") { label = lastLabel; delete input; return 0; }
623 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
624 set<string> labels; labels.insert(label);
625 set<string> processedLabels;
626 set<string> userLabels = labels;
628 //as long as you are not at the end of the file or done wih the lines you want
629 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
630 if (m->control_pressed) { delete input; return 0; }
632 if(labels.count(lookup[0]->getLabel()) == 1){
633 processedLabels.insert(lookup[0]->getLabel());
634 userLabels.erase(lookup[0]->getLabel());
638 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
639 string saveLabel = lookup[0]->getLabel();
641 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
642 lookup = input->getSharedRAbundVectors(lastLabel);
644 processedLabels.insert(lookup[0]->getLabel());
645 userLabels.erase(lookup[0]->getLabel());
647 //restore real lastlabel to save below
648 lookup[0]->setLabel(saveLabel);
652 lastLabel = lookup[0]->getLabel();
654 //get next line to process
655 //prevent memory leak
656 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
657 lookup = input->getSharedRAbundVectors();
661 if (m->control_pressed) { delete input; return 0; }
663 //output error messages about any remaining user labels
664 set<string>::iterator it;
665 bool needToRun = false;
666 for (it = userLabels.begin(); it != userLabels.end(); it++) {
667 m->mothurOut("Your file does not include the label " + *it);
668 if (processedLabels.count(lastLabel) != 1) {
669 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
672 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
676 //run last label if you need to
677 if (needToRun == true) {
678 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
679 lookup = input->getSharedRAbundVectors(lastLabel);
685 catch(exception& e) {
686 m->errorOut(e, "CorrAxesCommand", "getShared");
690 //**********************************************************************************************************************
691 int CorrAxesCommand::getSharedFloat(){
693 InputData* input = new InputData(relabundfile, "relabund");
694 lookupFloat = input->getSharedRAbundFloatVectors();
695 string lastLabel = lookupFloat[0]->getLabel();
697 if (label == "") { label = lastLabel; delete input; return 0; }
699 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
700 set<string> labels; labels.insert(label);
701 set<string> processedLabels;
702 set<string> userLabels = labels;
704 //as long as you are not at the end of the file or done wih the lines you want
705 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
707 if (m->control_pressed) { delete input; return 0; }
709 if(labels.count(lookupFloat[0]->getLabel()) == 1){
710 processedLabels.insert(lookupFloat[0]->getLabel());
711 userLabels.erase(lookupFloat[0]->getLabel());
715 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
716 string saveLabel = lookupFloat[0]->getLabel();
718 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
719 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
721 processedLabels.insert(lookupFloat[0]->getLabel());
722 userLabels.erase(lookupFloat[0]->getLabel());
724 //restore real lastlabel to save below
725 lookupFloat[0]->setLabel(saveLabel);
729 lastLabel = lookupFloat[0]->getLabel();
731 //get next line to process
732 //prevent memory leak
733 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
734 lookupFloat = input->getSharedRAbundFloatVectors();
738 if (m->control_pressed) { delete input; return 0; }
740 //output error messages about any remaining user labels
741 set<string>::iterator it;
742 bool needToRun = false;
743 for (it = userLabels.begin(); it != userLabels.end(); it++) {
744 m->mothurOut("Your file does not include the label " + *it);
745 if (processedLabels.count(lastLabel) != 1) {
746 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
749 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
753 //run last label if you need to
754 if (needToRun == true) {
755 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
756 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
762 catch(exception& e) {
763 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
767 /*****************************************************************/
768 int CorrAxesCommand::convertToRelabund(){
771 vector<SharedRAbundFloatVector*> newLookup;
772 for (int i = 0; i < lookup.size(); i++) {
773 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
774 temp->setLabel(lookup[i]->getLabel());
775 temp->setGroup(lookup[i]->getGroup());
776 newLookup.push_back(temp);
779 for (int i = 0; i < lookup.size(); i++) {
781 for (int j = 0; j < lookup[i]->getNumBins(); j++) {
783 if (m->control_pressed) { return 0; }
785 int abund = lookup[i]->getAbundance(j);
787 float relabund = abund / (float) lookup[i]->getNumSeqs();
789 newLookup[i]->push_back(relabund, lookup[i]->getGroup());
793 if (pickedGroups) { eliminateZeroOTUS(newLookup); }
795 lookupFloat = newLookup;
799 catch(exception& e) {
800 m->errorOut(e, "CorrAxesCommand", "convertToRelabund");
804 //**********************************************************************************************************************
805 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
808 vector<SharedRAbundFloatVector*> newLookup;
809 for (int i = 0; i < thislookup.size(); i++) {
810 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
811 temp->setLabel(thislookup[i]->getLabel());
812 temp->setGroup(thislookup[i]->getGroup());
813 newLookup.push_back(temp);
817 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
818 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
820 //look at each sharedRabund and make sure they are not all zero
822 for (int j = 0; j < thislookup.size(); j++) {
823 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
826 //if they are not all zero add this bin
828 for (int j = 0; j < thislookup.size(); j++) {
829 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
834 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
836 thislookup = newLookup;
841 catch(exception& e) {
842 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
846 /*****************************************************************/
847 map<string, vector<float> > CorrAxesCommand::readAxes(){
849 map<string, vector<float> > axes;
852 m->openInputFile(axesfile, in);
854 string headerLine = m->getline(in); m->gobble(in);
856 //count the number of axis you are reading
860 int pos = headerLine.find("axis");
861 if (pos != string::npos) {
863 headerLine = headerLine.substr(pos+4);
864 }else { done = true; }
867 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
871 if (m->control_pressed) { in.close(); return axes; }
874 in >> group; m->gobble(in);
876 vector<float> thisGroupsAxes;
877 for (int i = 0; i < count; i++) {
881 //only save the axis we want
882 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
885 //save group if its one the user selected
886 if (names.count(group) != 0) {
887 map<string, vector<float> >::iterator it = axes.find(group);
889 if (it == axes.end()) {
890 axes[group] = thisGroupsAxes;
892 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
902 catch(exception& e) {
903 m->errorOut(e, "CorrAxesCommand", "readAxes");
907 /*****************************************************************/
908 int CorrAxesCommand::getMetadata(){
910 vector<string> groupNames;
913 m->openInputFile(axesfile, in);
915 string headerLine = m->getline(in); m->gobble(in);
916 istringstream iss (headerLine,istringstream::in);
918 //read the first label, because it refers to the groups
920 iss >> columnLabel; m->gobble(iss);
922 //save names of columns you are reading
924 iss >> columnLabel; m->gobble(iss);
925 metadataLabels.push_back(columnLabel);
927 int count = metadataLabels.size();
932 if (m->control_pressed) { in.close(); return 0; }
935 in >> group; m->gobble(in);
936 groupNames.push_back(group);
938 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
939 tempLookup->setGroup(group);
940 tempLookup->setLabel("1");
942 for (int i = 0; i < count; i++) {
946 tempLookup->push_back(temp, group);
949 lookupFloat.push_back(tempLookup);
955 //remove any groups the user does not want, and set globaldata->groups with only valid groups
957 util = new SharedUtil();
959 util->setGroups(globaldata->Groups, groupNames);
961 for (int i = 0; i < lookupFloat.size(); i++) {
962 //if this sharedrabund is not from a group the user wants then delete it.
963 if (util->isValidGroup(lookupFloat[i]->getGroup(), globaldata->Groups) == false) {
964 delete lookupFloat[i]; lookupFloat[i] = NULL;
965 lookupFloat.erase(lookupFloat.begin()+i);
974 catch(exception& e) {
975 m->errorOut(e, "CorrAxesCommand", "getMetadata");
979 /*****************************************************************/