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 string CorrAxesCommand::getOutputFileNameTag(string type, string inputName=""){
61 string outputFileName = "";
62 map<string, vector<string> >::iterator it;
64 //is this a type this command creates
65 it = outputTypes.find(type);
66 if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
68 if (type == "corraxes") { outputFileName = "corr.axes"; }
69 else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
71 return outputFileName;
74 m->errorOut(e, "CorrAxesCommand", "getOutputFileNameTag");
78 //**********************************************************************************************************************
79 CorrAxesCommand::CorrAxesCommand(){
81 abort = true; calledHelp = true;
83 vector<string> tempOutNames;
84 outputTypes["corraxes"] = tempOutNames;
87 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
91 //**********************************************************************************************************************
92 CorrAxesCommand::CorrAxesCommand(string option) {
94 abort = false; calledHelp = false;
96 //allow user to run help
97 if(option == "help") { help(); abort = true; calledHelp = true; }
98 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
101 vector<string> myArray = setParameters();
103 OptionParser parser(option);
104 map<string, string> parameters = parser.getParameters();
106 ValidParameters validParameter;
107 map<string, string>::iterator it;
109 //check to make sure all parameters are valid for command
110 for (it = parameters.begin(); it != parameters.end(); it++) {
111 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
114 vector<string> tempOutNames;
115 outputTypes["corraxes"] = tempOutNames;
117 //if the user changes the input directory command factory will send this info to us in the output parameter
118 string inputDir = validParameter.validFile(parameters, "inputdir", false);
119 if (inputDir == "not found"){ inputDir = ""; }
122 it = parameters.find("axes");
123 //user has given a template file
124 if(it != parameters.end()){
125 path = m->hasPath(it->second);
126 //if the user has not given a path then, add inputdir. else leave path alone.
127 if (path == "") { parameters["axes"] = inputDir + it->second; }
130 it = parameters.find("shared");
131 //user has given a template file
132 if(it != parameters.end()){
133 path = m->hasPath(it->second);
134 //if the user has not given a path then, add inputdir. else leave path alone.
135 if (path == "") { parameters["shared"] = inputDir + it->second; }
138 it = parameters.find("relabund");
139 //user has given a template file
140 if(it != parameters.end()){
141 path = m->hasPath(it->second);
142 //if the user has not given a path then, add inputdir. else leave path alone.
143 if (path == "") { parameters["relabund"] = inputDir + it->second; }
146 it = parameters.find("metadata");
147 //user has given a template file
148 if(it != parameters.end()){
149 path = m->hasPath(it->second);
150 //if the user has not given a path then, add inputdir. else leave path alone.
151 if (path == "") { parameters["metadata"] = inputDir + it->second; }
156 //check for required parameters
157 axesfile = validParameter.validFile(parameters, "axes", true);
158 if (axesfile == "not open") { abort = true; }
159 else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true; }
161 sharedfile = validParameter.validFile(parameters, "shared", true);
162 if (sharedfile == "not open") { abort = true; }
163 else if (sharedfile == "not found") { sharedfile = ""; }
164 else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
166 relabundfile = validParameter.validFile(parameters, "relabund", true);
167 if (relabundfile == "not open") { abort = true; }
168 else if (relabundfile == "not found") { relabundfile = ""; }
169 else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
171 metadatafile = validParameter.validFile(parameters, "metadata", true);
172 if (metadatafile == "not open") { abort = true; }
173 else if (metadatafile == "not found") { metadatafile = ""; }
174 else { inputFileName = metadatafile; }
176 groups = validParameter.validFile(parameters, "groups", false);
177 if (groups == "not found") { groups = ""; pickedGroups = false; }
180 m->splitAtDash(groups, Groups);
182 m->setGroups(Groups);
184 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
186 label = validParameter.validFile(parameters, "label", false);
187 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=""; }
189 if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) {
190 //is there are current file available for any of these?
191 //give priority to shared, then relabund
192 //if there is a current shared file, use it
193 sharedfile = m->getSharedFile();
194 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
196 relabundfile = m->getRelAbundFile();
197 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
199 m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true;
204 if (metadatafile != "") {
205 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
207 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
210 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
211 m->mothurConvert(temp, numaxes);
213 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
215 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; }
219 catch(exception& e) {
220 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
224 //**********************************************************************************************************************
226 int CorrAxesCommand::execute(){
229 if (abort == true) { if (calledHelp) { return 0; } return 2; }
231 /*************************************************************************************/
232 // use smart distancing to get right sharedRabund and convert to relabund if needed //
233 /************************************************************************************/
234 if (sharedfile != "") {
235 InputData* input = new InputData(sharedfile, "sharedfile");
236 getSharedFloat(input);
239 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
240 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
242 }else if (relabundfile != "") {
243 InputData* input = new InputData(relabundfile, "relabund");
244 getSharedFloat(input);
247 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
248 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
250 }else if (metadatafile != "") {
251 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
252 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
253 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
255 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
257 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
259 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
261 //this is for a sanity check to make sure the axes file and shared file match
262 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
264 /*************************************************************************************/
265 // read axes file and check for file mismatches with shared or relabund file //
266 /************************************************************************************/
269 map<string, vector<float> > axes = readAxes();
271 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
273 //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
274 if (axes.size() != lookupFloat.size()) {
275 map<string, vector<float> >::iterator it;
276 for (int i = 0; i < lookupFloat.size(); i++) {
277 it = axes.find(lookupFloat[i]->getGroup());
278 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(); }
280 m->control_pressed = true;
283 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
285 /*************************************************************************************/
286 // calc the r values //
287 /************************************************************************************/
289 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + "." + getOutputFileNameTag("corraxes");
290 outputNames.push_back(outputFileName); outputTypes["corraxes"].push_back(outputFileName);
292 m->openOutputFile(outputFileName, out);
293 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
296 if (metadatafile == "") { out << "OTU"; }
297 else { out << "Feature"; }
299 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1) << "\tp-value"; }
300 out << "\tlength" << endl;
302 if (method == "pearson") { calcPearson(axes, out); }
303 else if (method == "spearman") { calcSpearman(axes, out); }
304 else if (method == "kendall") { calcKendall(axes, out); }
305 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
308 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
310 if (m->control_pressed) { return 0; }
312 m->mothurOutEndLine();
313 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
314 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
315 m->mothurOutEndLine();
319 catch(exception& e) {
320 m->errorOut(e, "CorrAxesCommand", "execute");
324 //**********************************************************************************************************************
325 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
328 LinearAlgebra linear;
330 //find average of each axis - X
331 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
332 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
333 vector<float> temp = it->second;
334 for (int i = 0; i < temp.size(); i++) {
335 averageAxes[i] += temp[i];
339 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
342 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
344 if (metadatafile == "") { out << m->currentBinLabels[i]; }
345 else { out << metadataLabels[i]; }
347 //find the averages this otu - Y
349 for (int j = 0; j < lookupFloat.size(); j++) {
350 sumOtu += lookupFloat[j]->getAbundance(i);
352 float Ybar = sumOtu / (float) lookupFloat.size();
354 vector<float> rValues(averageAxes.size());
356 //find r value for each axis
357 for (int k = 0; k < averageAxes.size(); k++) {
360 double numerator = 0.0;
361 double denomTerm1 = 0.0;
362 double denomTerm2 = 0.0;
363 for (int j = 0; j < lookupFloat.size(); j++) {
364 float Yi = lookupFloat[j]->getAbundance(i);
365 float Xi = axes[lookupFloat[j]->getGroup()][k];
367 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
368 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
369 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
372 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
374 r = numerator / denom;
376 if (isnan(r) || isinf(r)) { r = 0.0; }
381 double sig = linear.calcPearsonSig(lookupFloat.size(), r);
387 for(int k=0;k<rValues.size();k++){
388 sum += rValues[k] * rValues[k];
390 out << '\t' << sqrt(sum) << endl;
395 catch(exception& e) {
396 m->errorOut(e, "CorrAxesCommand", "calcPearson");
400 //**********************************************************************************************************************
401 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
404 LinearAlgebra linear;
408 vector< map<float, int> > tableX; tableX.resize(numaxes);
409 map<float, int>::iterator itTable;
410 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
411 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
412 vector<float> temp = it->second;
413 for (int i = 0; i < temp.size(); i++) {
414 spearmanRank member(it->first, temp[i]);
415 scores[i].push_back(member);
417 //count number of repeats
418 itTable = tableX[i].find(temp[i]);
419 if (itTable == tableX[i].end()) {
420 tableX[i][temp[i]] = 1;
422 tableX[i][temp[i]]++;
429 vector<double> Lx; Lx.resize(numaxes, 0.0);
430 for (int i = 0; i < numaxes; i++) {
431 for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) {
432 double tx = (double) itTable->second;
433 Lx[i] += ((pow(tx, 3.0) - tx) / 12.0);
438 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
440 //find ranks of xi in each axis
441 map<string, vector<float> > rankAxes;
442 for (int i = 0; i < numaxes; i++) {
444 vector<spearmanRank> ties;
447 for (int j = 0; j < scores[i].size(); j++) {
449 ties.push_back(scores[i][j]);
451 if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
452 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
454 for (int k = 0; k < ties.size(); k++) {
455 float thisrank = rankTotal / (float) ties.size();
456 rankAxes[ties[k].name].push_back(thisrank);
463 }else { // you are the last one
465 for (int k = 0; k < ties.size(); k++) {
466 float thisrank = rankTotal / (float) ties.size();
467 rankAxes[ties[k].name].push_back(thisrank);
472 sf.push_back(sfTemp);
477 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
479 if (metadatafile == "") { out << m->currentBinLabels[i]; }
480 else { out << metadataLabels[i]; }
482 //find the ranks of this otu - Y
483 vector<spearmanRank> otuScores;
484 map<float, int> tableY;
485 for (int j = 0; j < lookupFloat.size(); j++) {
486 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
487 otuScores.push_back(member);
489 itTable = tableY.find(member.score);
490 if (itTable == tableY.end()) {
491 tableY[member.score] = 1;
493 tableY[member.score]++;
500 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
501 double ty = (double) itTable->second;
502 Ly += ((pow(ty, 3.0) - ty) / 12.0);
505 sort(otuScores.begin(), otuScores.end(), compareSpearman);
508 map<string, float> rankOtus;
509 vector<spearmanRank> ties;
511 for (int j = 0; j < otuScores.size(); j++) {
513 ties.push_back(otuScores[j]);
515 if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
516 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
518 for (int k = 0; k < ties.size(); k++) {
519 float thisrank = rankTotal / (float) ties.size();
520 rankOtus[ties[k].name] = thisrank;
527 }else { // you are the last one
529 for (int k = 0; k < ties.size(); k++) {
530 float thisrank = rankTotal / (float) ties.size();
531 rankOtus[ties[k].name] = thisrank;
535 vector<double> pValues(numaxes);
537 //calc spearman ranks for each axis for this otu
538 for (int j = 0; j < numaxes; j++) {
541 for (int k = 0; k < lookupFloat.size(); k++) {
543 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
544 float yi = rankOtus[lookupFloat[k]->getGroup()];
546 di += ((xi - yi) * (xi - yi));
551 double n = (double) lookupFloat.size();
553 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
554 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
556 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
558 if (isnan(p) || isinf(p)) { p = 0.0; }
564 double sig = linear.calcSpearmanSig(n, sf[j], sg, di);
570 for(int k=0;k<numaxes;k++){
571 sum += pValues[k] * pValues[k];
573 out << '\t' << sqrt(sum) << endl;
578 catch(exception& e) {
579 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
583 //**********************************************************************************************************************
584 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
587 LinearAlgebra linear;
590 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
591 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
592 vector<float> temp = it->second;
593 for (int i = 0; i < temp.size(); i++) {
594 spearmanRank member(it->first, temp[i]);
595 scores[i].push_back(member);
600 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
602 //convert scores to ranks of xi in each axis
603 for (int i = 0; i < numaxes; i++) {
605 vector<spearmanRank*> ties;
607 for (int j = 0; j < scores[i].size(); j++) {
609 ties.push_back(&(scores[i][j]));
611 if (j != scores[i].size()-1) { // you are not the last so you can look ahead
612 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
613 for (int k = 0; k < ties.size(); k++) {
614 float thisrank = rankTotal / (float) ties.size();
615 (*ties[k]).score = thisrank;
620 }else { // you are the last one
621 for (int k = 0; k < ties.size(); k++) {
622 float thisrank = rankTotal / (float) ties.size();
623 (*ties[k]).score = thisrank;
630 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
632 if (metadatafile == "") { out << m->currentBinLabels[i]; }
633 else { out << metadataLabels[i]; }
635 //find the ranks of this otu - Y
636 vector<spearmanRank> otuScores;
637 for (int j = 0; j < lookupFloat.size(); j++) {
638 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
639 otuScores.push_back(member);
642 sort(otuScores.begin(), otuScores.end(), compareSpearman);
644 map<string, float> rankOtus;
645 vector<spearmanRank> ties;
647 for (int j = 0; j < otuScores.size(); j++) {
649 ties.push_back(otuScores[j]);
651 if (j != otuScores.size()-1) { // you are not the last so you can look ahead
652 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
653 for (int k = 0; k < ties.size(); k++) {
654 float thisrank = rankTotal / (float) ties.size();
655 rankOtus[ties[k].name] = thisrank;
660 }else { // you are the last one
661 for (int k = 0; k < ties.size(); k++) {
662 float thisrank = rankTotal / (float) ties.size();
663 rankOtus[ties[k].name] = thisrank;
669 vector<double> pValues(numaxes);
671 //calc spearman ranks for each axis for this otu
672 for (int j = 0; j < numaxes; j++) {
677 vector<spearmanRank> otus;
678 vector<spearmanRank> otusTemp;
679 for (int l = 0; l < scores[j].size(); l++) {
680 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
681 otus.push_back(member);
685 for (int l = 0; l < scores[j].size(); l++) {
687 int numWithHigherRank = 0;
688 int numWithLowerRank = 0;
689 float thisrank = otus[l].score;
691 for (int u = l+1; u < scores[j].size(); u++) {
692 if (otus[u].score > thisrank) { numWithHigherRank++; }
693 else if (otus[u].score < thisrank) { numWithLowerRank++; }
697 numCoor += numWithHigherRank;
698 numDisCoor += numWithLowerRank;
701 double p = (numCoor - numDisCoor) / (float) count;
702 if (isnan(p) || isinf(p)) { p = 0.0; }
707 double sig = linear.calcKendallSig(scores[j].size(), p);
713 for(int k=0;k<numaxes;k++){
714 sum += pValues[k] * pValues[k];
716 out << '\t' << sqrt(sum) << endl;
721 catch(exception& e) {
722 m->errorOut(e, "CorrAxesCommand", "calcKendall");
726 //**********************************************************************************************************************
727 int CorrAxesCommand::getSharedFloat(InputData* input){
729 lookupFloat = input->getSharedRAbundFloatVectors();
730 string lastLabel = lookupFloat[0]->getLabel();
732 if (label == "") { label = lastLabel; return 0; }
734 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
735 set<string> labels; labels.insert(label);
736 set<string> processedLabels;
737 set<string> userLabels = labels;
739 //as long as you are not at the end of the file or done wih the lines you want
740 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
742 if (m->control_pressed) { return 0; }
744 if(labels.count(lookupFloat[0]->getLabel()) == 1){
745 processedLabels.insert(lookupFloat[0]->getLabel());
746 userLabels.erase(lookupFloat[0]->getLabel());
750 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
751 string saveLabel = lookupFloat[0]->getLabel();
753 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
754 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
756 processedLabels.insert(lookupFloat[0]->getLabel());
757 userLabels.erase(lookupFloat[0]->getLabel());
759 //restore real lastlabel to save below
760 lookupFloat[0]->setLabel(saveLabel);
764 lastLabel = lookupFloat[0]->getLabel();
766 //get next line to process
767 //prevent memory leak
768 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
769 lookupFloat = input->getSharedRAbundFloatVectors();
773 if (m->control_pressed) { return 0; }
775 //output error messages about any remaining user labels
776 set<string>::iterator it;
777 bool needToRun = false;
778 for (it = userLabels.begin(); it != userLabels.end(); it++) {
779 m->mothurOut("Your file does not include the label " + *it);
780 if (processedLabels.count(lastLabel) != 1) {
781 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
784 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
788 //run last label if you need to
789 if (needToRun == true) {
790 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
791 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
796 catch(exception& e) {
797 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
801 //**********************************************************************************************************************
802 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
805 vector<SharedRAbundFloatVector*> newLookup;
806 for (int i = 0; i < thislookup.size(); i++) {
807 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
808 temp->setLabel(thislookup[i]->getLabel());
809 temp->setGroup(thislookup[i]->getGroup());
810 newLookup.push_back(temp);
814 vector<string> newBinLabels;
815 string snumBins = toString(thislookup[0]->getNumBins());
816 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
817 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
819 //look at each sharedRabund and make sure they are not all zero
821 for (int j = 0; j < thislookup.size(); j++) {
822 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
825 //if they are not all zero add this bin
827 for (int j = 0; j < thislookup.size(); j++) {
828 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
831 //if there is a bin label use it otherwise make one
832 string binLabel = "Otu";
833 string sbinNumber = toString(i+1);
834 if (sbinNumber.length() < snumBins.length()) {
835 int diff = snumBins.length() - sbinNumber.length();
836 for (int h = 0; h < diff; h++) { binLabel += "0"; }
838 binLabel += sbinNumber;
839 if (i < m->currentBinLabels.size()) { binLabel = m->currentBinLabels[i]; }
841 newBinLabels.push_back(binLabel);
845 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
847 thislookup = newLookup;
848 m->currentBinLabels = newBinLabels;
853 catch(exception& e) {
854 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
858 /*****************************************************************/
859 map<string, vector<float> > CorrAxesCommand::readAxes(){
861 map<string, vector<float> > axes;
864 m->openInputFile(axesfile, in);
866 string headerLine = m->getline(in); m->gobble(in);
868 //count the number of axis you are reading
872 int pos = headerLine.find("axis");
873 if (pos != string::npos) {
875 headerLine = headerLine.substr(pos+4);
876 }else { done = true; }
879 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
883 if (m->control_pressed) { in.close(); return axes; }
886 in >> group; m->gobble(in);
888 vector<float> thisGroupsAxes;
889 for (int i = 0; i < count; i++) {
893 //only save the axis we want
894 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
897 //save group if its one the user selected
898 if (names.count(group) != 0) {
899 map<string, vector<float> >::iterator it = axes.find(group);
901 if (it == axes.end()) {
902 axes[group] = thisGroupsAxes;
904 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
914 catch(exception& e) {
915 m->errorOut(e, "CorrAxesCommand", "readAxes");
919 /*****************************************************************/
920 int CorrAxesCommand::getMetadata(){
922 vector<string> groupNames;
925 m->openInputFile(metadatafile, in);
927 string headerLine = m->getline(in); m->gobble(in);
928 vector<string> pieces = m->splitWhiteSpace(headerLine);
930 //save names of columns you are reading
931 for (int i = 1; i < pieces.size(); i++) {
932 metadataLabels.push_back(pieces[i]);
934 int count = metadataLabels.size();
939 if (m->control_pressed) { in.close(); return 0; }
942 in >> group; m->gobble(in);
943 groupNames.push_back(group);
945 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
946 tempLookup->setGroup(group);
947 tempLookup->setLabel("1");
949 for (int i = 0; i < count; i++) {
952 tempLookup->push_back(temp, group);
955 lookupFloat.push_back(tempLookup);
961 //remove any groups the user does not want, and set globaldata->groups with only valid groups
963 util = new SharedUtil();
964 Groups = m->getGroups();
965 util->setGroups(Groups, groupNames);
966 m->setGroups(Groups);
968 for (int i = 0; i < lookupFloat.size(); i++) {
969 //if this sharedrabund is not from a group the user wants then delete it.
970 if (util->isValidGroup(lookupFloat[i]->getGroup(), m->getGroups()) == false) {
971 delete lookupFloat[i]; lookupFloat[i] = NULL;
972 lookupFloat.erase(lookupFloat.begin()+i);
981 catch(exception& e) {
982 m->errorOut(e, "CorrAxesCommand", "getMetadata");
986 /*****************************************************************/