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","corraxes",false,true,true); parameters.push_back(paxes);
18 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none","",false,false,true); parameters.push_back(pshared);
19 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none","",false,false,true); 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::getOutputPattern(string type) {
63 if (type == "corraxes") { pattern = "[filename],[tag],corr.axes"; }
64 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
69 m->errorOut(e, "CorrAxesCommand", "getOutputPattern");
74 //**********************************************************************************************************************
75 CorrAxesCommand::CorrAxesCommand(){
77 abort = true; calledHelp = true;
79 vector<string> tempOutNames;
80 outputTypes["corraxes"] = tempOutNames;
83 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
87 //**********************************************************************************************************************
88 CorrAxesCommand::CorrAxesCommand(string option) {
90 abort = false; calledHelp = false;
92 //allow user to run help
93 if(option == "help") { help(); abort = true; calledHelp = true; }
94 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
97 vector<string> myArray = setParameters();
99 OptionParser parser(option);
100 map<string, string> parameters = parser.getParameters();
102 ValidParameters validParameter;
103 map<string, string>::iterator it;
105 //check to make sure all parameters are valid for command
106 for (it = parameters.begin(); it != parameters.end(); it++) {
107 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
110 vector<string> tempOutNames;
111 outputTypes["corraxes"] = tempOutNames;
113 //if the user changes the input directory command factory will send this info to us in the output parameter
114 string inputDir = validParameter.validFile(parameters, "inputdir", false);
115 if (inputDir == "not found"){ inputDir = ""; }
118 it = parameters.find("axes");
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["axes"] = inputDir + it->second; }
126 it = parameters.find("shared");
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["shared"] = inputDir + it->second; }
134 it = parameters.find("relabund");
135 //user has given a template file
136 if(it != parameters.end()){
137 path = m->hasPath(it->second);
138 //if the user has not given a path then, add inputdir. else leave path alone.
139 if (path == "") { parameters["relabund"] = inputDir + it->second; }
142 it = parameters.find("metadata");
143 //user has given a template file
144 if(it != parameters.end()){
145 path = m->hasPath(it->second);
146 //if the user has not given a path then, add inputdir. else leave path alone.
147 if (path == "") { parameters["metadata"] = inputDir + it->second; }
152 //check for required parameters
153 axesfile = validParameter.validFile(parameters, "axes", true);
154 if (axesfile == "not open") { abort = true; }
155 else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true; }
157 sharedfile = validParameter.validFile(parameters, "shared", true);
158 if (sharedfile == "not open") { abort = true; }
159 else if (sharedfile == "not found") { sharedfile = ""; }
160 else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
162 relabundfile = validParameter.validFile(parameters, "relabund", true);
163 if (relabundfile == "not open") { abort = true; }
164 else if (relabundfile == "not found") { relabundfile = ""; }
165 else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
167 metadatafile = validParameter.validFile(parameters, "metadata", true);
168 if (metadatafile == "not open") { abort = true; }
169 else if (metadatafile == "not found") { metadatafile = ""; }
170 else { inputFileName = metadatafile; }
172 groups = validParameter.validFile(parameters, "groups", false);
173 if (groups == "not found") { groups = ""; pickedGroups = false; }
176 m->splitAtDash(groups, Groups);
178 m->setGroups(Groups);
180 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
182 label = validParameter.validFile(parameters, "label", false);
183 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=""; }
185 if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) {
186 //is there are current file available for any of these?
187 //give priority to shared, then relabund
188 //if there is a current shared file, use it
189 sharedfile = m->getSharedFile();
190 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
192 relabundfile = m->getRelAbundFile();
193 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
195 m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true;
200 if (metadatafile != "") {
201 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
203 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
206 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
207 m->mothurConvert(temp, numaxes);
209 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
211 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; }
215 catch(exception& e) {
216 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
220 //**********************************************************************************************************************
222 int CorrAxesCommand::execute(){
225 if (abort == true) { if (calledHelp) { return 0; } return 2; }
227 /*************************************************************************************/
228 // use smart distancing to get right sharedRabund and convert to relabund if needed //
229 /************************************************************************************/
230 if (sharedfile != "") {
231 InputData* input = new InputData(sharedfile, "sharedfile");
232 getSharedFloat(input);
235 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
236 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
238 }else if (relabundfile != "") {
239 InputData* input = new InputData(relabundfile, "relabund");
240 getSharedFloat(input);
243 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
244 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
246 }else if (metadatafile != "") {
247 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
248 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
249 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
251 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
253 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
255 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
257 //this is for a sanity check to make sure the axes file and shared file match
258 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
260 /*************************************************************************************/
261 // read axes file and check for file mismatches with shared or relabund file //
262 /************************************************************************************/
265 map<string, vector<float> > axes = readAxes();
267 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
269 //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
270 if (axes.size() != lookupFloat.size()) {
271 map<string, vector<float> >::iterator it;
272 for (int i = 0; i < lookupFloat.size(); i++) {
273 it = axes.find(lookupFloat[i]->getGroup());
274 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(); }
276 m->control_pressed = true;
279 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
281 /*************************************************************************************/
282 // calc the r values //
283 /************************************************************************************/
284 map<string, string> variables;
285 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(inputFileName));
286 variables["[tag]"] = method;
287 string outputFileName = getOutputFileName("corraxes", variables);
288 outputNames.push_back(outputFileName); outputTypes["corraxes"].push_back(outputFileName);
290 m->openOutputFile(outputFileName, out);
291 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
294 if (metadatafile == "") { out << "OTU"; }
295 else { out << "Feature"; }
297 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1) << "\tp-value"; }
298 out << "\tlength" << endl;
300 if (method == "pearson") { calcPearson(axes, out); }
301 else if (method == "spearman") { calcSpearman(axes, out); }
302 else if (method == "kendall") { calcKendall(axes, out); }
303 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
306 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
308 if (m->control_pressed) { return 0; }
310 m->mothurOutEndLine();
311 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
312 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
313 m->mothurOutEndLine();
317 catch(exception& e) {
318 m->errorOut(e, "CorrAxesCommand", "execute");
322 //**********************************************************************************************************************
323 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
326 LinearAlgebra linear;
328 //find average of each axis - X
329 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
330 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
331 vector<float> temp = it->second;
332 for (int i = 0; i < temp.size(); i++) {
333 averageAxes[i] += temp[i];
337 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
340 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
342 if (metadatafile == "") { out << m->currentBinLabels[i]; }
343 else { out << metadataLabels[i]; }
345 //find the averages this otu - Y
347 for (int j = 0; j < lookupFloat.size(); j++) {
348 sumOtu += lookupFloat[j]->getAbundance(i);
350 float Ybar = sumOtu / (float) lookupFloat.size();
352 vector<float> rValues(averageAxes.size());
354 //find r value for each axis
355 for (int k = 0; k < averageAxes.size(); k++) {
358 double numerator = 0.0;
359 double denomTerm1 = 0.0;
360 double denomTerm2 = 0.0;
361 for (int j = 0; j < lookupFloat.size(); j++) {
362 float Yi = lookupFloat[j]->getAbundance(i);
363 float Xi = axes[lookupFloat[j]->getGroup()][k];
365 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
366 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
367 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
370 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
372 r = numerator / denom;
374 if (isnan(r) || isinf(r)) { r = 0.0; }
379 double sig = linear.calcPearsonSig(lookupFloat.size(), r);
385 for(int k=0;k<rValues.size();k++){
386 sum += rValues[k] * rValues[k];
388 out << '\t' << sqrt(sum) << endl;
393 catch(exception& e) {
394 m->errorOut(e, "CorrAxesCommand", "calcPearson");
398 //**********************************************************************************************************************
399 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
402 LinearAlgebra linear;
406 vector< map<float, int> > tableX; tableX.resize(numaxes);
407 map<float, int>::iterator itTable;
408 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
409 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
410 vector<float> temp = it->second;
411 for (int i = 0; i < temp.size(); i++) {
412 spearmanRank member(it->first, temp[i]);
413 scores[i].push_back(member);
415 //count number of repeats
416 itTable = tableX[i].find(temp[i]);
417 if (itTable == tableX[i].end()) {
418 tableX[i][temp[i]] = 1;
420 tableX[i][temp[i]]++;
427 vector<double> Lx; Lx.resize(numaxes, 0.0);
428 for (int i = 0; i < numaxes; i++) {
429 for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) {
430 double tx = (double) itTable->second;
431 Lx[i] += ((pow(tx, 3.0) - tx) / 12.0);
436 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
438 //find ranks of xi in each axis
439 map<string, vector<float> > rankAxes;
440 for (int i = 0; i < numaxes; i++) {
442 vector<spearmanRank> ties;
445 for (int j = 0; j < scores[i].size(); j++) {
447 ties.push_back(scores[i][j]);
449 if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
450 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
452 for (int k = 0; k < ties.size(); k++) {
453 float thisrank = rankTotal / (float) ties.size();
454 rankAxes[ties[k].name].push_back(thisrank);
461 }else { // you are the last one
463 for (int k = 0; k < ties.size(); k++) {
464 float thisrank = rankTotal / (float) ties.size();
465 rankAxes[ties[k].name].push_back(thisrank);
470 sf.push_back(sfTemp);
475 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
477 if (metadatafile == "") { out << m->currentBinLabels[i]; }
478 else { out << metadataLabels[i]; }
480 //find the ranks of this otu - Y
481 vector<spearmanRank> otuScores;
482 map<float, int> tableY;
483 for (int j = 0; j < lookupFloat.size(); j++) {
484 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
485 otuScores.push_back(member);
487 itTable = tableY.find(member.score);
488 if (itTable == tableY.end()) {
489 tableY[member.score] = 1;
491 tableY[member.score]++;
498 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
499 double ty = (double) itTable->second;
500 Ly += ((pow(ty, 3.0) - ty) / 12.0);
503 sort(otuScores.begin(), otuScores.end(), compareSpearman);
506 map<string, float> rankOtus;
507 vector<spearmanRank> ties;
509 for (int j = 0; j < otuScores.size(); j++) {
511 ties.push_back(otuScores[j]);
513 if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
514 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
516 for (int k = 0; k < ties.size(); k++) {
517 float thisrank = rankTotal / (float) ties.size();
518 rankOtus[ties[k].name] = thisrank;
525 }else { // you are the last one
527 for (int k = 0; k < ties.size(); k++) {
528 float thisrank = rankTotal / (float) ties.size();
529 rankOtus[ties[k].name] = thisrank;
533 vector<double> pValues(numaxes);
535 //calc spearman ranks for each axis for this otu
536 for (int j = 0; j < numaxes; j++) {
539 for (int k = 0; k < lookupFloat.size(); k++) {
541 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
542 float yi = rankOtus[lookupFloat[k]->getGroup()];
544 di += ((xi - yi) * (xi - yi));
549 double n = (double) lookupFloat.size();
551 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
552 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
554 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
556 if (isnan(p) || isinf(p)) { p = 0.0; }
562 double sig = linear.calcSpearmanSig(n, sf[j], sg, di);
568 for(int k=0;k<numaxes;k++){
569 sum += pValues[k] * pValues[k];
571 out << '\t' << sqrt(sum) << endl;
576 catch(exception& e) {
577 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
581 //**********************************************************************************************************************
582 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
585 LinearAlgebra linear;
588 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
589 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
590 vector<float> temp = it->second;
591 for (int i = 0; i < temp.size(); i++) {
592 spearmanRank member(it->first, temp[i]);
593 scores[i].push_back(member);
598 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
600 //convert scores to ranks of xi in each axis
601 for (int i = 0; i < numaxes; i++) {
603 vector<spearmanRank*> ties;
605 for (int j = 0; j < scores[i].size(); j++) {
607 ties.push_back(&(scores[i][j]));
609 if (j != scores[i].size()-1) { // you are not the last so you can look ahead
610 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
611 for (int k = 0; k < ties.size(); k++) {
612 float thisrank = rankTotal / (float) ties.size();
613 (*ties[k]).score = thisrank;
618 }else { // you are the last one
619 for (int k = 0; k < ties.size(); k++) {
620 float thisrank = rankTotal / (float) ties.size();
621 (*ties[k]).score = thisrank;
628 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
630 if (metadatafile == "") { out << m->currentBinLabels[i]; }
631 else { out << metadataLabels[i]; }
633 //find the ranks of this otu - Y
634 vector<spearmanRank> otuScores;
635 for (int j = 0; j < lookupFloat.size(); j++) {
636 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
637 otuScores.push_back(member);
640 sort(otuScores.begin(), otuScores.end(), compareSpearman);
642 map<string, float> rankOtus;
643 vector<spearmanRank> ties;
645 for (int j = 0; j < otuScores.size(); j++) {
647 ties.push_back(otuScores[j]);
649 if (j != otuScores.size()-1) { // you are not the last so you can look ahead
650 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
651 for (int k = 0; k < ties.size(); k++) {
652 float thisrank = rankTotal / (float) ties.size();
653 rankOtus[ties[k].name] = thisrank;
658 }else { // you are the last one
659 for (int k = 0; k < ties.size(); k++) {
660 float thisrank = rankTotal / (float) ties.size();
661 rankOtus[ties[k].name] = thisrank;
667 vector<double> pValues(numaxes);
669 //calc spearman ranks for each axis for this otu
670 for (int j = 0; j < numaxes; j++) {
675 vector<spearmanRank> otus;
676 vector<spearmanRank> otusTemp;
677 for (int l = 0; l < scores[j].size(); l++) {
678 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
679 otus.push_back(member);
683 for (int l = 0; l < scores[j].size(); l++) {
685 int numWithHigherRank = 0;
686 int numWithLowerRank = 0;
687 float thisrank = otus[l].score;
689 for (int u = l+1; u < scores[j].size(); u++) {
690 if (otus[u].score > thisrank) { numWithHigherRank++; }
691 else if (otus[u].score < thisrank) { numWithLowerRank++; }
695 numCoor += numWithHigherRank;
696 numDisCoor += numWithLowerRank;
699 double p = (numCoor - numDisCoor) / (float) count;
700 if (isnan(p) || isinf(p)) { p = 0.0; }
705 double sig = linear.calcKendallSig(scores[j].size(), p);
711 for(int k=0;k<numaxes;k++){
712 sum += pValues[k] * pValues[k];
714 out << '\t' << sqrt(sum) << endl;
719 catch(exception& e) {
720 m->errorOut(e, "CorrAxesCommand", "calcKendall");
724 //**********************************************************************************************************************
725 int CorrAxesCommand::getSharedFloat(InputData* input){
727 lookupFloat = input->getSharedRAbundFloatVectors();
728 string lastLabel = lookupFloat[0]->getLabel();
730 if (label == "") { label = lastLabel; return 0; }
732 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
733 set<string> labels; labels.insert(label);
734 set<string> processedLabels;
735 set<string> userLabels = labels;
737 //as long as you are not at the end of the file or done wih the lines you want
738 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
740 if (m->control_pressed) { return 0; }
742 if(labels.count(lookupFloat[0]->getLabel()) == 1){
743 processedLabels.insert(lookupFloat[0]->getLabel());
744 userLabels.erase(lookupFloat[0]->getLabel());
748 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
749 string saveLabel = lookupFloat[0]->getLabel();
751 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
752 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
754 processedLabels.insert(lookupFloat[0]->getLabel());
755 userLabels.erase(lookupFloat[0]->getLabel());
757 //restore real lastlabel to save below
758 lookupFloat[0]->setLabel(saveLabel);
762 lastLabel = lookupFloat[0]->getLabel();
764 //get next line to process
765 //prevent memory leak
766 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
767 lookupFloat = input->getSharedRAbundFloatVectors();
771 if (m->control_pressed) { return 0; }
773 //output error messages about any remaining user labels
774 set<string>::iterator it;
775 bool needToRun = false;
776 for (it = userLabels.begin(); it != userLabels.end(); it++) {
777 m->mothurOut("Your file does not include the label " + *it);
778 if (processedLabels.count(lastLabel) != 1) {
779 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
782 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
786 //run last label if you need to
787 if (needToRun == true) {
788 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
789 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
794 catch(exception& e) {
795 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
799 //**********************************************************************************************************************
800 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
803 vector<SharedRAbundFloatVector*> newLookup;
804 for (int i = 0; i < thislookup.size(); i++) {
805 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
806 temp->setLabel(thislookup[i]->getLabel());
807 temp->setGroup(thislookup[i]->getGroup());
808 newLookup.push_back(temp);
812 vector<string> newBinLabels;
813 string snumBins = toString(thislookup[0]->getNumBins());
814 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
815 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
817 //look at each sharedRabund and make sure they are not all zero
819 for (int j = 0; j < thislookup.size(); j++) {
820 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
823 //if they are not all zero add this bin
825 for (int j = 0; j < thislookup.size(); j++) {
826 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
829 //if there is a bin label use it otherwise make one
830 string binLabel = "Otu";
831 string sbinNumber = toString(i+1);
832 if (sbinNumber.length() < snumBins.length()) {
833 int diff = snumBins.length() - sbinNumber.length();
834 for (int h = 0; h < diff; h++) { binLabel += "0"; }
836 binLabel += sbinNumber;
837 if (i < m->currentBinLabels.size()) { binLabel = m->currentBinLabels[i]; }
839 newBinLabels.push_back(binLabel);
843 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
845 thislookup = newLookup;
846 m->currentBinLabels = newBinLabels;
851 catch(exception& e) {
852 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
856 /*****************************************************************/
857 map<string, vector<float> > CorrAxesCommand::readAxes(){
859 map<string, vector<float> > axes;
862 m->openInputFile(axesfile, in);
864 string headerLine = m->getline(in); m->gobble(in);
866 //count the number of axis you are reading
870 int pos = headerLine.find("axis");
871 if (pos != string::npos) {
873 headerLine = headerLine.substr(pos+4);
874 }else { done = true; }
877 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
881 if (m->control_pressed) { in.close(); return axes; }
884 in >> group; m->gobble(in);
886 vector<float> thisGroupsAxes;
887 for (int i = 0; i < count; i++) {
891 //only save the axis we want
892 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
895 //save group if its one the user selected
896 if (names.count(group) != 0) {
897 map<string, vector<float> >::iterator it = axes.find(group);
899 if (it == axes.end()) {
900 axes[group] = thisGroupsAxes;
902 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
912 catch(exception& e) {
913 m->errorOut(e, "CorrAxesCommand", "readAxes");
917 /*****************************************************************/
918 int CorrAxesCommand::getMetadata(){
920 vector<string> groupNames;
923 m->openInputFile(metadatafile, in);
925 string headerLine = m->getline(in); m->gobble(in);
926 vector<string> pieces = m->splitWhiteSpace(headerLine);
928 //save names of columns you are reading
929 for (int i = 1; i < pieces.size(); i++) {
930 metadataLabels.push_back(pieces[i]);
932 int count = metadataLabels.size();
937 if (m->control_pressed) { in.close(); return 0; }
940 in >> group; m->gobble(in);
941 groupNames.push_back(group);
943 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
944 tempLookup->setGroup(group);
945 tempLookup->setLabel("1");
947 for (int i = 0; i < count; i++) {
950 tempLookup->push_back(temp, group);
953 lookupFloat.push_back(tempLookup);
959 //remove any groups the user does not want, and set globaldata->groups with only valid groups
961 util = new SharedUtil();
962 Groups = m->getGroups();
963 util->setGroups(Groups, groupNames);
964 m->setGroups(Groups);
966 for (int i = 0; i < lookupFloat.size(); i++) {
967 //if this sharedrabund is not from a group the user wants then delete it.
968 if (util->isValidGroup(lookupFloat[i]->getGroup(), m->getGroups()) == false) {
969 delete lookupFloat[i]; lookupFloat[i] = NULL;
970 lookupFloat.erase(lookupFloat.begin()+i);
979 catch(exception& e) {
980 m->errorOut(e, "CorrAxesCommand", "getMetadata");
984 /*****************************************************************/