5 * Created by westcott on 12/22/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "corraxescommand.h"
11 #include "sharedutilities.h"
13 //**********************************************************************************************************************
14 vector<string> CorrAxesCommand::setParameters(){
16 CommandParameter paxes("axes", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(paxes);
17 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pshared);
18 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(prelabund);
19 CommandParameter pmetadata("metadata", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pmetadata);
20 CommandParameter pnumaxes("numaxes", "Number", "", "3", "", "", "",false,false); parameters.push_back(pnumaxes);
21 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
22 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
23 CommandParameter pmethod("method", "Multiple", "pearson-spearman-kendall", "pearson", "", "", "",false,false); parameters.push_back(pmethod);
24 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
25 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
27 vector<string> myArray;
28 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
32 m->errorOut(e, "CorrAxesCommand", "setParameters");
36 //**********************************************************************************************************************
37 string CorrAxesCommand::getHelpString(){
39 string helpString = "";
40 helpString += "The corr.axes command reads a shared, relabund or metadata file as well as an axes file and calculates the correlation coefficient.\n";
41 helpString += "The corr.axes command parameters are shared, relabund, axes, metadata, groups, method, numaxes and label. The shared, relabund or metadata and axes parameters are required. If shared is given the relative abundance is calculated.\n";
42 helpString += "The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes.\n";
43 helpString += "The label parameter allows you to select what distance level you would like used, if none is given the first distance is used.\n";
44 helpString += "The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n";
45 helpString += "The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n";
46 helpString += "The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n";
47 helpString += "Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n";
48 helpString += "The corr.axes command outputs a .corr.axes file.\n";
49 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
53 m->errorOut(e, "CorrAxesCommand", "getHelpString");
57 //**********************************************************************************************************************
58 CorrAxesCommand::CorrAxesCommand(){
60 abort = true; calledHelp = true;
62 vector<string> tempOutNames;
63 outputTypes["corr.axes"] = tempOutNames;
66 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
70 //**********************************************************************************************************************
71 CorrAxesCommand::CorrAxesCommand(string option) {
73 abort = false; calledHelp = false;
75 //allow user to run help
76 if(option == "help") { help(); abort = true; calledHelp = true; }
77 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
80 vector<string> myArray = setParameters();
82 OptionParser parser(option);
83 map<string, string> parameters = parser.getParameters();
85 ValidParameters validParameter;
86 map<string, string>::iterator it;
88 //check to make sure all parameters are valid for command
89 for (it = parameters.begin(); it != parameters.end(); it++) {
90 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
93 vector<string> tempOutNames;
94 outputTypes["corr.axes"] = tempOutNames;
96 //if the user changes the input directory command factory will send this info to us in the output parameter
97 string inputDir = validParameter.validFile(parameters, "inputdir", false);
98 if (inputDir == "not found"){ inputDir = ""; }
101 it = parameters.find("axes");
102 //user has given a template file
103 if(it != parameters.end()){
104 path = m->hasPath(it->second);
105 //if the user has not given a path then, add inputdir. else leave path alone.
106 if (path == "") { parameters["axes"] = inputDir + it->second; }
109 it = parameters.find("shared");
110 //user has given a template file
111 if(it != parameters.end()){
112 path = m->hasPath(it->second);
113 //if the user has not given a path then, add inputdir. else leave path alone.
114 if (path == "") { parameters["shared"] = inputDir + it->second; }
117 it = parameters.find("relabund");
118 //user has given a template file
119 if(it != parameters.end()){
120 path = m->hasPath(it->second);
121 //if the user has not given a path then, add inputdir. else leave path alone.
122 if (path == "") { parameters["relabund"] = inputDir + it->second; }
125 it = parameters.find("metadata");
126 //user has given a template file
127 if(it != parameters.end()){
128 path = m->hasPath(it->second);
129 //if the user has not given a path then, add inputdir. else leave path alone.
130 if (path == "") { parameters["metadata"] = inputDir + it->second; }
135 //check for required parameters
136 axesfile = validParameter.validFile(parameters, "axes", true);
137 if (axesfile == "not open") { abort = true; }
138 else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true; }
140 sharedfile = validParameter.validFile(parameters, "shared", true);
141 if (sharedfile == "not open") { abort = true; }
142 else if (sharedfile == "not found") { sharedfile = ""; }
143 else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
145 relabundfile = validParameter.validFile(parameters, "relabund", true);
146 if (relabundfile == "not open") { abort = true; }
147 else if (relabundfile == "not found") { relabundfile = ""; }
148 else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
150 metadatafile = validParameter.validFile(parameters, "metadata", true);
151 if (metadatafile == "not open") { abort = true; }
152 else if (metadatafile == "not found") { metadatafile = ""; }
153 else { inputFileName = metadatafile; }
155 groups = validParameter.validFile(parameters, "groups", false);
156 if (groups == "not found") { groups = ""; pickedGroups = false; }
159 m->splitAtDash(groups, Groups);
161 m->setGroups(Groups);
163 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
165 label = validParameter.validFile(parameters, "label", false);
166 if (label == "not found") { label = ""; m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); label=""; }
168 if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) {
169 //is there are current file available for any of these?
170 //give priority to shared, then relabund
171 //if there is a current shared file, use it
172 sharedfile = m->getSharedFile();
173 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
175 relabundfile = m->getRelAbundFile();
176 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
178 m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true;
183 if (metadatafile != "") {
184 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
186 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
189 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
190 m->mothurConvert(temp, numaxes);
192 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
194 if ((method != "pearson") && (method != "spearman") && (method != "kendall")) { m->mothurOut(method + " is not a valid method. Valid methods are pearson, spearman, and kendall."); m->mothurOutEndLine(); abort = true; }
198 catch(exception& e) {
199 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
203 //**********************************************************************************************************************
205 int CorrAxesCommand::execute(){
208 if (abort == true) { if (calledHelp) { return 0; } return 2; }
210 /*************************************************************************************/
211 // use smart distancing to get right sharedRabund and convert to relabund if needed //
212 /************************************************************************************/
213 if (sharedfile != "") {
214 InputData* input = new InputData(sharedfile, "sharedfile");
215 getSharedFloat(input);
218 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
219 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
221 }else if (relabundfile != "") {
222 InputData* input = new InputData(relabundfile, "relabund");
223 getSharedFloat(input);
226 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
227 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
229 }else if (metadatafile != "") {
230 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
231 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
232 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
234 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
236 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
238 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
240 //this is for a sanity check to make sure the axes file and shared file match
241 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
243 /*************************************************************************************/
244 // read axes file and check for file mismatches with shared or relabund file //
245 /************************************************************************************/
248 map<string, vector<float> > axes = readAxes();
250 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
252 //sanity check, the read only adds groups that are in the shared or relabund file, but we want to make sure the axes file isn't missing anyone
253 if (axes.size() != lookupFloat.size()) {
254 map<string, vector<float> >::iterator it;
255 for (int i = 0; i < lookupFloat.size(); i++) {
256 it = axes.find(lookupFloat[i]->getGroup());
257 if (it == axes.end()) { m->mothurOut(lookupFloat[i]->getGroup() + " is in your shared of relabund file but not in your axes file, please correct."); m->mothurOutEndLine(); }
259 m->control_pressed = true;
262 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
264 /*************************************************************************************/
265 // calc the r values //
266 /************************************************************************************/
268 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
269 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
271 m->openOutputFile(outputFileName, out);
272 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
275 if (metadatafile == "") { out << "OTU"; }
276 else { out << "Feature"; }
278 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1) << "\tp-value"; }
279 out << "\tlength" << endl;
281 if (method == "pearson") { calcPearson(axes, out); }
282 else if (method == "spearman") { calcSpearman(axes, out); }
283 else if (method == "kendall") { calcKendall(axes, out); }
284 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
287 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
289 if (m->control_pressed) { return 0; }
291 m->mothurOutEndLine();
292 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
293 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
294 m->mothurOutEndLine();
298 catch(exception& e) {
299 m->errorOut(e, "CorrAxesCommand", "execute");
303 //**********************************************************************************************************************
304 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
307 //find average of each axis - X
308 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
309 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
310 vector<float> temp = it->second;
311 for (int i = 0; i < temp.size(); i++) {
312 averageAxes[i] += temp[i];
316 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
319 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
321 if (metadatafile == "") { out << i+1; }
322 else { out << metadataLabels[i]; }
324 //find the averages this otu - Y
326 for (int j = 0; j < lookupFloat.size(); j++) {
327 sumOtu += lookupFloat[j]->getAbundance(i);
329 float Ybar = sumOtu / (float) lookupFloat.size();
331 vector<float> rValues(averageAxes.size());
333 //find r value for each axis
334 for (int k = 0; k < averageAxes.size(); k++) {
337 double numerator = 0.0;
338 double denomTerm1 = 0.0;
339 double denomTerm2 = 0.0;
340 for (int j = 0; j < lookupFloat.size(); j++) {
341 float Yi = lookupFloat[j]->getAbundance(i);
342 float Xi = axes[lookupFloat[j]->getGroup()][k];
344 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
345 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
346 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
349 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
351 r = numerator / denom;
353 if (isnan(r) || isinf(r)) { r = 0.0; }
358 //signifigance calc - http://faculty.vassar.edu/lowry/ch4apx.html
359 double temp = (1- (r*r)) / (double) (lookupFloat.size()-2);
361 double sig = r / temp;
362 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
368 for(int k=0;k<rValues.size();k++){
369 sum += rValues[k] * rValues[k];
371 out << '\t' << sqrt(sum) << endl;
376 catch(exception& e) {
377 m->errorOut(e, "CorrAxesCommand", "calcPearson");
381 //**********************************************************************************************************************
382 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
386 vector< map<float, int> > tableX; tableX.resize(numaxes);
387 map<float, int>::iterator itTable;
388 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
389 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
390 vector<float> temp = it->second;
391 for (int i = 0; i < temp.size(); i++) {
392 spearmanRank member(it->first, temp[i]);
393 scores[i].push_back(member);
395 //count number of repeats
396 itTable = tableX[i].find(temp[i]);
397 if (itTable == tableX[i].end()) {
398 tableX[i][temp[i]] = 1;
400 tableX[i][temp[i]]++;
407 vector<double> Lx; Lx.resize(numaxes, 0.0);
408 for (int i = 0; i < numaxes; i++) {
409 for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) {
410 double tx = (double) itTable->second;
411 Lx[i] += ((pow(tx, 3.0) - tx) / 12.0);
416 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
418 //find ranks of xi in each axis
419 map<string, vector<float> > rankAxes;
420 for (int i = 0; i < numaxes; i++) {
422 vector<spearmanRank> ties;
424 for (int j = 0; j < scores[i].size(); j++) {
426 ties.push_back(scores[i][j]);
428 if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
429 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
431 for (int k = 0; k < ties.size(); k++) {
432 float thisrank = rankTotal / (float) ties.size();
433 rankAxes[ties[k].name].push_back(thisrank);
438 }else { // you are the last one
440 for (int k = 0; k < ties.size(); k++) {
441 float thisrank = rankTotal / (float) ties.size();
442 rankAxes[ties[k].name].push_back(thisrank);
451 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
453 if (metadatafile == "") { out << i+1; }
454 else { out << metadataLabels[i]; }
456 //find the ranks of this otu - Y
457 vector<spearmanRank> otuScores;
458 map<float, int> tableY;
459 for (int j = 0; j < lookupFloat.size(); j++) {
460 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
461 otuScores.push_back(member);
463 itTable = tableY.find(member.score);
464 if (itTable == tableY.end()) {
465 tableY[member.score] = 1;
467 tableY[member.score]++;
474 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
475 double ty = (double) itTable->second;
476 Ly += ((pow(ty, 3.0) - ty) / 12.0);
479 sort(otuScores.begin(), otuScores.end(), compareSpearman);
481 map<string, float> rankOtus;
482 vector<spearmanRank> ties;
484 for (int j = 0; j < otuScores.size(); j++) {
486 ties.push_back(otuScores[j]);
488 if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
489 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
491 for (int k = 0; k < ties.size(); k++) {
492 float thisrank = rankTotal / (float) ties.size();
493 rankOtus[ties[k].name] = thisrank;
498 }else { // you are the last one
500 for (int k = 0; k < ties.size(); k++) {
501 float thisrank = rankTotal / (float) ties.size();
502 rankOtus[ties[k].name] = thisrank;
506 vector<double> pValues(numaxes);
508 //calc spearman ranks for each axis for this otu
509 for (int j = 0; j < numaxes; j++) {
512 for (int k = 0; k < lookupFloat.size(); k++) {
514 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
515 float yi = rankOtus[lookupFloat[k]->getGroup()];
517 di += ((xi - yi) * (xi - yi));
522 double n = (double) lookupFloat.size();
524 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
525 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
527 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
529 if (isnan(p) || isinf(p)) { p = 0.0; }
535 //signifigance calc - http://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient
536 double temp = (lookupFloat.size()-2) / (double) (1- (p*p));
539 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
546 for(int k=0;k<numaxes;k++){
547 sum += pValues[k] * pValues[k];
549 out << '\t' << sqrt(sum) << endl;
554 catch(exception& e) {
555 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
559 //**********************************************************************************************************************
560 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
564 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
565 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
566 vector<float> temp = it->second;
567 for (int i = 0; i < temp.size(); i++) {
568 spearmanRank member(it->first, temp[i]);
569 scores[i].push_back(member);
574 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
576 //convert scores to ranks of xi in each axis
577 for (int i = 0; i < numaxes; i++) {
579 vector<spearmanRank*> ties;
581 for (int j = 0; j < scores[i].size(); j++) {
583 ties.push_back(&(scores[i][j]));
585 if (j != scores[i].size()-1) { // you are not the last so you can look ahead
586 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
587 for (int k = 0; k < ties.size(); k++) {
588 float thisrank = rankTotal / (float) ties.size();
589 (*ties[k]).score = thisrank;
594 }else { // you are the last one
595 for (int k = 0; k < ties.size(); k++) {
596 float thisrank = rankTotal / (float) ties.size();
597 (*ties[k]).score = thisrank;
604 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
606 if (metadatafile == "") { out << i+1; }
607 else { out << metadataLabels[i]; }
609 //find the ranks of this otu - Y
610 vector<spearmanRank> otuScores;
611 for (int j = 0; j < lookupFloat.size(); j++) {
612 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
613 otuScores.push_back(member);
616 sort(otuScores.begin(), otuScores.end(), compareSpearman);
618 map<string, float> rankOtus;
619 vector<spearmanRank> ties;
621 for (int j = 0; j < otuScores.size(); j++) {
623 ties.push_back(otuScores[j]);
625 if (j != otuScores.size()-1) { // you are not the last so you can look ahead
626 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
627 for (int k = 0; k < ties.size(); k++) {
628 float thisrank = rankTotal / (float) ties.size();
629 rankOtus[ties[k].name] = thisrank;
634 }else { // you are the last one
635 for (int k = 0; k < ties.size(); k++) {
636 float thisrank = rankTotal / (float) ties.size();
637 rankOtus[ties[k].name] = thisrank;
643 vector<double> pValues(numaxes);
645 //calc spearman ranks for each axis for this otu
646 for (int j = 0; j < numaxes; j++) {
651 vector<spearmanRank> otus;
652 vector<spearmanRank> otusTemp;
653 for (int l = 0; l < scores[j].size(); l++) {
654 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
655 otus.push_back(member);
659 for (int l = 0; l < scores[j].size(); l++) {
661 int numWithHigherRank = 0;
662 int numWithLowerRank = 0;
663 float thisrank = otus[l].score;
665 for (int u = l+1; u < scores[j].size(); u++) {
666 if (otus[u].score > thisrank) { numWithHigherRank++; }
667 else if (otus[u].score < thisrank) { numWithLowerRank++; }
671 numCoor += numWithHigherRank;
672 numDisCoor += numWithLowerRank;
675 double p = (numCoor - numDisCoor) / (float) count;
676 if (isnan(p) || isinf(p)) { p = 0.0; }
681 //calc signif - zA - http://en.wikipedia.org/wiki/Kendall_tau_rank_correlation_coefficient#Significance_tests
682 double numer = 3.0 * (numCoor - numDisCoor);
683 int n = scores[j].size();
684 double denom = n * (n-1) * (2*n + 5) / (double) 2.0;
686 double sig = numer / denom;
688 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
694 for(int k=0;k<numaxes;k++){
695 sum += pValues[k] * pValues[k];
697 out << '\t' << sqrt(sum) << endl;
702 catch(exception& e) {
703 m->errorOut(e, "CorrAxesCommand", "calcKendall");
707 //**********************************************************************************************************************
708 int CorrAxesCommand::getSharedFloat(InputData* input){
710 lookupFloat = input->getSharedRAbundFloatVectors();
711 string lastLabel = lookupFloat[0]->getLabel();
713 if (label == "") { label = lastLabel; return 0; }
715 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
716 set<string> labels; labels.insert(label);
717 set<string> processedLabels;
718 set<string> userLabels = labels;
720 //as long as you are not at the end of the file or done wih the lines you want
721 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
723 if (m->control_pressed) { return 0; }
725 if(labels.count(lookupFloat[0]->getLabel()) == 1){
726 processedLabels.insert(lookupFloat[0]->getLabel());
727 userLabels.erase(lookupFloat[0]->getLabel());
731 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
732 string saveLabel = lookupFloat[0]->getLabel();
734 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
735 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
737 processedLabels.insert(lookupFloat[0]->getLabel());
738 userLabels.erase(lookupFloat[0]->getLabel());
740 //restore real lastlabel to save below
741 lookupFloat[0]->setLabel(saveLabel);
745 lastLabel = lookupFloat[0]->getLabel();
747 //get next line to process
748 //prevent memory leak
749 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
750 lookupFloat = input->getSharedRAbundFloatVectors();
754 if (m->control_pressed) { return 0; }
756 //output error messages about any remaining user labels
757 set<string>::iterator it;
758 bool needToRun = false;
759 for (it = userLabels.begin(); it != userLabels.end(); it++) {
760 m->mothurOut("Your file does not include the label " + *it);
761 if (processedLabels.count(lastLabel) != 1) {
762 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
765 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
769 //run last label if you need to
770 if (needToRun == true) {
771 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
772 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
777 catch(exception& e) {
778 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
782 //**********************************************************************************************************************
783 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
786 vector<SharedRAbundFloatVector*> newLookup;
787 for (int i = 0; i < thislookup.size(); i++) {
788 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
789 temp->setLabel(thislookup[i]->getLabel());
790 temp->setGroup(thislookup[i]->getGroup());
791 newLookup.push_back(temp);
795 vector<string> newBinLabels;
796 string snumBins = toString(thislookup[0]->getNumBins());
797 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
798 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
800 //look at each sharedRabund and make sure they are not all zero
802 for (int j = 0; j < thislookup.size(); j++) {
803 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
806 //if they are not all zero add this bin
808 for (int j = 0; j < thislookup.size(); j++) {
809 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
812 //if there is a bin label use it otherwise make one
813 string binLabel = "Otu";
814 string sbinNumber = toString(i+1);
815 if (sbinNumber.length() < snumBins.length()) {
816 int diff = snumBins.length() - sbinNumber.length();
817 for (int h = 0; h < diff; h++) { binLabel += "0"; }
819 binLabel += sbinNumber;
820 if (i < m->currentBinLabels.size()) { binLabel = m->currentBinLabels[i]; }
822 newBinLabels.push_back(binLabel);
826 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
828 thislookup = newLookup;
829 m->currentBinLabels = newBinLabels;
834 catch(exception& e) {
835 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
839 /*****************************************************************/
840 map<string, vector<float> > CorrAxesCommand::readAxes(){
842 map<string, vector<float> > axes;
845 m->openInputFile(axesfile, in);
847 string headerLine = m->getline(in); m->gobble(in);
849 //count the number of axis you are reading
853 int pos = headerLine.find("axis");
854 if (pos != string::npos) {
856 headerLine = headerLine.substr(pos+4);
857 }else { done = true; }
860 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
864 if (m->control_pressed) { in.close(); return axes; }
867 in >> group; m->gobble(in);
869 vector<float> thisGroupsAxes;
870 for (int i = 0; i < count; i++) {
874 //only save the axis we want
875 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
878 //save group if its one the user selected
879 if (names.count(group) != 0) {
880 map<string, vector<float> >::iterator it = axes.find(group);
882 if (it == axes.end()) {
883 axes[group] = thisGroupsAxes;
885 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
895 catch(exception& e) {
896 m->errorOut(e, "CorrAxesCommand", "readAxes");
900 /*****************************************************************/
901 int CorrAxesCommand::getMetadata(){
903 vector<string> groupNames;
906 m->openInputFile(metadatafile, in);
908 string headerLine = m->getline(in); m->gobble(in);
909 istringstream iss (headerLine,istringstream::in);
911 //read the first label, because it refers to the groups
913 iss >> columnLabel; m->gobble(iss);
915 //save names of columns you are reading
917 iss >> columnLabel; m->gobble(iss);
918 metadataLabels.push_back(columnLabel);
920 int count = metadataLabels.size();
925 if (m->control_pressed) { in.close(); return 0; }
928 in >> group; m->gobble(in);
929 groupNames.push_back(group);
931 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
932 tempLookup->setGroup(group);
933 tempLookup->setLabel("1");
935 for (int i = 0; i < count; i++) {
938 tempLookup->push_back(temp, group);
941 lookupFloat.push_back(tempLookup);
947 //remove any groups the user does not want, and set globaldata->groups with only valid groups
949 util = new SharedUtil();
950 Groups = m->getGroups();
951 util->setGroups(Groups, groupNames);
952 m->setGroups(Groups);
954 for (int i = 0; i < lookupFloat.size(); i++) {
955 //if this sharedrabund is not from a group the user wants then delete it.
956 if (util->isValidGroup(lookupFloat[i]->getGroup(), m->getGroups()) == false) {
957 delete lookupFloat[i]; lookupFloat[i] = NULL;
958 lookupFloat.erase(lookupFloat.begin()+i);
967 catch(exception& e) {
968 m->errorOut(e, "CorrAxesCommand", "getMetadata");
972 /*****************************************************************/