5 * Created by westcott on 12/22/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "corraxescommand.h"
11 #include "sharedutilities.h"
12 #include "linearalgebra.h"
14 //**********************************************************************************************************************
15 vector<string> CorrAxesCommand::setParameters(){
17 CommandParameter paxes("axes", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(paxes);
18 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pshared);
19 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(prelabund);
20 CommandParameter pmetadata("metadata", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pmetadata);
21 CommandParameter pnumaxes("numaxes", "Number", "", "3", "", "", "",false,false); parameters.push_back(pnumaxes);
22 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
23 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
24 CommandParameter pmethod("method", "Multiple", "pearson-spearman-kendall", "pearson", "", "", "",false,false); parameters.push_back(pmethod);
25 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
26 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
28 vector<string> myArray;
29 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
33 m->errorOut(e, "CorrAxesCommand", "setParameters");
37 //**********************************************************************************************************************
38 string CorrAxesCommand::getHelpString(){
40 string helpString = "";
41 helpString += "The corr.axes command reads a shared, relabund or metadata file as well as an axes file and calculates the correlation coefficient.\n";
42 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";
43 helpString += "The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes.\n";
44 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";
45 helpString += "The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n";
46 helpString += "The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n";
47 helpString += "The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n";
48 helpString += "Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n";
49 helpString += "The corr.axes command outputs a .corr.axes file.\n";
50 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
54 m->errorOut(e, "CorrAxesCommand", "getHelpString");
58 //**********************************************************************************************************************
59 CorrAxesCommand::CorrAxesCommand(){
61 abort = true; calledHelp = true;
63 vector<string> tempOutNames;
64 outputTypes["corr.axes"] = tempOutNames;
67 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
71 //**********************************************************************************************************************
72 CorrAxesCommand::CorrAxesCommand(string option) {
74 abort = false; calledHelp = false;
76 //allow user to run help
77 if(option == "help") { help(); abort = true; calledHelp = true; }
78 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
81 vector<string> myArray = setParameters();
83 OptionParser parser(option);
84 map<string, string> parameters = parser.getParameters();
86 ValidParameters validParameter;
87 map<string, string>::iterator it;
89 //check to make sure all parameters are valid for command
90 for (it = parameters.begin(); it != parameters.end(); it++) {
91 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
94 vector<string> tempOutNames;
95 outputTypes["corr.axes"] = tempOutNames;
97 //if the user changes the input directory command factory will send this info to us in the output parameter
98 string inputDir = validParameter.validFile(parameters, "inputdir", false);
99 if (inputDir == "not found"){ inputDir = ""; }
102 it = parameters.find("axes");
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["axes"] = inputDir + it->second; }
110 it = parameters.find("shared");
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["shared"] = inputDir + it->second; }
118 it = parameters.find("relabund");
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["relabund"] = inputDir + it->second; }
126 it = parameters.find("metadata");
127 //user has given a template file
128 if(it != parameters.end()){
129 path = m->hasPath(it->second);
130 //if the user has not given a path then, add inputdir. else leave path alone.
131 if (path == "") { parameters["metadata"] = inputDir + it->second; }
136 //check for required parameters
137 axesfile = validParameter.validFile(parameters, "axes", true);
138 if (axesfile == "not open") { abort = true; }
139 else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true; }
141 sharedfile = validParameter.validFile(parameters, "shared", true);
142 if (sharedfile == "not open") { abort = true; }
143 else if (sharedfile == "not found") { sharedfile = ""; }
144 else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
146 relabundfile = validParameter.validFile(parameters, "relabund", true);
147 if (relabundfile == "not open") { abort = true; }
148 else if (relabundfile == "not found") { relabundfile = ""; }
149 else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
151 metadatafile = validParameter.validFile(parameters, "metadata", true);
152 if (metadatafile == "not open") { abort = true; }
153 else if (metadatafile == "not found") { metadatafile = ""; }
154 else { inputFileName = metadatafile; }
156 groups = validParameter.validFile(parameters, "groups", false);
157 if (groups == "not found") { groups = ""; pickedGroups = false; }
160 m->splitAtDash(groups, Groups);
162 m->setGroups(Groups);
164 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
166 label = validParameter.validFile(parameters, "label", false);
167 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=""; }
169 if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) {
170 //is there are current file available for any of these?
171 //give priority to shared, then relabund
172 //if there is a current shared file, use it
173 sharedfile = m->getSharedFile();
174 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
176 relabundfile = m->getRelAbundFile();
177 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
179 m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true;
184 if (metadatafile != "") {
185 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
187 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
190 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
191 m->mothurConvert(temp, numaxes);
193 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
195 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; }
199 catch(exception& e) {
200 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
204 //**********************************************************************************************************************
206 int CorrAxesCommand::execute(){
209 if (abort == true) { if (calledHelp) { return 0; } return 2; }
211 /*************************************************************************************/
212 // use smart distancing to get right sharedRabund and convert to relabund if needed //
213 /************************************************************************************/
214 if (sharedfile != "") {
215 InputData* input = new InputData(sharedfile, "sharedfile");
216 getSharedFloat(input);
219 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
220 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
222 }else if (relabundfile != "") {
223 InputData* input = new InputData(relabundfile, "relabund");
224 getSharedFloat(input);
227 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
228 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
230 }else if (metadatafile != "") {
231 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
232 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
233 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
235 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
237 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
239 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
241 //this is for a sanity check to make sure the axes file and shared file match
242 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
244 /*************************************************************************************/
245 // read axes file and check for file mismatches with shared or relabund file //
246 /************************************************************************************/
249 map<string, vector<float> > axes = readAxes();
251 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
253 //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
254 if (axes.size() != lookupFloat.size()) {
255 map<string, vector<float> >::iterator it;
256 for (int i = 0; i < lookupFloat.size(); i++) {
257 it = axes.find(lookupFloat[i]->getGroup());
258 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(); }
260 m->control_pressed = true;
263 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
265 /*************************************************************************************/
266 // calc the r values //
267 /************************************************************************************/
269 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
270 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
272 m->openOutputFile(outputFileName, out);
273 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
276 if (metadatafile == "") { out << "OTU"; }
277 else { out << "Feature"; }
279 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1) << "\tp-value"; }
280 out << "\tlength" << endl;
282 if (method == "pearson") { calcPearson(axes, out); }
283 else if (method == "spearman") { calcSpearman(axes, out); }
284 else if (method == "kendall") { calcKendall(axes, out); }
285 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
288 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
290 if (m->control_pressed) { return 0; }
292 m->mothurOutEndLine();
293 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
294 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
295 m->mothurOutEndLine();
299 catch(exception& e) {
300 m->errorOut(e, "CorrAxesCommand", "execute");
304 //**********************************************************************************************************************
305 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
308 LinearAlgebra linear;
310 //find average of each axis - X
311 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
312 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
313 vector<float> temp = it->second;
314 for (int i = 0; i < temp.size(); i++) {
315 averageAxes[i] += temp[i];
319 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
322 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
324 if (metadatafile == "") { out << m->currentBinLabels[i]; }
325 else { out << metadataLabels[i]; }
327 //find the averages this otu - Y
329 for (int j = 0; j < lookupFloat.size(); j++) {
330 sumOtu += lookupFloat[j]->getAbundance(i);
332 float Ybar = sumOtu / (float) lookupFloat.size();
334 vector<float> rValues(averageAxes.size());
336 //find r value for each axis
337 for (int k = 0; k < averageAxes.size(); k++) {
340 double numerator = 0.0;
341 double denomTerm1 = 0.0;
342 double denomTerm2 = 0.0;
343 for (int j = 0; j < lookupFloat.size(); j++) {
344 float Yi = lookupFloat[j]->getAbundance(i);
345 float Xi = axes[lookupFloat[j]->getGroup()][k];
347 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
348 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
349 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
352 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
354 r = numerator / denom;
356 if (isnan(r) || isinf(r)) { r = 0.0; }
361 double sig = linear.calcPearsonSig(lookupFloat.size(), r);
367 for(int k=0;k<rValues.size();k++){
368 sum += rValues[k] * rValues[k];
370 out << '\t' << sqrt(sum) << endl;
375 catch(exception& e) {
376 m->errorOut(e, "CorrAxesCommand", "calcPearson");
380 //**********************************************************************************************************************
381 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
384 LinearAlgebra linear;
388 vector< map<float, int> > tableX; tableX.resize(numaxes);
389 map<float, int>::iterator itTable;
390 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
391 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
392 vector<float> temp = it->second;
393 for (int i = 0; i < temp.size(); i++) {
394 spearmanRank member(it->first, temp[i]);
395 scores[i].push_back(member);
397 //count number of repeats
398 itTable = tableX[i].find(temp[i]);
399 if (itTable == tableX[i].end()) {
400 tableX[i][temp[i]] = 1;
402 tableX[i][temp[i]]++;
409 vector<double> Lx; Lx.resize(numaxes, 0.0);
410 for (int i = 0; i < numaxes; i++) {
411 for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) {
412 double tx = (double) itTable->second;
413 Lx[i] += ((pow(tx, 3.0) - tx) / 12.0);
418 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
420 //find ranks of xi in each axis
421 map<string, vector<float> > rankAxes;
422 for (int i = 0; i < numaxes; i++) {
424 vector<spearmanRank> ties;
427 for (int j = 0; j < scores[i].size(); j++) {
429 ties.push_back(scores[i][j]);
431 if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
432 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
434 for (int k = 0; k < ties.size(); k++) {
435 float thisrank = rankTotal / (float) ties.size();
436 rankAxes[ties[k].name].push_back(thisrank);
443 }else { // you are the last one
445 for (int k = 0; k < ties.size(); k++) {
446 float thisrank = rankTotal / (float) ties.size();
447 rankAxes[ties[k].name].push_back(thisrank);
452 sf.push_back(sfTemp);
457 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
459 if (metadatafile == "") { out << m->currentBinLabels[i]; }
460 else { out << metadataLabels[i]; }
462 //find the ranks of this otu - Y
463 vector<spearmanRank> otuScores;
464 map<float, int> tableY;
465 for (int j = 0; j < lookupFloat.size(); j++) {
466 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
467 otuScores.push_back(member);
469 itTable = tableY.find(member.score);
470 if (itTable == tableY.end()) {
471 tableY[member.score] = 1;
473 tableY[member.score]++;
480 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
481 double ty = (double) itTable->second;
482 Ly += ((pow(ty, 3.0) - ty) / 12.0);
485 sort(otuScores.begin(), otuScores.end(), compareSpearman);
488 map<string, float> rankOtus;
489 vector<spearmanRank> ties;
491 for (int j = 0; j < otuScores.size(); j++) {
493 ties.push_back(otuScores[j]);
495 if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
496 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
498 for (int k = 0; k < ties.size(); k++) {
499 float thisrank = rankTotal / (float) ties.size();
500 rankOtus[ties[k].name] = thisrank;
507 }else { // you are the last one
509 for (int k = 0; k < ties.size(); k++) {
510 float thisrank = rankTotal / (float) ties.size();
511 rankOtus[ties[k].name] = thisrank;
515 vector<double> pValues(numaxes);
517 //calc spearman ranks for each axis for this otu
518 for (int j = 0; j < numaxes; j++) {
521 for (int k = 0; k < lookupFloat.size(); k++) {
523 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
524 float yi = rankOtus[lookupFloat[k]->getGroup()];
526 di += ((xi - yi) * (xi - yi));
531 double n = (double) lookupFloat.size();
533 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
534 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
536 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
538 if (isnan(p) || isinf(p)) { p = 0.0; }
544 double sig = linear.calcSpearmanSig(n, sf[j], sg, di);
550 for(int k=0;k<numaxes;k++){
551 sum += pValues[k] * pValues[k];
553 out << '\t' << sqrt(sum) << endl;
558 catch(exception& e) {
559 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
563 //**********************************************************************************************************************
564 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
567 LinearAlgebra linear;
570 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
571 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
572 vector<float> temp = it->second;
573 for (int i = 0; i < temp.size(); i++) {
574 spearmanRank member(it->first, temp[i]);
575 scores[i].push_back(member);
580 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
582 //convert scores to ranks of xi in each axis
583 for (int i = 0; i < numaxes; i++) {
585 vector<spearmanRank*> ties;
587 for (int j = 0; j < scores[i].size(); j++) {
589 ties.push_back(&(scores[i][j]));
591 if (j != scores[i].size()-1) { // you are not the last so you can look ahead
592 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
593 for (int k = 0; k < ties.size(); k++) {
594 float thisrank = rankTotal / (float) ties.size();
595 (*ties[k]).score = thisrank;
600 }else { // you are the last one
601 for (int k = 0; k < ties.size(); k++) {
602 float thisrank = rankTotal / (float) ties.size();
603 (*ties[k]).score = thisrank;
610 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
612 if (metadatafile == "") { out << m->currentBinLabels[i]; }
613 else { out << metadataLabels[i]; }
615 //find the ranks of this otu - Y
616 vector<spearmanRank> otuScores;
617 for (int j = 0; j < lookupFloat.size(); j++) {
618 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
619 otuScores.push_back(member);
622 sort(otuScores.begin(), otuScores.end(), compareSpearman);
624 map<string, float> rankOtus;
625 vector<spearmanRank> ties;
627 for (int j = 0; j < otuScores.size(); j++) {
629 ties.push_back(otuScores[j]);
631 if (j != otuScores.size()-1) { // you are not the last so you can look ahead
632 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
633 for (int k = 0; k < ties.size(); k++) {
634 float thisrank = rankTotal / (float) ties.size();
635 rankOtus[ties[k].name] = thisrank;
640 }else { // you are the last one
641 for (int k = 0; k < ties.size(); k++) {
642 float thisrank = rankTotal / (float) ties.size();
643 rankOtus[ties[k].name] = thisrank;
649 vector<double> pValues(numaxes);
651 //calc spearman ranks for each axis for this otu
652 for (int j = 0; j < numaxes; j++) {
657 vector<spearmanRank> otus;
658 vector<spearmanRank> otusTemp;
659 for (int l = 0; l < scores[j].size(); l++) {
660 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
661 otus.push_back(member);
665 for (int l = 0; l < scores[j].size(); l++) {
667 int numWithHigherRank = 0;
668 int numWithLowerRank = 0;
669 float thisrank = otus[l].score;
671 for (int u = l+1; u < scores[j].size(); u++) {
672 if (otus[u].score > thisrank) { numWithHigherRank++; }
673 else if (otus[u].score < thisrank) { numWithLowerRank++; }
677 numCoor += numWithHigherRank;
678 numDisCoor += numWithLowerRank;
681 double p = (numCoor - numDisCoor) / (float) count;
682 if (isnan(p) || isinf(p)) { p = 0.0; }
687 double sig = linear.calcKendallSig(scores[j].size(), p);
693 for(int k=0;k<numaxes;k++){
694 sum += pValues[k] * pValues[k];
696 out << '\t' << sqrt(sum) << endl;
701 catch(exception& e) {
702 m->errorOut(e, "CorrAxesCommand", "calcKendall");
706 //**********************************************************************************************************************
707 int CorrAxesCommand::getSharedFloat(InputData* input){
709 lookupFloat = input->getSharedRAbundFloatVectors();
710 string lastLabel = lookupFloat[0]->getLabel();
712 if (label == "") { label = lastLabel; return 0; }
714 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
715 set<string> labels; labels.insert(label);
716 set<string> processedLabels;
717 set<string> userLabels = labels;
719 //as long as you are not at the end of the file or done wih the lines you want
720 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
722 if (m->control_pressed) { return 0; }
724 if(labels.count(lookupFloat[0]->getLabel()) == 1){
725 processedLabels.insert(lookupFloat[0]->getLabel());
726 userLabels.erase(lookupFloat[0]->getLabel());
730 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
731 string saveLabel = lookupFloat[0]->getLabel();
733 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
734 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
736 processedLabels.insert(lookupFloat[0]->getLabel());
737 userLabels.erase(lookupFloat[0]->getLabel());
739 //restore real lastlabel to save below
740 lookupFloat[0]->setLabel(saveLabel);
744 lastLabel = lookupFloat[0]->getLabel();
746 //get next line to process
747 //prevent memory leak
748 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
749 lookupFloat = input->getSharedRAbundFloatVectors();
753 if (m->control_pressed) { return 0; }
755 //output error messages about any remaining user labels
756 set<string>::iterator it;
757 bool needToRun = false;
758 for (it = userLabels.begin(); it != userLabels.end(); it++) {
759 m->mothurOut("Your file does not include the label " + *it);
760 if (processedLabels.count(lastLabel) != 1) {
761 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
764 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
768 //run last label if you need to
769 if (needToRun == true) {
770 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
771 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
776 catch(exception& e) {
777 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
781 //**********************************************************************************************************************
782 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
785 vector<SharedRAbundFloatVector*> newLookup;
786 for (int i = 0; i < thislookup.size(); i++) {
787 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
788 temp->setLabel(thislookup[i]->getLabel());
789 temp->setGroup(thislookup[i]->getGroup());
790 newLookup.push_back(temp);
794 vector<string> newBinLabels;
795 string snumBins = toString(thislookup[0]->getNumBins());
796 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
797 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
799 //look at each sharedRabund and make sure they are not all zero
801 for (int j = 0; j < thislookup.size(); j++) {
802 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
805 //if they are not all zero add this bin
807 for (int j = 0; j < thislookup.size(); j++) {
808 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
811 //if there is a bin label use it otherwise make one
812 string binLabel = "Otu";
813 string sbinNumber = toString(i+1);
814 if (sbinNumber.length() < snumBins.length()) {
815 int diff = snumBins.length() - sbinNumber.length();
816 for (int h = 0; h < diff; h++) { binLabel += "0"; }
818 binLabel += sbinNumber;
819 if (i < m->currentBinLabels.size()) { binLabel = m->currentBinLabels[i]; }
821 newBinLabels.push_back(binLabel);
825 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
827 thislookup = newLookup;
828 m->currentBinLabels = newBinLabels;
833 catch(exception& e) {
834 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
838 /*****************************************************************/
839 map<string, vector<float> > CorrAxesCommand::readAxes(){
841 map<string, vector<float> > axes;
844 m->openInputFile(axesfile, in);
846 string headerLine = m->getline(in); m->gobble(in);
848 //count the number of axis you are reading
852 int pos = headerLine.find("axis");
853 if (pos != string::npos) {
855 headerLine = headerLine.substr(pos+4);
856 }else { done = true; }
859 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
863 if (m->control_pressed) { in.close(); return axes; }
866 in >> group; m->gobble(in);
868 vector<float> thisGroupsAxes;
869 for (int i = 0; i < count; i++) {
873 //only save the axis we want
874 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
877 //save group if its one the user selected
878 if (names.count(group) != 0) {
879 map<string, vector<float> >::iterator it = axes.find(group);
881 if (it == axes.end()) {
882 axes[group] = thisGroupsAxes;
884 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
894 catch(exception& e) {
895 m->errorOut(e, "CorrAxesCommand", "readAxes");
899 /*****************************************************************/
900 int CorrAxesCommand::getMetadata(){
902 vector<string> groupNames;
905 m->openInputFile(metadatafile, in);
907 string headerLine = m->getline(in); m->gobble(in);
908 istringstream iss (headerLine,istringstream::in);
910 //read the first label, because it refers to the groups
912 iss >> columnLabel; m->gobble(iss);
914 //save names of columns you are reading
916 iss >> columnLabel; m->gobble(iss);
917 metadataLabels.push_back(columnLabel);
919 int count = metadataLabels.size();
924 if (m->control_pressed) { in.close(); return 0; }
927 in >> group; m->gobble(in);
928 groupNames.push_back(group);
930 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
931 tempLookup->setGroup(group);
932 tempLookup->setLabel("1");
934 for (int i = 0; i < count; i++) {
937 tempLookup->push_back(temp, group);
940 lookupFloat.push_back(tempLookup);
946 //remove any groups the user does not want, and set globaldata->groups with only valid groups
948 util = new SharedUtil();
949 Groups = m->getGroups();
950 util->setGroups(Groups, groupNames);
951 m->setGroups(Groups);
953 for (int i = 0; i < lookupFloat.size(); i++) {
954 //if this sharedrabund is not from a group the user wants then delete it.
955 if (util->isValidGroup(lookupFloat[i]->getGroup(), m->getGroups()) == false) {
956 delete lookupFloat[i]; lookupFloat[i] = NULL;
957 lookupFloat.erase(lookupFloat.begin()+i);
966 catch(exception& e) {
967 m->errorOut(e, "CorrAxesCommand", "getMetadata");
971 /*****************************************************************/