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; }
79 vector<string> myArray = setParameters();
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);
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 == "")) {
168 //is there are current file available for any of these?
169 //give priority to shared, then relabund
170 //if there is a current shared file, use it
171 sharedfile = m->getSharedFile();
172 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
174 relabundfile = m->getRelAbundFile();
175 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
177 m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true;
182 if (metadatafile != "") {
183 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
185 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
188 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
189 convert(temp, numaxes);
191 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
193 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; }
197 catch(exception& e) {
198 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
202 //**********************************************************************************************************************
204 int CorrAxesCommand::execute(){
207 if (abort == true) { if (calledHelp) { return 0; } return 2; }
209 /*************************************************************************************/
210 // use smart distancing to get right sharedRabund and convert to relabund if needed //
211 /************************************************************************************/
212 if (sharedfile != "") {
213 InputData* input = new InputData(sharedfile, "sharedfile");
214 getSharedFloat(input);
217 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
218 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
220 }else if (relabundfile != "") {
221 InputData* input = new InputData(relabundfile, "relabund");
222 getSharedFloat(input);
225 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
226 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
228 }else if (metadatafile != "") {
229 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
230 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
231 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
233 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
235 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
237 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
239 //this is for a sanity check to make sure the axes file and shared file match
240 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
242 /*************************************************************************************/
243 // read axes file and check for file mismatches with shared or relabund file //
244 /************************************************************************************/
247 map<string, vector<float> > axes = readAxes();
249 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
251 //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
252 if (axes.size() != lookupFloat.size()) {
253 map<string, vector<float> >::iterator it;
254 for (int i = 0; i < lookupFloat.size(); i++) {
255 it = axes.find(lookupFloat[i]->getGroup());
256 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(); }
258 m->control_pressed = true;
261 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
263 /*************************************************************************************/
264 // calc the r values //
265 /************************************************************************************/
267 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
268 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
270 m->openOutputFile(outputFileName, out);
271 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
274 if (metadatafile == "") { out << "OTU"; }
275 else { out << "Feature"; }
277 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1); }
278 out << "\tlength" << endl;
280 if (method == "pearson") { calcPearson(axes, out); }
281 else if (method == "spearman") { calcSpearman(axes, out); }
282 else if (method == "kendall") { calcKendall(axes, out); }
283 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
286 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
288 if (m->control_pressed) { return 0; }
290 m->mothurOutEndLine();
291 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
292 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
293 m->mothurOutEndLine();
297 catch(exception& e) {
298 m->errorOut(e, "CorrAxesCommand", "execute");
302 //**********************************************************************************************************************
303 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
306 //find average of each axis - X
307 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
308 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
309 vector<float> temp = it->second;
310 for (int i = 0; i < temp.size(); i++) {
311 averageAxes[i] += temp[i];
315 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
318 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
320 if (metadatafile == "") { out << i+1; }
321 else { out << metadataLabels[i]; }
323 //find the averages this otu - Y
325 for (int j = 0; j < lookupFloat.size(); j++) {
326 sumOtu += lookupFloat[j]->getAbundance(i);
328 float Ybar = sumOtu / (float) lookupFloat.size();
330 vector<float> rValues(averageAxes.size());
332 //find r value for each axis
333 for (int k = 0; k < averageAxes.size(); k++) {
336 double numerator = 0.0;
337 double denomTerm1 = 0.0;
338 double denomTerm2 = 0.0;
339 for (int j = 0; j < lookupFloat.size(); j++) {
340 float Yi = lookupFloat[j]->getAbundance(i);
341 float Xi = axes[lookupFloat[j]->getGroup()][k];
343 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
344 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
345 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
348 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
350 r = numerator / denom;
356 for(int k=0;k<rValues.size();k++){
357 sum += rValues[k] * rValues[k];
359 out << '\t' << sqrt(sum) << endl;
364 catch(exception& e) {
365 m->errorOut(e, "CorrAxesCommand", "calcPearson");
369 //**********************************************************************************************************************
370 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
374 vector< map<float, int> > tableX; tableX.resize(numaxes);
375 map<float, int>::iterator itTable;
376 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
377 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
378 vector<float> temp = it->second;
379 for (int i = 0; i < temp.size(); i++) {
380 spearmanRank member(it->first, temp[i]);
381 scores[i].push_back(member);
383 //count number of repeats
384 itTable = tableX[i].find(temp[i]);
385 if (itTable == tableX[i].end()) {
386 tableX[i][temp[i]] = 1;
388 tableX[i][temp[i]]++;
395 vector<double> Lx; Lx.resize(numaxes, 0.0);
396 for (int i = 0; i < numaxes; i++) {
397 for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) {
398 double tx = (double) itTable->second;
399 Lx[i] += ((pow(tx, 3.0) - tx) / 12.0);
404 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
406 //find ranks of xi in each axis
407 map<string, vector<float> > rankAxes;
408 for (int i = 0; i < numaxes; i++) {
410 vector<spearmanRank> ties;
412 for (int j = 0; j < scores[i].size(); j++) {
414 ties.push_back(scores[i][j]);
416 if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
417 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
419 for (int k = 0; k < ties.size(); k++) {
420 float thisrank = rankTotal / (float) ties.size();
421 rankAxes[ties[k].name].push_back(thisrank);
426 }else { // you are the last one
428 for (int k = 0; k < ties.size(); k++) {
429 float thisrank = rankTotal / (float) ties.size();
430 rankAxes[ties[k].name].push_back(thisrank);
439 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
441 if (metadatafile == "") { out << i+1; }
442 else { out << metadataLabels[i]; }
444 //find the ranks of this otu - Y
445 vector<spearmanRank> otuScores;
446 map<float, int> tableY;
447 for (int j = 0; j < lookupFloat.size(); j++) {
448 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
449 otuScores.push_back(member);
451 itTable = tableY.find(member.score);
452 if (itTable == tableY.end()) {
453 tableY[member.score] = 1;
455 tableY[member.score]++;
462 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
463 double ty = (double) itTable->second;
464 Ly += ((pow(ty, 3.0) - ty) / 12.0);
467 sort(otuScores.begin(), otuScores.end(), compareSpearman);
469 map<string, float> rankOtus;
470 vector<spearmanRank> ties;
472 for (int j = 0; j < otuScores.size(); j++) {
474 ties.push_back(otuScores[j]);
476 if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
477 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
479 for (int k = 0; k < ties.size(); k++) {
480 float thisrank = rankTotal / (float) ties.size();
481 rankOtus[ties[k].name] = thisrank;
486 }else { // you are the last one
488 for (int k = 0; k < ties.size(); k++) {
489 float thisrank = rankTotal / (float) ties.size();
490 rankOtus[ties[k].name] = thisrank;
494 vector<double> pValues(numaxes);
496 //calc spearman ranks for each axis for this otu
497 for (int j = 0; j < numaxes; j++) {
500 for (int k = 0; k < lookupFloat.size(); k++) {
502 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
503 float yi = rankOtus[lookupFloat[k]->getGroup()];
505 di += ((xi - yi) * (xi - yi));
510 double n = (double) lookupFloat.size();
512 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
513 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
515 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
523 for(int k=0;k<numaxes;k++){
524 sum += pValues[k] * pValues[k];
526 out << '\t' << sqrt(sum) << endl;
531 catch(exception& e) {
532 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
536 //**********************************************************************************************************************
537 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
541 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
542 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
543 vector<float> temp = it->second;
544 for (int i = 0; i < temp.size(); i++) {
545 spearmanRank member(it->first, temp[i]);
546 scores[i].push_back(member);
551 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
553 //convert scores to ranks of xi in each axis
554 for (int i = 0; i < numaxes; i++) {
556 vector<spearmanRank*> ties;
558 for (int j = 0; j < scores[i].size(); j++) {
560 ties.push_back(&(scores[i][j]));
562 if (j != scores[i].size()-1) { // you are not the last so you can look ahead
563 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
564 for (int k = 0; k < ties.size(); k++) {
565 float thisrank = rankTotal / (float) ties.size();
566 (*ties[k]).score = thisrank;
571 }else { // you are the last one
572 for (int k = 0; k < ties.size(); k++) {
573 float thisrank = rankTotal / (float) ties.size();
574 (*ties[k]).score = thisrank;
581 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
583 if (metadatafile == "") { out << i+1; }
584 else { out << metadataLabels[i]; }
586 //find the ranks of this otu - Y
587 vector<spearmanRank> otuScores;
588 for (int j = 0; j < lookupFloat.size(); j++) {
589 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
590 otuScores.push_back(member);
593 sort(otuScores.begin(), otuScores.end(), compareSpearman);
595 map<string, float> rankOtus;
596 vector<spearmanRank> ties;
598 for (int j = 0; j < otuScores.size(); j++) {
600 ties.push_back(otuScores[j]);
602 if (j != otuScores.size()-1) { // you are not the last so you can look ahead
603 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
604 for (int k = 0; k < ties.size(); k++) {
605 float thisrank = rankTotal / (float) ties.size();
606 rankOtus[ties[k].name] = thisrank;
611 }else { // you are the last one
612 for (int k = 0; k < ties.size(); k++) {
613 float thisrank = rankTotal / (float) ties.size();
614 rankOtus[ties[k].name] = thisrank;
620 vector<double> pValues(numaxes);
622 //calc spearman ranks for each axis for this otu
623 for (int j = 0; j < numaxes; j++) {
628 vector<spearmanRank> otus;
629 vector<spearmanRank> otusTemp;
630 for (int l = 0; l < scores[j].size(); l++) {
631 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
632 otus.push_back(member);
636 for (int l = 0; l < scores[j].size(); l++) {
638 int numWithHigherRank = 0;
639 int numWithLowerRank = 0;
640 float thisrank = otus[l].score;
642 for (int u = l+1; u < scores[j].size(); u++) {
643 if (otus[u].score > thisrank) { numWithHigherRank++; }
644 else if (otus[u].score < thisrank) { numWithLowerRank++; }
648 numCoor += numWithHigherRank;
649 numDisCoor += numWithLowerRank;
652 double p = (numCoor - numDisCoor) / (float) count;
660 for(int k=0;k<numaxes;k++){
661 sum += pValues[k] * pValues[k];
663 out << '\t' << sqrt(sum) << endl;
668 catch(exception& e) {
669 m->errorOut(e, "CorrAxesCommand", "calcKendall");
673 //**********************************************************************************************************************
674 int CorrAxesCommand::getSharedFloat(InputData* input){
676 lookupFloat = input->getSharedRAbundFloatVectors();
677 string lastLabel = lookupFloat[0]->getLabel();
679 if (label == "") { label = lastLabel; return 0; }
681 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
682 set<string> labels; labels.insert(label);
683 set<string> processedLabels;
684 set<string> userLabels = labels;
686 //as long as you are not at the end of the file or done wih the lines you want
687 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
689 if (m->control_pressed) { return 0; }
691 if(labels.count(lookupFloat[0]->getLabel()) == 1){
692 processedLabels.insert(lookupFloat[0]->getLabel());
693 userLabels.erase(lookupFloat[0]->getLabel());
697 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
698 string saveLabel = lookupFloat[0]->getLabel();
700 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
701 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
703 processedLabels.insert(lookupFloat[0]->getLabel());
704 userLabels.erase(lookupFloat[0]->getLabel());
706 //restore real lastlabel to save below
707 lookupFloat[0]->setLabel(saveLabel);
711 lastLabel = lookupFloat[0]->getLabel();
713 //get next line to process
714 //prevent memory leak
715 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
716 lookupFloat = input->getSharedRAbundFloatVectors();
720 if (m->control_pressed) { return 0; }
722 //output error messages about any remaining user labels
723 set<string>::iterator it;
724 bool needToRun = false;
725 for (it = userLabels.begin(); it != userLabels.end(); it++) {
726 m->mothurOut("Your file does not include the label " + *it);
727 if (processedLabels.count(lastLabel) != 1) {
728 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
731 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
735 //run last label if you need to
736 if (needToRun == true) {
737 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
738 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
743 catch(exception& e) {
744 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
748 //**********************************************************************************************************************
749 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
752 vector<SharedRAbundFloatVector*> newLookup;
753 for (int i = 0; i < thislookup.size(); i++) {
754 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
755 temp->setLabel(thislookup[i]->getLabel());
756 temp->setGroup(thislookup[i]->getGroup());
757 newLookup.push_back(temp);
761 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
762 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
764 //look at each sharedRabund and make sure they are not all zero
766 for (int j = 0; j < thislookup.size(); j++) {
767 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
770 //if they are not all zero add this bin
772 for (int j = 0; j < thislookup.size(); j++) {
773 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
778 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
780 thislookup = newLookup;
785 catch(exception& e) {
786 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
790 /*****************************************************************/
791 map<string, vector<float> > CorrAxesCommand::readAxes(){
793 map<string, vector<float> > axes;
796 m->openInputFile(axesfile, in);
798 string headerLine = m->getline(in); m->gobble(in);
800 //count the number of axis you are reading
804 int pos = headerLine.find("axis");
805 if (pos != string::npos) {
807 headerLine = headerLine.substr(pos+4);
808 }else { done = true; }
811 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
815 if (m->control_pressed) { in.close(); return axes; }
818 in >> group; m->gobble(in);
820 vector<float> thisGroupsAxes;
821 for (int i = 0; i < count; i++) {
825 //only save the axis we want
826 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
829 //save group if its one the user selected
830 if (names.count(group) != 0) {
831 map<string, vector<float> >::iterator it = axes.find(group);
833 if (it == axes.end()) {
834 axes[group] = thisGroupsAxes;
836 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
846 catch(exception& e) {
847 m->errorOut(e, "CorrAxesCommand", "readAxes");
851 /*****************************************************************/
852 int CorrAxesCommand::getMetadata(){
854 vector<string> groupNames;
857 m->openInputFile(metadatafile, in);
859 string headerLine = m->getline(in); m->gobble(in);
860 istringstream iss (headerLine,istringstream::in);
862 //read the first label, because it refers to the groups
864 iss >> columnLabel; m->gobble(iss);
866 //save names of columns you are reading
868 iss >> columnLabel; m->gobble(iss);
869 metadataLabels.push_back(columnLabel);
871 int count = metadataLabels.size();
876 if (m->control_pressed) { in.close(); return 0; }
879 in >> group; m->gobble(in);
880 groupNames.push_back(group);
882 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
883 tempLookup->setGroup(group);
884 tempLookup->setLabel("1");
886 for (int i = 0; i < count; i++) {
890 tempLookup->push_back(temp, group);
893 lookupFloat.push_back(tempLookup);
899 //remove any groups the user does not want, and set globaldata->groups with only valid groups
901 util = new SharedUtil();
903 util->setGroups(m->Groups, groupNames);
905 for (int i = 0; i < lookupFloat.size(); i++) {
906 //if this sharedrabund is not from a group the user wants then delete it.
907 if (util->isValidGroup(lookupFloat[i]->getGroup(), m->Groups) == false) {
908 delete lookupFloat[i]; lookupFloat[i] = NULL;
909 lookupFloat.erase(lookupFloat.begin()+i);
918 catch(exception& e) {
919 m->errorOut(e, "CorrAxesCommand", "getMetadata");
923 /*****************************************************************/