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::setParameters(){
16 CommandParameter paxes("axes", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(paxes);
17 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pshared);
18 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(prelabund);
19 CommandParameter pmetadata("metadata", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pmetadata);
20 CommandParameter pnumaxes("numaxes", "Number", "", "3", "", "", "",false,false); parameters.push_back(pnumaxes);
21 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
22 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
23 CommandParameter pmethod("method", "Multiple", "pearson-spearman-kendall", "pearson", "", "", "",false,false); parameters.push_back(pmethod);
24 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
25 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
27 vector<string> myArray;
28 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
32 m->errorOut(e, "CorrAxesCommand", "setParameters");
36 //**********************************************************************************************************************
37 string CorrAxesCommand::getHelpString(){
39 string helpString = "";
40 helpString += "The corr.axes command reads a shared, relabund or metadata file as well as an axes file and calculates the correlation coefficient.\n";
41 helpString += "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";
42 helpString += "The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes.\n";
43 helpString += "The label parameter allows you to select what distance level you would like used, if none is given the first distance is used.\n";
44 helpString += "The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n";
45 helpString += "The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n";
46 helpString += "The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n";
47 helpString += "Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n";
48 helpString += "The corr.axes command outputs a .corr.axes file.\n";
49 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
53 m->errorOut(e, "CorrAxesCommand", "getHelpString");
57 //**********************************************************************************************************************
58 CorrAxesCommand::CorrAxesCommand(){
60 abort = true; calledHelp = true;
62 vector<string> tempOutNames;
63 outputTypes["corr.axes"] = tempOutNames;
66 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
70 //**********************************************************************************************************************
71 CorrAxesCommand::CorrAxesCommand(string option) {
73 abort = false; calledHelp = false;
75 //allow user to run help
76 if(option == "help") { help(); abort = true; calledHelp = true; }
77 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
80 vector<string> myArray = setParameters();
82 OptionParser parser(option);
83 map<string, string> parameters = parser.getParameters();
85 ValidParameters validParameter;
86 map<string, string>::iterator it;
88 //check to make sure all parameters are valid for command
89 for (it = parameters.begin(); it != parameters.end(); it++) {
90 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
93 vector<string> tempOutNames;
94 outputTypes["corr.axes"] = tempOutNames;
96 //if the user changes the input directory command factory will send this info to us in the output parameter
97 string inputDir = validParameter.validFile(parameters, "inputdir", false);
98 if (inputDir == "not found"){ inputDir = ""; }
101 it = parameters.find("axes");
102 //user has given a template file
103 if(it != parameters.end()){
104 path = m->hasPath(it->second);
105 //if the user has not given a path then, add inputdir. else leave path alone.
106 if (path == "") { parameters["axes"] = inputDir + it->second; }
109 it = parameters.find("shared");
110 //user has given a template file
111 if(it != parameters.end()){
112 path = m->hasPath(it->second);
113 //if the user has not given a path then, add inputdir. else leave path alone.
114 if (path == "") { parameters["shared"] = inputDir + it->second; }
117 it = parameters.find("relabund");
118 //user has given a template file
119 if(it != parameters.end()){
120 path = m->hasPath(it->second);
121 //if the user has not given a path then, add inputdir. else leave path alone.
122 if (path == "") { parameters["relabund"] = inputDir + it->second; }
125 it = parameters.find("metadata");
126 //user has given a template file
127 if(it != parameters.end()){
128 path = m->hasPath(it->second);
129 //if the user has not given a path then, add inputdir. else leave path alone.
130 if (path == "") { parameters["metadata"] = inputDir + it->second; }
135 //check for required parameters
136 axesfile = validParameter.validFile(parameters, "axes", true);
137 if (axesfile == "not open") { abort = true; }
138 else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true; }
140 sharedfile = validParameter.validFile(parameters, "shared", true);
141 if (sharedfile == "not open") { abort = true; }
142 else if (sharedfile == "not found") { sharedfile = ""; }
143 else { inputFileName = sharedfile; }
145 relabundfile = validParameter.validFile(parameters, "relabund", true);
146 if (relabundfile == "not open") { abort = true; }
147 else if (relabundfile == "not found") { relabundfile = ""; }
148 else { inputFileName = relabundfile; }
150 metadatafile = validParameter.validFile(parameters, "metadata", true);
151 if (metadatafile == "not open") { abort = true; }
152 else if (metadatafile == "not found") { metadatafile = ""; }
153 else { inputFileName = metadatafile; }
155 groups = validParameter.validFile(parameters, "groups", false);
156 if (groups == "not found") { groups = ""; pickedGroups = false; }
159 m->splitAtDash(groups, Groups);
163 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
165 label = validParameter.validFile(parameters, "label", false);
166 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=""; }
168 if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) {
169 //is there are current file available for any of these?
170 //give priority to shared, then relabund
171 //if there is a current shared file, use it
172 sharedfile = m->getSharedFile();
173 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
175 relabundfile = m->getRelAbundFile();
176 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
178 m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true;
183 if (metadatafile != "") {
184 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
186 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
189 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
190 convert(temp, numaxes);
192 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
194 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; }
198 catch(exception& e) {
199 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
203 //**********************************************************************************************************************
205 int CorrAxesCommand::execute(){
208 if (abort == true) { if (calledHelp) { return 0; } return 2; }
210 /*************************************************************************************/
211 // use smart distancing to get right sharedRabund and convert to relabund if needed //
212 /************************************************************************************/
213 if (sharedfile != "") {
214 InputData* input = new InputData(sharedfile, "sharedfile");
215 getSharedFloat(input);
218 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
219 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
221 }else if (relabundfile != "") {
222 InputData* input = new InputData(relabundfile, "relabund");
223 getSharedFloat(input);
226 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
227 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
229 }else if (metadatafile != "") {
230 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
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 metadata file."); m->mothurOutEndLine(); return 0; }
234 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
236 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
238 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
240 //this is for a sanity check to make sure the axes file and shared file match
241 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
243 /*************************************************************************************/
244 // read axes file and check for file mismatches with shared or relabund file //
245 /************************************************************************************/
248 map<string, vector<float> > axes = readAxes();
250 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
252 //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
253 if (axes.size() != lookupFloat.size()) {
254 map<string, vector<float> >::iterator it;
255 for (int i = 0; i < lookupFloat.size(); i++) {
256 it = axes.find(lookupFloat[i]->getGroup());
257 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(); }
259 m->control_pressed = true;
262 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
264 /*************************************************************************************/
265 // calc the r values //
266 /************************************************************************************/
268 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
269 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
271 m->openOutputFile(outputFileName, out);
272 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
275 if (metadatafile == "") { out << "OTU"; }
276 else { out << "Feature"; }
278 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1); }
279 out << "\tlength" << endl;
281 if (method == "pearson") { calcPearson(axes, out); }
282 else if (method == "spearman") { calcSpearman(axes, out); }
283 else if (method == "kendall") { calcKendall(axes, out); }
284 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
287 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
289 if (m->control_pressed) { return 0; }
291 m->mothurOutEndLine();
292 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
293 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
294 m->mothurOutEndLine();
298 catch(exception& e) {
299 m->errorOut(e, "CorrAxesCommand", "execute");
303 //**********************************************************************************************************************
304 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
307 //find average of each axis - X
308 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
309 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
310 vector<float> temp = it->second;
311 for (int i = 0; i < temp.size(); i++) {
312 averageAxes[i] += temp[i];
316 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
319 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
321 if (metadatafile == "") { out << i+1; }
322 else { out << metadataLabels[i]; }
324 //find the averages this otu - Y
326 for (int j = 0; j < lookupFloat.size(); j++) {
327 sumOtu += lookupFloat[j]->getAbundance(i);
329 float Ybar = sumOtu / (float) lookupFloat.size();
331 vector<float> rValues(averageAxes.size());
333 //find r value for each axis
334 for (int k = 0; k < averageAxes.size(); k++) {
337 double numerator = 0.0;
338 double denomTerm1 = 0.0;
339 double denomTerm2 = 0.0;
340 for (int j = 0; j < lookupFloat.size(); j++) {
341 float Yi = lookupFloat[j]->getAbundance(i);
342 float Xi = axes[lookupFloat[j]->getGroup()][k];
344 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
345 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
346 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
349 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
351 r = numerator / denom;
357 for(int k=0;k<rValues.size();k++){
358 sum += rValues[k] * rValues[k];
360 out << '\t' << sqrt(sum) << endl;
365 catch(exception& e) {
366 m->errorOut(e, "CorrAxesCommand", "calcPearson");
370 //**********************************************************************************************************************
371 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
375 vector< map<float, int> > tableX; tableX.resize(numaxes);
376 map<float, int>::iterator itTable;
377 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
378 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
379 vector<float> temp = it->second;
380 for (int i = 0; i < temp.size(); i++) {
381 spearmanRank member(it->first, temp[i]);
382 scores[i].push_back(member);
384 //count number of repeats
385 itTable = tableX[i].find(temp[i]);
386 if (itTable == tableX[i].end()) {
387 tableX[i][temp[i]] = 1;
389 tableX[i][temp[i]]++;
396 vector<double> Lx; Lx.resize(numaxes, 0.0);
397 for (int i = 0; i < numaxes; i++) {
398 for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) {
399 double tx = (double) itTable->second;
400 Lx[i] += ((pow(tx, 3.0) - tx) / 12.0);
405 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
407 //find ranks of xi in each axis
408 map<string, vector<float> > rankAxes;
409 for (int i = 0; i < numaxes; i++) {
411 vector<spearmanRank> ties;
413 for (int j = 0; j < scores[i].size(); j++) {
415 ties.push_back(scores[i][j]);
417 if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
418 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
420 for (int k = 0; k < ties.size(); k++) {
421 float thisrank = rankTotal / (float) ties.size();
422 rankAxes[ties[k].name].push_back(thisrank);
427 }else { // you are the last one
429 for (int k = 0; k < ties.size(); k++) {
430 float thisrank = rankTotal / (float) ties.size();
431 rankAxes[ties[k].name].push_back(thisrank);
440 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
442 if (metadatafile == "") { out << i+1; }
443 else { out << metadataLabels[i]; }
445 //find the ranks of this otu - Y
446 vector<spearmanRank> otuScores;
447 map<float, int> tableY;
448 for (int j = 0; j < lookupFloat.size(); j++) {
449 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
450 otuScores.push_back(member);
452 itTable = tableY.find(member.score);
453 if (itTable == tableY.end()) {
454 tableY[member.score] = 1;
456 tableY[member.score]++;
463 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
464 double ty = (double) itTable->second;
465 Ly += ((pow(ty, 3.0) - ty) / 12.0);
468 sort(otuScores.begin(), otuScores.end(), compareSpearman);
470 map<string, float> rankOtus;
471 vector<spearmanRank> ties;
473 for (int j = 0; j < otuScores.size(); j++) {
475 ties.push_back(otuScores[j]);
477 if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
478 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
480 for (int k = 0; k < ties.size(); k++) {
481 float thisrank = rankTotal / (float) ties.size();
482 rankOtus[ties[k].name] = thisrank;
487 }else { // you are the last one
489 for (int k = 0; k < ties.size(); k++) {
490 float thisrank = rankTotal / (float) ties.size();
491 rankOtus[ties[k].name] = thisrank;
495 vector<double> pValues(numaxes);
497 //calc spearman ranks for each axis for this otu
498 for (int j = 0; j < numaxes; j++) {
501 for (int k = 0; k < lookupFloat.size(); k++) {
503 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
504 float yi = rankOtus[lookupFloat[k]->getGroup()];
506 di += ((xi - yi) * (xi - yi));
511 double n = (double) lookupFloat.size();
513 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
514 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
516 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
524 for(int k=0;k<numaxes;k++){
525 sum += pValues[k] * pValues[k];
527 out << '\t' << sqrt(sum) << endl;
532 catch(exception& e) {
533 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
537 //**********************************************************************************************************************
538 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
542 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
543 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
544 vector<float> temp = it->second;
545 for (int i = 0; i < temp.size(); i++) {
546 spearmanRank member(it->first, temp[i]);
547 scores[i].push_back(member);
552 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
554 //convert scores to ranks of xi in each axis
555 for (int i = 0; i < numaxes; i++) {
557 vector<spearmanRank*> ties;
559 for (int j = 0; j < scores[i].size(); j++) {
561 ties.push_back(&(scores[i][j]));
563 if (j != scores[i].size()-1) { // you are not the last so you can look ahead
564 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
565 for (int k = 0; k < ties.size(); k++) {
566 float thisrank = rankTotal / (float) ties.size();
567 (*ties[k]).score = thisrank;
572 }else { // you are the last one
573 for (int k = 0; k < ties.size(); k++) {
574 float thisrank = rankTotal / (float) ties.size();
575 (*ties[k]).score = thisrank;
582 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
584 if (metadatafile == "") { out << i+1; }
585 else { out << metadataLabels[i]; }
587 //find the ranks of this otu - Y
588 vector<spearmanRank> otuScores;
589 for (int j = 0; j < lookupFloat.size(); j++) {
590 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
591 otuScores.push_back(member);
594 sort(otuScores.begin(), otuScores.end(), compareSpearman);
596 map<string, float> rankOtus;
597 vector<spearmanRank> ties;
599 for (int j = 0; j < otuScores.size(); j++) {
601 ties.push_back(otuScores[j]);
603 if (j != otuScores.size()-1) { // you are not the last so you can look ahead
604 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
605 for (int k = 0; k < ties.size(); k++) {
606 float thisrank = rankTotal / (float) ties.size();
607 rankOtus[ties[k].name] = thisrank;
612 }else { // you are the last one
613 for (int k = 0; k < ties.size(); k++) {
614 float thisrank = rankTotal / (float) ties.size();
615 rankOtus[ties[k].name] = thisrank;
621 vector<double> pValues(numaxes);
623 //calc spearman ranks for each axis for this otu
624 for (int j = 0; j < numaxes; j++) {
629 vector<spearmanRank> otus;
630 vector<spearmanRank> otusTemp;
631 for (int l = 0; l < scores[j].size(); l++) {
632 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
633 otus.push_back(member);
637 for (int l = 0; l < scores[j].size(); l++) {
639 int numWithHigherRank = 0;
640 int numWithLowerRank = 0;
641 float thisrank = otus[l].score;
643 for (int u = l+1; u < scores[j].size(); u++) {
644 if (otus[u].score > thisrank) { numWithHigherRank++; }
645 else if (otus[u].score < thisrank) { numWithLowerRank++; }
649 numCoor += numWithHigherRank;
650 numDisCoor += numWithLowerRank;
653 double p = (numCoor - numDisCoor) / (float) count;
661 for(int k=0;k<numaxes;k++){
662 sum += pValues[k] * pValues[k];
664 out << '\t' << sqrt(sum) << endl;
669 catch(exception& e) {
670 m->errorOut(e, "CorrAxesCommand", "calcKendall");
674 //**********************************************************************************************************************
675 int CorrAxesCommand::getSharedFloat(InputData* input){
677 lookupFloat = input->getSharedRAbundFloatVectors();
678 string lastLabel = lookupFloat[0]->getLabel();
680 if (label == "") { label = lastLabel; return 0; }
682 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
683 set<string> labels; labels.insert(label);
684 set<string> processedLabels;
685 set<string> userLabels = labels;
687 //as long as you are not at the end of the file or done wih the lines you want
688 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
690 if (m->control_pressed) { return 0; }
692 if(labels.count(lookupFloat[0]->getLabel()) == 1){
693 processedLabels.insert(lookupFloat[0]->getLabel());
694 userLabels.erase(lookupFloat[0]->getLabel());
698 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
699 string saveLabel = lookupFloat[0]->getLabel();
701 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
702 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
704 processedLabels.insert(lookupFloat[0]->getLabel());
705 userLabels.erase(lookupFloat[0]->getLabel());
707 //restore real lastlabel to save below
708 lookupFloat[0]->setLabel(saveLabel);
712 lastLabel = lookupFloat[0]->getLabel();
714 //get next line to process
715 //prevent memory leak
716 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
717 lookupFloat = input->getSharedRAbundFloatVectors();
721 if (m->control_pressed) { return 0; }
723 //output error messages about any remaining user labels
724 set<string>::iterator it;
725 bool needToRun = false;
726 for (it = userLabels.begin(); it != userLabels.end(); it++) {
727 m->mothurOut("Your file does not include the label " + *it);
728 if (processedLabels.count(lastLabel) != 1) {
729 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
732 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
736 //run last label if you need to
737 if (needToRun == true) {
738 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
739 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
744 catch(exception& e) {
745 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
749 //**********************************************************************************************************************
750 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
753 vector<SharedRAbundFloatVector*> newLookup;
754 for (int i = 0; i < thislookup.size(); i++) {
755 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
756 temp->setLabel(thislookup[i]->getLabel());
757 temp->setGroup(thislookup[i]->getGroup());
758 newLookup.push_back(temp);
762 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
763 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
765 //look at each sharedRabund and make sure they are not all zero
767 for (int j = 0; j < thislookup.size(); j++) {
768 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
771 //if they are not all zero add this bin
773 for (int j = 0; j < thislookup.size(); j++) {
774 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
779 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
781 thislookup = newLookup;
786 catch(exception& e) {
787 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
791 /*****************************************************************/
792 map<string, vector<float> > CorrAxesCommand::readAxes(){
794 map<string, vector<float> > axes;
797 m->openInputFile(axesfile, in);
799 string headerLine = m->getline(in); m->gobble(in);
801 //count the number of axis you are reading
805 int pos = headerLine.find("axis");
806 if (pos != string::npos) {
808 headerLine = headerLine.substr(pos+4);
809 }else { done = true; }
812 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
816 if (m->control_pressed) { in.close(); return axes; }
819 in >> group; m->gobble(in);
821 vector<float> thisGroupsAxes;
822 for (int i = 0; i < count; i++) {
826 //only save the axis we want
827 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
830 //save group if its one the user selected
831 if (names.count(group) != 0) {
832 map<string, vector<float> >::iterator it = axes.find(group);
834 if (it == axes.end()) {
835 axes[group] = thisGroupsAxes;
837 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
847 catch(exception& e) {
848 m->errorOut(e, "CorrAxesCommand", "readAxes");
852 /*****************************************************************/
853 int CorrAxesCommand::getMetadata(){
855 vector<string> groupNames;
858 m->openInputFile(metadatafile, in);
860 string headerLine = m->getline(in); m->gobble(in);
861 istringstream iss (headerLine,istringstream::in);
863 //read the first label, because it refers to the groups
865 iss >> columnLabel; m->gobble(iss);
867 //save names of columns you are reading
869 iss >> columnLabel; m->gobble(iss);
870 metadataLabels.push_back(columnLabel);
872 int count = metadataLabels.size();
877 if (m->control_pressed) { in.close(); return 0; }
880 in >> group; m->gobble(in);
881 groupNames.push_back(group);
883 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
884 tempLookup->setGroup(group);
885 tempLookup->setLabel("1");
887 for (int i = 0; i < count; i++) {
891 tempLookup->push_back(temp, group);
894 lookupFloat.push_back(tempLookup);
900 //remove any groups the user does not want, and set globaldata->groups with only valid groups
902 util = new SharedUtil();
904 util->setGroups(m->Groups, groupNames);
906 for (int i = 0; i < lookupFloat.size(); i++) {
907 //if this sharedrabund is not from a group the user wants then delete it.
908 if (util->isValidGroup(lookupFloat[i]->getGroup(), m->Groups) == false) {
909 delete lookupFloat[i]; lookupFloat[i] = NULL;
910 lookupFloat.erase(lookupFloat.begin()+i);
919 catch(exception& e) {
920 m->errorOut(e, "CorrAxesCommand", "getMetadata");
924 /*****************************************************************/