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 //sorts highest to lowest
15 inline bool compareSpearman(spearmanRank left, spearmanRank right){
16 return (left.score > right.score);
18 //********************************************************************************************************************
19 //sorts lowest to highest
20 inline bool compareSpearmanReverse(spearmanRank left, spearmanRank right){
21 return (left.score < right.score);
23 //**********************************************************************************************************************
24 vector<string> CorrAxesCommand::getValidParameters(){
26 string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
27 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
31 m->errorOut(e, "CorrAxesCommand", "getValidParameters");
35 //**********************************************************************************************************************
36 vector<string> CorrAxesCommand::getRequiredParameters(){
38 string Array[] = {"axes"};
39 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
43 m->errorOut(e, "CorrAxesCommand", "getRequiredParameters");
47 //**********************************************************************************************************************
48 CorrAxesCommand::CorrAxesCommand(){
51 //initialize outputTypes
52 vector<string> tempOutNames;
53 outputTypes["corr.axes"] = tempOutNames;
56 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
61 //**********************************************************************************************************************
62 vector<string> CorrAxesCommand::getRequiredFiles(){
64 vector<string> myArray;
68 m->errorOut(e, "CorrAxesCommand", "getRequiredFiles");
72 //**********************************************************************************************************************
73 CorrAxesCommand::CorrAxesCommand(string option) {
76 globaldata = GlobalData::getInstance();
78 //allow user to run help
79 if(option == "help") { help(); abort = true; }
82 //valid paramters for this command
83 string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
84 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
86 OptionParser parser(option);
87 map<string, string> parameters = parser.getParameters();
89 ValidParameters validParameter;
90 map<string, string>::iterator it;
92 //check to make sure all parameters are valid for command
93 for (it = parameters.begin(); it != parameters.end(); it++) {
94 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
97 vector<string> tempOutNames;
98 outputTypes["corr.axes"] = tempOutNames;
100 //if the user changes the input directory command factory will send this info to us in the output parameter
101 string inputDir = validParameter.validFile(parameters, "inputdir", false);
102 if (inputDir == "not found"){ inputDir = ""; }
105 it = parameters.find("axes");
106 //user has given a template file
107 if(it != parameters.end()){
108 path = m->hasPath(it->second);
109 //if the user has not given a path then, add inputdir. else leave path alone.
110 if (path == "") { parameters["axes"] = inputDir + it->second; }
113 it = parameters.find("shared");
114 //user has given a template file
115 if(it != parameters.end()){
116 path = m->hasPath(it->second);
117 //if the user has not given a path then, add inputdir. else leave path alone.
118 if (path == "") { parameters["shared"] = inputDir + it->second; }
121 it = parameters.find("relabund");
122 //user has given a template file
123 if(it != parameters.end()){
124 path = m->hasPath(it->second);
125 //if the user has not given a path then, add inputdir. else leave path alone.
126 if (path == "") { parameters["relabund"] = inputDir + it->second; }
129 it = parameters.find("metadata");
130 //user has given a template file
131 if(it != parameters.end()){
132 path = m->hasPath(it->second);
133 //if the user has not given a path then, add inputdir. else leave path alone.
134 if (path == "") { parameters["metadata"] = inputDir + it->second; }
139 //check for required parameters
140 axesfile = validParameter.validFile(parameters, "axes", true);
141 if (axesfile == "not open") { abort = true; }
142 else if (axesfile == "not found") { axesfile = ""; m->mothurOut("axes is a required parameter for the corr.axes command."); m->mothurOutEndLine(); abort = true; }
144 sharedfile = validParameter.validFile(parameters, "shared", true);
145 if (sharedfile == "not open") { abort = true; }
146 else if (sharedfile == "not found") { sharedfile = ""; }
147 else { inputFileName = sharedfile; }
149 relabundfile = validParameter.validFile(parameters, "relabund", true);
150 if (relabundfile == "not open") { abort = true; }
151 else if (relabundfile == "not found") { relabundfile = ""; }
152 else { inputFileName = relabundfile; }
154 metadatafile = validParameter.validFile(parameters, "metadata", true);
155 if (metadatafile == "not open") { abort = true; }
156 else if (metadatafile == "not found") { metadatafile = ""; }
157 else { inputFileName = metadatafile; }
159 groups = validParameter.validFile(parameters, "groups", false);
160 if (groups == "not found") { groups = ""; pickedGroups = false; }
163 m->splitAtDash(groups, Groups);
165 globaldata->Groups = Groups;
167 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(inputFileName); }
169 label = validParameter.validFile(parameters, "label", false);
170 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=""; }
172 if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) { m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true; }
174 if (metadatafile != "") {
175 if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
177 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true; }
180 temp = validParameter.validFile(parameters, "numaxes", false); if (temp == "not found"){ temp = "3"; }
181 convert(temp, numaxes);
183 method = validParameter.validFile(parameters, "method", false); if (method == "not found"){ method = "pearson"; }
185 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; }
189 catch(exception& e) {
190 m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
194 //**********************************************************************************************************************
196 void CorrAxesCommand::help(){
198 m->mothurOut("The corr.axes command reads a shared, relabund or metadata file as well as an axes file and calculates the correlation coefficient.\n");
199 m->mothurOut("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");
200 m->mothurOut("The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes.\n");
201 m->mothurOut("The label parameter allows you to select what distance level you would like used, if none is given the first distance is used.\n");
202 m->mothurOut("The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n");
203 m->mothurOut("The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n");
204 m->mothurOut("The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n");
205 m->mothurOut("Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n");
206 m->mothurOut("The corr.axes command outputs a .corr.axes file.\n");
207 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
209 catch(exception& e) {
210 m->errorOut(e, "CorrAxesCommand", "help");
215 //**********************************************************************************************************************
217 CorrAxesCommand::~CorrAxesCommand(){}
219 //**********************************************************************************************************************
221 int CorrAxesCommand::execute(){
224 if (abort == true) { return 0; }
226 /*************************************************************************************/
227 // use smart distancing to get right sharedRabund and convert to relabund if needed //
228 /************************************************************************************/
229 if (sharedfile != "") {
230 InputData* input = new InputData(sharedfile, "sharedfile");
231 getSharedFloat(input);
234 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
235 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
237 }else if (relabundfile != "") {
238 InputData* input = new InputData(relabundfile, "relabund");
239 getSharedFloat(input);
242 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
243 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
245 }else if (metadatafile != "") {
246 getMetadata(); //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
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 metadata file."); m->mothurOutEndLine(); return 0; }
250 if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
252 }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
254 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
256 //this is for a sanity check to make sure the axes file and shared file match
257 for (int i = 0; i < lookupFloat.size(); i++) { names.insert(lookupFloat[i]->getGroup()); }
259 /*************************************************************************************/
260 // read axes file and check for file mismatches with shared or relabund file //
261 /************************************************************************************/
264 map<string, vector<float> > axes = readAxes();
266 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
268 //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
269 if (axes.size() != lookupFloat.size()) {
270 map<string, vector<float> >::iterator it;
271 for (int i = 0; i < lookupFloat.size(); i++) {
272 it = axes.find(lookupFloat[i]->getGroup());
273 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(); }
275 m->control_pressed = true;
278 if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
280 /*************************************************************************************/
281 // calc the r values //
282 /************************************************************************************/
284 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
285 outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
287 m->openOutputFile(outputFileName, out);
288 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
291 if (metadatafile == "") { out << "OTU"; }
292 else { out << "Feature"; }
294 for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1); }
295 out << "\tlength" << endl;
297 if (method == "pearson") { calcPearson(axes, out); }
298 else if (method == "spearman") { calcSpearman(axes, out); }
299 else if (method == "kendall") { calcKendall(axes, out); }
300 else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
303 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
305 if (m->control_pressed) { return 0; }
307 m->mothurOutEndLine();
308 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
309 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
310 m->mothurOutEndLine();
314 catch(exception& e) {
315 m->errorOut(e, "CorrAxesCommand", "execute");
319 //**********************************************************************************************************************
320 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
323 //find average of each axis - X
324 vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
325 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
326 vector<float> temp = it->second;
327 for (int i = 0; i < temp.size(); i++) {
328 averageAxes[i] += temp[i];
332 for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
335 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
337 if (metadatafile == "") { out << i+1; }
338 else { out << metadataLabels[i]; }
340 //find the averages this otu - Y
342 for (int j = 0; j < lookupFloat.size(); j++) {
343 sumOtu += lookupFloat[j]->getAbundance(i);
345 float Ybar = sumOtu / (float) lookupFloat.size();
347 vector<float> rValues(averageAxes.size());
349 //find r value for each axis
350 for (int k = 0; k < averageAxes.size(); k++) {
353 double numerator = 0.0;
354 double denomTerm1 = 0.0;
355 double denomTerm2 = 0.0;
356 for (int j = 0; j < lookupFloat.size(); j++) {
357 float Yi = lookupFloat[j]->getAbundance(i);
358 float Xi = axes[lookupFloat[j]->getGroup()][k];
360 numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
361 denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
362 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
365 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
367 r = numerator / denom;
373 for(int k=0;k<rValues.size();k++){
374 sum += rValues[k] * rValues[k];
376 out << '\t' << sqrt(sum) << endl;
381 catch(exception& e) {
382 m->errorOut(e, "CorrAxesCommand", "calcPearson");
386 //**********************************************************************************************************************
387 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
391 axes["1a"].push_back(1);
392 axes["1b"].push_back(1);
393 axes["1c"].push_back(1);
394 axes["2a"].push_back(2);
395 axes["2b"].push_back(2);
396 axes["2c"].push_back(2);
397 axes["3a"].push_back(3);
398 axes["3b"].push_back(3);
399 axes["3c"].push_back(3);
400 axes["4a"].push_back(4);
401 axes["4b"].push_back(4);
402 axes["4c"].push_back(4);
403 axes["5a"].push_back(5);
404 axes["5b"].push_back(5);
405 axes["5c"].push_back(5);*/
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;
446 for (int j = 0; j < scores[i].size(); j++) {
448 ties.push_back(scores[i][j]);
450 if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
451 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
453 for (int k = 0; k < ties.size(); k++) {
454 float thisrank = rankTotal / (float) ties.size();
455 rankAxes[ties[k].name].push_back(thisrank);
460 }else { // you are the last one
462 for (int k = 0; k < ties.size(); k++) {
463 float thisrank = rankTotal / (float) ties.size();
464 rankAxes[ties[k].name].push_back(thisrank);
473 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
475 if (metadatafile == "") { out << i+1; }
476 else { out << metadataLabels[i]; }
478 //find the ranks of this otu - Y
479 vector<spearmanRank> otuScores;
480 map<float, int> tableY;
481 for (int j = 0; j < lookupFloat.size(); j++) {
482 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
483 otuScores.push_back(member);
485 itTable = tableY.find(member.score);
486 if (itTable == tableY.end()) {
487 tableY[member.score] = 1;
489 tableY[member.score]++;
496 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
497 double ty = (double) itTable->second;
498 Ly += ((pow(ty, 3.0) - ty) / 12.0);
501 sort(otuScores.begin(), otuScores.end(), compareSpearman);
503 map<string, float> rankOtus;
504 vector<spearmanRank> ties;
506 for (int j = 0; j < otuScores.size(); j++) {
508 ties.push_back(otuScores[j]);
510 if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
511 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
513 for (int k = 0; k < ties.size(); k++) {
514 float thisrank = rankTotal / (float) ties.size();
515 rankOtus[ties[k].name] = thisrank;
520 }else { // you are the last one
522 for (int k = 0; k < ties.size(); k++) {
523 float thisrank = rankTotal / (float) ties.size();
524 rankOtus[ties[k].name] = thisrank;
528 vector<double> pValues(numaxes);
530 //calc spearman ranks for each axis for this otu
531 for (int j = 0; j < numaxes; j++) {
534 for (int k = 0; k < lookupFloat.size(); k++) {
536 float xi = rankAxes[lookupFloat[k]->getGroup()][j];
537 float yi = rankOtus[lookupFloat[k]->getGroup()];
539 di += ((xi - yi) * (xi - yi));
544 double n = (double) lookupFloat.size();
546 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
547 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
549 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
557 for(int k=0;k<numaxes;k++){
558 sum += pValues[k] * pValues[k];
560 out << '\t' << sqrt(sum) << endl;
565 catch(exception& e) {
566 m->errorOut(e, "CorrAxesCommand", "calcSpearman");
570 //**********************************************************************************************************************
571 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
576 axes["1a"].push_back(1);
577 axes["1b"].push_back(1);
578 axes["1c"].push_back(1);
579 axes["2a"].push_back(2);
580 axes["2b"].push_back(2);
581 axes["2c"].push_back(2);
582 axes["3a"].push_back(3);
583 axes["3b"].push_back(3);
584 axes["3c"].push_back(3);
585 axes["4a"].push_back(4);
586 axes["4b"].push_back(4);
587 axes["4c"].push_back(4);
588 axes["5a"].push_back(5);
589 axes["5b"].push_back(5);
590 axes["5c"].push_back(5);
594 vector< vector<spearmanRank> > scores; scores.resize(numaxes);
595 for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
596 vector<float> temp = it->second;
597 for (int i = 0; i < temp.size(); i++) {
598 spearmanRank member(it->first, temp[i]);
599 scores[i].push_back(member);
604 for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
606 //convert scores to ranks of xi in each axis
607 for (int i = 0; i < numaxes; i++) {
609 vector<spearmanRank*> ties;
611 for (int j = 0; j < scores[i].size(); j++) {
613 ties.push_back(&(scores[i][j]));
615 if (j != scores[i].size()-1) { // you are not the last so you can look ahead
616 if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
617 for (int k = 0; k < ties.size(); k++) {
618 float thisrank = rankTotal / (float) ties.size();
619 (*ties[k]).score = thisrank;
624 }else { // you are the last one
625 for (int k = 0; k < ties.size(); k++) {
626 float thisrank = rankTotal / (float) ties.size();
627 (*ties[k]).score = thisrank;
632 cout << "axes ranks = ";
633 for(int i = 0; i < scores[0].size(); i++) { cout << scores[0][i].score << '\t'; }
636 lookupFloat.resize(15);
637 for (int i = 0; i < lookupFloat.size(); i++) {
638 lookupFloat[i] = new SharedRAbundFloatVector();
640 lookupFloat[0]->push_back(0.2288227, "1a"); lookupFloat[0]->setGroup("1a");
641 lookupFloat[1]->push_back(0.7394062, "1b"); lookupFloat[1]->setGroup("1b");
642 lookupFloat[2]->push_back(0.4521187, "1c"); lookupFloat[2]->setGroup("1c");
643 lookupFloat[3]->push_back(0.1598630, "2a"); lookupFloat[3]->setGroup("2a");
644 lookupFloat[4]->push_back(0.09588156, "2b"); lookupFloat[4]->setGroup("2b");
646 lookupFloat[5]->push_back(0.933174, "2c"); lookupFloat[5]->setGroup("2c");
647 lookupFloat[6]->push_back(0.3958304, "3a"); lookupFloat[6]->setGroup("3a");;
648 lookupFloat[7]->push_back(0.2364419, "3b"); lookupFloat[7]->setGroup("3b");
649 lookupFloat[8]->push_back(0.1697712, "3c"); lookupFloat[8]->setGroup("3c");
650 lookupFloat[9]->push_back(0.4077173, "4a"); lookupFloat[9]->setGroup("4a");
652 lookupFloat[10]->push_back(0.6116547, "4b"); lookupFloat[10]->setGroup("4b");
653 lookupFloat[11]->push_back(0.9374322, "4c"); lookupFloat[11]->setGroup("4c");
654 lookupFloat[12]->push_back(0.852184, "5a"); lookupFloat[12]->setGroup("5a");
655 lookupFloat[13]->push_back(0.845094, "5b"); lookupFloat[13]->setGroup("5b");
656 lookupFloat[14]->push_back(0.5795778, "5c"); lookupFloat[14]->setGroup("5c");
660 for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
662 if (metadatafile == "") { out << i+1; }
663 else { out << metadataLabels[i]; }
665 //find the ranks of this otu - Y
666 vector<spearmanRank> otuScores;
667 for (int j = 0; j < lookupFloat.size(); j++) {
668 spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
669 otuScores.push_back(member);
672 sort(otuScores.begin(), otuScores.end(), compareSpearman);
674 map<string, float> rankOtus;
675 vector<spearmanRank> ties;
677 for (int j = 0; j < otuScores.size(); j++) {
679 ties.push_back(otuScores[j]);
681 if (j != otuScores.size()-1) { // you are not the last so you can look ahead
682 if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
683 for (int k = 0; k < ties.size(); k++) {
684 float thisrank = rankTotal / (float) ties.size();
685 rankOtus[ties[k].name] = thisrank;
690 }else { // you are the last one
691 for (int k = 0; k < ties.size(); k++) {
692 float thisrank = rankTotal / (float) ties.size();
693 rankOtus[ties[k].name] = thisrank;
699 vector<double> pValues(numaxes);
701 //calc spearman ranks for each axis for this otu
702 for (int j = 0; j < numaxes; j++) {
707 //assemble otus ranks in same order as axis ranks, if ties, put in the best order possible
708 //NOTE: after this ordering the scores[j] indexes may not match the otus indexes within tied sections
709 //since we do not use the axes ranks except to order the otu ranks, I did not take the time to reorder the axes ranks
710 vector<spearmanRank> otus;
711 vector<spearmanRank> otusTemp;
712 for (int l = 0; l < scores[j].size(); l++) {
713 spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]);
714 otusTemp.push_back(member);
716 if (l != (scores[j].size()-1)) { // you are not the last so you can look ahead
717 if (scores[j][l].score != scores[j][l+1].score) { // you are done with ties, order them and continue
718 //order otus within tied section in the best way possible to make coor pairs
719 sort(otusTemp.begin(), otusTemp.end(), compareSpearmanReverse);
722 for (int h = 0; h < otusTemp.size(); h++) { ; otus.push_back(otusTemp[h]); }
726 }else { // you are the last one
727 //order otus within tied section in the best way possible to make coor pairs
728 sort(otusTemp.begin(), otusTemp.end(), compareSpearmanReverse);
731 for (int h = 0; h < otusTemp.size(); h++) { otus.push_back(otusTemp[h]); }
735 cout << "otu ranks = ";
736 for (int h = 0; h < otus.size(); h++ ) { cout << otus[h].score << '\t'; }
740 for (int l = 0; l < scores[j].size(); l++) {
742 int numWithHigherRank = 0;
743 int numWithLowerRank = 0;
744 float thisrank = otus[l].score;
746 for (int u = l; u < scores[j].size(); u++) {
747 if (otus[u].score > thisrank) { numWithHigherRank++; }
748 else if (otus[u].score < thisrank) { numWithLowerRank++; }
752 numCoor += numWithHigherRank;
753 numDisCoor += numWithLowerRank;
757 //int n = lookupFloat.size();
758 //comparing to yourself
759 count -= lookupFloat.size();
761 //double p = ( (4 * P) / (float) ((n * (n - 1)) - numTies) ) - 1.0;
762 double p = (numCoor - numDisCoor) / (float) count;
763 cout << "numCoor = " << numCoor << " numDisCoor = " << numDisCoor << " p = " << p << " count = " << count << endl;
770 for(int k=0;k<numaxes;k++){
771 sum += pValues[k] * pValues[k];
773 out << '\t' << sqrt(sum) << endl;
778 catch(exception& e) {
779 m->errorOut(e, "CorrAxesCommand", "calcKendall");
783 //**********************************************************************************************************************
784 int CorrAxesCommand::getSharedFloat(InputData* input){
786 lookupFloat = input->getSharedRAbundFloatVectors();
787 string lastLabel = lookupFloat[0]->getLabel();
789 if (label == "") { label = lastLabel; return 0; }
791 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
792 set<string> labels; labels.insert(label);
793 set<string> processedLabels;
794 set<string> userLabels = labels;
796 //as long as you are not at the end of the file or done wih the lines you want
797 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
799 if (m->control_pressed) { return 0; }
801 if(labels.count(lookupFloat[0]->getLabel()) == 1){
802 processedLabels.insert(lookupFloat[0]->getLabel());
803 userLabels.erase(lookupFloat[0]->getLabel());
807 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
808 string saveLabel = lookupFloat[0]->getLabel();
810 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
811 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
813 processedLabels.insert(lookupFloat[0]->getLabel());
814 userLabels.erase(lookupFloat[0]->getLabel());
816 //restore real lastlabel to save below
817 lookupFloat[0]->setLabel(saveLabel);
821 lastLabel = lookupFloat[0]->getLabel();
823 //get next line to process
824 //prevent memory leak
825 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
826 lookupFloat = input->getSharedRAbundFloatVectors();
830 if (m->control_pressed) { return 0; }
832 //output error messages about any remaining user labels
833 set<string>::iterator it;
834 bool needToRun = false;
835 for (it = userLabels.begin(); it != userLabels.end(); it++) {
836 m->mothurOut("Your file does not include the label " + *it);
837 if (processedLabels.count(lastLabel) != 1) {
838 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
841 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
845 //run last label if you need to
846 if (needToRun == true) {
847 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
848 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
853 catch(exception& e) {
854 m->errorOut(e, "CorrAxesCommand", "getSharedFloat");
858 //**********************************************************************************************************************
859 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
862 vector<SharedRAbundFloatVector*> newLookup;
863 for (int i = 0; i < thislookup.size(); i++) {
864 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
865 temp->setLabel(thislookup[i]->getLabel());
866 temp->setGroup(thislookup[i]->getGroup());
867 newLookup.push_back(temp);
871 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
872 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
874 //look at each sharedRabund and make sure they are not all zero
876 for (int j = 0; j < thislookup.size(); j++) {
877 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
880 //if they are not all zero add this bin
882 for (int j = 0; j < thislookup.size(); j++) {
883 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
888 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
890 thislookup = newLookup;
895 catch(exception& e) {
896 m->errorOut(e, "CorrAxesCommand", "eliminateZeroOTUS");
900 /*****************************************************************/
901 map<string, vector<float> > CorrAxesCommand::readAxes(){
903 map<string, vector<float> > axes;
906 m->openInputFile(axesfile, in);
908 string headerLine = m->getline(in); m->gobble(in);
910 //count the number of axis you are reading
914 int pos = headerLine.find("axis");
915 if (pos != string::npos) {
917 headerLine = headerLine.substr(pos+4);
918 }else { done = true; }
921 if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
925 if (m->control_pressed) { in.close(); return axes; }
928 in >> group; m->gobble(in);
930 vector<float> thisGroupsAxes;
931 for (int i = 0; i < count; i++) {
935 //only save the axis we want
936 if (i < numaxes) { thisGroupsAxes.push_back(temp); }
939 //save group if its one the user selected
940 if (names.count(group) != 0) {
941 map<string, vector<float> >::iterator it = axes.find(group);
943 if (it == axes.end()) {
944 axes[group] = thisGroupsAxes;
946 m->mothurOut(group + " is already in your axes file, using first definition."); m->mothurOutEndLine();
956 catch(exception& e) {
957 m->errorOut(e, "CorrAxesCommand", "readAxes");
961 /*****************************************************************/
962 int CorrAxesCommand::getMetadata(){
964 vector<string> groupNames;
967 m->openInputFile(axesfile, in);
969 string headerLine = m->getline(in); m->gobble(in);
970 istringstream iss (headerLine,istringstream::in);
972 //read the first label, because it refers to the groups
974 iss >> columnLabel; m->gobble(iss);
976 //save names of columns you are reading
978 iss >> columnLabel; m->gobble(iss);
979 metadataLabels.push_back(columnLabel);
981 int count = metadataLabels.size();
986 if (m->control_pressed) { in.close(); return 0; }
989 in >> group; m->gobble(in);
990 groupNames.push_back(group);
992 SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
993 tempLookup->setGroup(group);
994 tempLookup->setLabel("1");
996 for (int i = 0; i < count; i++) {
1000 tempLookup->push_back(temp, group);
1003 lookupFloat.push_back(tempLookup);
1009 //remove any groups the user does not want, and set globaldata->groups with only valid groups
1011 util = new SharedUtil();
1013 util->setGroups(globaldata->Groups, groupNames);
1015 for (int i = 0; i < lookupFloat.size(); i++) {
1016 //if this sharedrabund is not from a group the user wants then delete it.
1017 if (util->isValidGroup(lookupFloat[i]->getGroup(), globaldata->Groups) == false) {
1018 delete lookupFloat[i]; lookupFloat[i] = NULL;
1019 lookupFloat.erase(lookupFloat.begin()+i);
1028 catch(exception& e) {
1029 m->errorOut(e, "CorrAxesCommand", "getMetadata");
1033 /*****************************************************************/