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 vector<string> CorrAxesCommand::getValidParameters(){
16 string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
17 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
21 m->errorOut(e, "CorrAxesCommand", "getValidParameters");
25 //**********************************************************************************************************************
26 vector<string> CorrAxesCommand::getRequiredParameters(){
28 string Array[] = {"axes"};
29 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
33 m->errorOut(e, "CorrAxesCommand", "getRequiredParameters");
37 //**********************************************************************************************************************
38 CorrAxesCommand::CorrAxesCommand(){
40 abort = true; calledHelp = true;
41 vector<string> tempOutNames;
42 outputTypes["corr.axes"] = tempOutNames;
45 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
50 //**********************************************************************************************************************
51 vector<string> CorrAxesCommand::getRequiredFiles(){
53 vector<string> myArray;
57 m->errorOut(e, "CorrAxesCommand", "getRequiredFiles");
61 //**********************************************************************************************************************
62 CorrAxesCommand::CorrAxesCommand(string option) {
64 abort = false; calledHelp = false;
65 globaldata = GlobalData::getInstance();
67 //allow user to run help
68 if(option == "help") { help(); abort = true; calledHelp = true; }
71 //valid paramters for this command
72 string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
73 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
75 OptionParser parser(option);
76 map<string, string> parameters = parser.getParameters();
78 ValidParameters validParameter;
79 map<string, string>::iterator it;
81 //check to make sure all parameters are valid for command
82 for (it = parameters.begin(); it != parameters.end(); it++) {
83 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
86 vector<string> tempOutNames;
87 outputTypes["corr.axes"] = tempOutNames;
89 //if the user changes the input directory command factory will send this info to us in the output parameter
90 string inputDir = validParameter.validFile(parameters, "inputdir", false);
91 if (inputDir == "not found"){ inputDir = ""; }
94 it = parameters.find("axes");
95 //user has given a template file
96 if(it != parameters.end()){
97 path = m->hasPath(it->second);
98 //if the user has not given a path then, add inputdir. else leave path alone.
99 if (path == "") { parameters["axes"] = inputDir + it->second; }
102 it = parameters.find("shared");
103 //user has given a template file
104 if(it != parameters.end()){
105 path = m->hasPath(it->second);
106 //if the user has not given a path then, add inputdir. else leave path alone.
107 if (path == "") { parameters["shared"] = inputDir + it->second; }
110 it = parameters.find("relabund");
111 //user has given a template file
112 if(it != parameters.end()){
113 path = m->hasPath(it->second);
114 //if the user has not given a path then, add inputdir. else leave path alone.
115 if (path == "") { parameters["relabund"] = inputDir + it->second; }
118 it = parameters.find("metadata");
119 //user has given a template file
120 if(it != parameters.end()){
121 path = m->hasPath(it->second);
122 //if the user has not given a path then, add inputdir. else leave path alone.
123 if (path == "") { parameters["metadata"] = inputDir + it->second; }
128 //check for required parameters
129 axesfile = validParameter.validFile(parameters, "axes", true);
130 if (axesfile == "not open") { abort = true; }
131 else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true; }
133 sharedfile = validParameter.validFile(parameters, "shared", true);
134 if (sharedfile == "not open") { abort = true; }
135 else if (sharedfile == "not found") { sharedfile = ""; }
136 else { inputFileName = sharedfile; }
138 relabundfile = validParameter.validFile(parameters, "relabund", true);
139 if (relabundfile == "not open") { abort = true; }
140 else if (relabundfile == "not found") { relabundfile = ""; }
141 else { inputFileName = relabundfile; }
143 metadatafile = validParameter.validFile(parameters, "metadata", true);
144 if (metadatafile == "not open") { abort = true; }
145 else if (metadatafile == "not found") { metadatafile = ""; }
146 else { inputFileName = metadatafile; }
148 groups = validParameter.validFile(parameters, "groups", false);
149 if (groups == "not found") { groups = ""; pickedGroups = false; }
152 m->splitAtDash(groups, Groups);
154 globaldata->Groups = Groups;
156 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
158 label = validParameter.validFile(parameters, "label", false);
159 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=""; }
161 if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) { m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true; }
163 if (metadatafile != "") {
164 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
166 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
169 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
170 convert(temp, numaxes);
172 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
174 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; }
178 catch(exception& e) {
179 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
183 //**********************************************************************************************************************
185 void CorrAxesCommand::help(){
187 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");
188 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");
189 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");
190 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");
191 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");
192 m->mothurOut("The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n");
193 m->mothurOut("The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n");
194 m->mothurOut("Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n");
195 m->mothurOut("The corr.axes command outputs a .corr.axes file.\n");
196 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
198 catch(exception& e) {
199 m->errorOut(e, "CorrAxesCommand", "help");
204 //**********************************************************************************************************************
206 CorrAxesCommand::~CorrAxesCommand(){}
208 //**********************************************************************************************************************
210 int CorrAxesCommand::execute(){
213 if (abort == true) { if (calledHelp) { return 0; } return 2; }
215 /*************************************************************************************/
216 // use smart distancing to get right sharedRabund and convert to relabund if needed //
217 /************************************************************************************/
218 if (sharedfile != "") {
219 InputData* input = new InputData(sharedfile, "sharedfile");
220 getSharedFloat(input);
223 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
224 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
226 }else if (relabundfile != "") {
227 InputData* input = new InputData(relabundfile, "relabund");
228 getSharedFloat(input);
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 }else if (metadatafile != "") {
235 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
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 metadata file."); m->mothurOutEndLine(); return 0; }
239 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
241 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
243 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
245 //this is for a sanity check to make sure the axes file and shared file match
246 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
248 /*************************************************************************************/
249 // read axes file and check for file mismatches with shared or relabund file //
250 /************************************************************************************/
253 map<string, vector<float> > axes = readAxes();
255 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
257 //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
258 if (axes.size() != lookupFloat.size()) {
259 map<string, vector<float> >::iterator it;
260 for (int i = 0; i < lookupFloat.size(); i++) {
261 it = axes.find(lookupFloat[i]->getGroup());
262 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(); }
264 m->control_pressed = true;
267 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
269 /*************************************************************************************/
270 // calc the r values //
271 /************************************************************************************/
273 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
274 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
276 m->openOutputFile(outputFileName, out);
277 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
280 if (metadatafile == "") { out << "OTU"; }
281 else { out << "Feature"; }
283 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1); }
284 out << "\tlength" << endl;
286 if (method == "pearson") { calcPearson(axes, out); }
287 else if (method == "spearman") { calcSpearman(axes, out); }
288 else if (method == "kendall") { calcKendall(axes, out); }
289 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
292 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
294 if (m->control_pressed) { return 0; }
296 m->mothurOutEndLine();
297 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
298 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
299 m->mothurOutEndLine();
303 catch(exception& e) {
304 m->errorOut(e, "CorrAxesCommand", "execute");
308 //**********************************************************************************************************************
309 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
312 //find average of each axis - X
313 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
314 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
315 vector<float> temp = it->second;
316 for (int i = 0; i < temp.size(); i++) {
317 averageAxes[i] += temp[i];
321 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
324 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
326 if (metadatafile == "") { out << i+1; }
327 else { out << metadataLabels[i]; }
329 //find the averages this otu - Y
331 for (int j = 0; j < lookupFloat.size(); j++) {
332 sumOtu += lookupFloat[j]->getAbundance(i);
334 float Ybar = sumOtu / (float) lookupFloat.size();
336 vector<float> rValues(averageAxes.size());
338 //find r value for each axis
339 for (int k = 0; k < averageAxes.size(); k++) {
342 double numerator = 0.0;
343 double denomTerm1 = 0.0;
344 double denomTerm2 = 0.0;
345 for (int j = 0; j < lookupFloat.size(); j++) {
346 float Yi = lookupFloat[j]->getAbundance(i);
347 float Xi = axes[lookupFloat[j]->getGroup()][k];
349 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
350 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
351 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
354 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
356 r = numerator / denom;
362 for(int k=0;k<rValues.size();k++){
363 sum += rValues[k] * rValues[k];
365 out << '\t' << sqrt(sum) << endl;
370 catch(exception& e) {
371 m->errorOut(e, "CorrAxesCommand", "calcPearson");
375 //**********************************************************************************************************************
376 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
380 vector< map<float, int> > tableX; tableX.resize(numaxes);
381 map<float, int>::iterator itTable;
382 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
383 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
384 vector<float> temp = it->second;
385 for (int i = 0; i < temp.size(); i++) {
386 spearmanRank member(it->first, temp[i]);
387 scores[i].push_back(member);
389 //count number of repeats
390 itTable = tableX[i].find(temp[i]);
391 if (itTable == tableX[i].end()) {
392 tableX[i][temp[i]] = 1;
394 tableX[i][temp[i]]++;
401 vector<double> Lx; Lx.resize(numaxes, 0.0);
402 for (int i = 0; i < numaxes; i++) {
403 for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) {
404 double tx = (double) itTable->second;
405 Lx[i] += ((pow(tx, 3.0) - tx) / 12.0);
410 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
412 //find ranks of xi in each axis
413 map<string, vector<float> > rankAxes;
414 for (int i = 0; i < numaxes; i++) {
416 vector<spearmanRank> ties;
418 for (int j = 0; j < scores[i].size(); j++) {
420 ties.push_back(scores[i][j]);
422 if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
423 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
425 for (int k = 0; k < ties.size(); k++) {
426 float thisrank = rankTotal / (float) ties.size();
427 rankAxes[ties[k].name].push_back(thisrank);
432 }else { // you are the last one
434 for (int k = 0; k < ties.size(); k++) {
435 float thisrank = rankTotal / (float) ties.size();
436 rankAxes[ties[k].name].push_back(thisrank);
445 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
447 if (metadatafile == "") { out << i+1; }
448 else { out << metadataLabels[i]; }
450 //find the ranks of this otu - Y
451 vector<spearmanRank> otuScores;
452 map<float, int> tableY;
453 for (int j = 0; j < lookupFloat.size(); j++) {
454 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
455 otuScores.push_back(member);
457 itTable = tableY.find(member.score);
458 if (itTable == tableY.end()) {
459 tableY[member.score] = 1;
461 tableY[member.score]++;
468 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
469 double ty = (double) itTable->second;
470 Ly += ((pow(ty, 3.0) - ty) / 12.0);
473 sort(otuScores.begin(), otuScores.end(), compareSpearman);
475 map<string, float> rankOtus;
476 vector<spearmanRank> ties;
478 for (int j = 0; j < otuScores.size(); j++) {
480 ties.push_back(otuScores[j]);
482 if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
483 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
485 for (int k = 0; k < ties.size(); k++) {
486 float thisrank = rankTotal / (float) ties.size();
487 rankOtus[ties[k].name] = thisrank;
492 }else { // you are the last one
494 for (int k = 0; k < ties.size(); k++) {
495 float thisrank = rankTotal / (float) ties.size();
496 rankOtus[ties[k].name] = thisrank;
500 vector<double> pValues(numaxes);
502 //calc spearman ranks for each axis for this otu
503 for (int j = 0; j < numaxes; j++) {
506 for (int k = 0; k < lookupFloat.size(); k++) {
508 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
509 float yi = rankOtus[lookupFloat[k]->getGroup()];
511 di += ((xi - yi) * (xi - yi));
516 double n = (double) lookupFloat.size();
518 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
519 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
521 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
529 for(int k=0;k<numaxes;k++){
530 sum += pValues[k] * pValues[k];
532 out << '\t' << sqrt(sum) << endl;
537 catch(exception& e) {
538 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
542 //**********************************************************************************************************************
543 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
547 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
548 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
549 vector<float> temp = it->second;
550 for (int i = 0; i < temp.size(); i++) {
551 spearmanRank member(it->first, temp[i]);
552 scores[i].push_back(member);
557 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
559 //convert scores to ranks of xi in each axis
560 for (int i = 0; i < numaxes; i++) {
562 vector<spearmanRank*> ties;
564 for (int j = 0; j < scores[i].size(); j++) {
566 ties.push_back(&(scores[i][j]));
568 if (j != scores[i].size()-1) { // you are not the last so you can look ahead
569 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
570 for (int k = 0; k < ties.size(); k++) {
571 float thisrank = rankTotal / (float) ties.size();
572 (*ties[k]).score = thisrank;
577 }else { // you are the last one
578 for (int k = 0; k < ties.size(); k++) {
579 float thisrank = rankTotal / (float) ties.size();
580 (*ties[k]).score = thisrank;
587 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
589 if (metadatafile == "") { out << i+1; }
590 else { out << metadataLabels[i]; }
592 //find the ranks of this otu - Y
593 vector<spearmanRank> otuScores;
594 for (int j = 0; j < lookupFloat.size(); j++) {
595 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
596 otuScores.push_back(member);
599 sort(otuScores.begin(), otuScores.end(), compareSpearman);
601 map<string, float> rankOtus;
602 vector<spearmanRank> ties;
604 for (int j = 0; j < otuScores.size(); j++) {
606 ties.push_back(otuScores[j]);
608 if (j != otuScores.size()-1) { // you are not the last so you can look ahead
609 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
610 for (int k = 0; k < ties.size(); k++) {
611 float thisrank = rankTotal / (float) ties.size();
612 rankOtus[ties[k].name] = thisrank;
617 }else { // you are the last one
618 for (int k = 0; k < ties.size(); k++) {
619 float thisrank = rankTotal / (float) ties.size();
620 rankOtus[ties[k].name] = thisrank;
626 vector<double> pValues(numaxes);
628 //calc spearman ranks for each axis for this otu
629 for (int j = 0; j < numaxes; j++) {
634 vector<spearmanRank> otus;
635 vector<spearmanRank> otusTemp;
636 for (int l = 0; l < scores[j].size(); l++) {
637 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
638 otus.push_back(member);
642 for (int l = 0; l < scores[j].size(); l++) {
644 int numWithHigherRank = 0;
645 int numWithLowerRank = 0;
646 float thisrank = otus[l].score;
648 for (int u = l+1; u < scores[j].size(); u++) {
649 if (otus[u].score > thisrank) { numWithHigherRank++; }
650 else if (otus[u].score < thisrank) { numWithLowerRank++; }
654 numCoor += numWithHigherRank;
655 numDisCoor += numWithLowerRank;
658 double p = (numCoor - numDisCoor) / (float) count;
666 for(int k=0;k<numaxes;k++){
667 sum += pValues[k] * pValues[k];
669 out << '\t' << sqrt(sum) << endl;
674 catch(exception& e) {
675 m->errorOut(e, "CorrAxesCommand", "calcKendall");
679 //**********************************************************************************************************************
680 int CorrAxesCommand::getSharedFloat(InputData* input){
682 lookupFloat = input->getSharedRAbundFloatVectors();
683 string lastLabel = lookupFloat[0]->getLabel();
685 if (label == "") { label = lastLabel; return 0; }
687 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
688 set<string> labels; labels.insert(label);
689 set<string> processedLabels;
690 set<string> userLabels = labels;
692 //as long as you are not at the end of the file or done wih the lines you want
693 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
695 if (m->control_pressed) { return 0; }
697 if(labels.count(lookupFloat[0]->getLabel()) == 1){
698 processedLabels.insert(lookupFloat[0]->getLabel());
699 userLabels.erase(lookupFloat[0]->getLabel());
703 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
704 string saveLabel = lookupFloat[0]->getLabel();
706 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
707 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
709 processedLabels.insert(lookupFloat[0]->getLabel());
710 userLabels.erase(lookupFloat[0]->getLabel());
712 //restore real lastlabel to save below
713 lookupFloat[0]->setLabel(saveLabel);
717 lastLabel = lookupFloat[0]->getLabel();
719 //get next line to process
720 //prevent memory leak
721 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
722 lookupFloat = input->getSharedRAbundFloatVectors();
726 if (m->control_pressed) { return 0; }
728 //output error messages about any remaining user labels
729 set<string>::iterator it;
730 bool needToRun = false;
731 for (it = userLabels.begin(); it != userLabels.end(); it++) {
732 m->mothurOut("Your file does not include the label " + *it);
733 if (processedLabels.count(lastLabel) != 1) {
734 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
737 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
741 //run last label if you need to
742 if (needToRun == true) {
743 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
744 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
749 catch(exception& e) {
750 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
754 //**********************************************************************************************************************
755 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
758 vector<SharedRAbundFloatVector*> newLookup;
759 for (int i = 0; i < thislookup.size(); i++) {
760 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
761 temp->setLabel(thislookup[i]->getLabel());
762 temp->setGroup(thislookup[i]->getGroup());
763 newLookup.push_back(temp);
767 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
768 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
770 //look at each sharedRabund and make sure they are not all zero
772 for (int j = 0; j < thislookup.size(); j++) {
773 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
776 //if they are not all zero add this bin
778 for (int j = 0; j < thislookup.size(); j++) {
779 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
784 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
786 thislookup = newLookup;
791 catch(exception& e) {
792 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
796 /*****************************************************************/
797 map<string, vector<float> > CorrAxesCommand::readAxes(){
799 map<string, vector<float> > axes;
802 m->openInputFile(axesfile, in);
804 string headerLine = m->getline(in); m->gobble(in);
806 //count the number of axis you are reading
810 int pos = headerLine.find("axis");
811 if (pos != string::npos) {
813 headerLine = headerLine.substr(pos+4);
814 }else { done = true; }
817 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
821 if (m->control_pressed) { in.close(); return axes; }
824 in >> group; m->gobble(in);
826 vector<float> thisGroupsAxes;
827 for (int i = 0; i < count; i++) {
831 //only save the axis we want
832 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
835 //save group if its one the user selected
836 if (names.count(group) != 0) {
837 map<string, vector<float> >::iterator it = axes.find(group);
839 if (it == axes.end()) {
840 axes[group] = thisGroupsAxes;
842 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
852 catch(exception& e) {
853 m->errorOut(e, "CorrAxesCommand", "readAxes");
857 /*****************************************************************/
858 int CorrAxesCommand::getMetadata(){
860 vector<string> groupNames;
863 m->openInputFile(metadatafile, in);
865 string headerLine = m->getline(in); m->gobble(in);
866 istringstream iss (headerLine,istringstream::in);
868 //read the first label, because it refers to the groups
870 iss >> columnLabel; m->gobble(iss);
872 //save names of columns you are reading
874 iss >> columnLabel; m->gobble(iss);
875 metadataLabels.push_back(columnLabel);
877 int count = metadataLabels.size();
882 if (m->control_pressed) { in.close(); return 0; }
885 in >> group; m->gobble(in);
886 groupNames.push_back(group);
888 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
889 tempLookup->setGroup(group);
890 tempLookup->setLabel("1");
892 for (int i = 0; i < count; i++) {
896 tempLookup->push_back(temp, group);
899 lookupFloat.push_back(tempLookup);
905 //remove any groups the user does not want, and set globaldata->groups with only valid groups
907 util = new SharedUtil();
909 util->setGroups(globaldata->Groups, groupNames);
911 for (int i = 0; i < lookupFloat.size(); i++) {
912 //if this sharedrabund is not from a group the user wants then delete it.
913 if (util->isValidGroup(lookupFloat[i]->getGroup(), globaldata->Groups) == false) {
914 delete lookupFloat[i]; lookupFloat[i] = NULL;
915 lookupFloat.erase(lookupFloat.begin()+i);
924 catch(exception& e) {
925 m->errorOut(e, "CorrAxesCommand", "getMetadata");
929 /*****************************************************************/