X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=corraxescommand.cpp;h=b6c1a705a39a891072f147997f942708fcc52272;hb=5d6d303e481489e226fdf8d6c5385b99b50718bc;hp=0d421b5aba7d6be8e8f9bf13ae00d39eecbac926;hpb=37eac2026d91179acda0494c4dcca22f176551b9;p=mothur.git diff --git a/corraxescommand.cpp b/corraxescommand.cpp index 0d421b5..b6c1a70 100644 --- a/corraxescommand.cpp +++ b/corraxescommand.cpp @@ -190,7 +190,7 @@ CorrAxesCommand::CorrAxesCommand(string option) { void CorrAxesCommand::help(){ try { - m->mothurOut("The corr.axes command reads a shared or relabund file as well as a pcoa file and calculates the correlation coefficient.\n"); + 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"); 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"); 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"); 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"); @@ -283,11 +283,11 @@ int CorrAxesCommand::execute(){ out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); //output headings - if (metadatafile == "") { out << "OTU\t"; } - else { out << "Feature\t"; } + if (metadatafile == "") { out << "OTU"; } + else { out << "Feature"; } - for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; } - out << endl; + for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1); } + out << "\tlength" << endl; if (method == "pearson") { calcPearson(axes, out); } else if (method == "spearman") { calcSpearman(axes, out); } @@ -329,9 +329,9 @@ int CorrAxesCommand::calcPearson(map >& axes, ofstream& ou //for each otu for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) { - if (metadatafile == "") { out << i+1 << '\t'; } - else { out << metadataLabels[i] << '\t'; } - + if (metadatafile == "") { out << i+1; } + else { out << metadataLabels[i]; } + //find the averages this otu - Y float sumOtu = 0.0; for (int j = 0; j < lookupFloat.size(); j++) { @@ -339,6 +339,8 @@ int CorrAxesCommand::calcPearson(map >& axes, ofstream& ou } float Ybar = sumOtu / (float) lookupFloat.size(); + vector rValues(averageAxes.size()); + //find r value for each axis for (int k = 0; k < averageAxes.size(); k++) { @@ -358,11 +360,15 @@ int CorrAxesCommand::calcPearson(map >& axes, ofstream& ou double denom = (sqrt(denomTerm1) * sqrt(denomTerm2)); r = numerator / denom; - - out << r << '\t'; + rValues[k] = r; + out << '\t' << r; } - out << endl; + double sum = 0; + for(int k=0;k >& axes, ofstream& ou int CorrAxesCommand::calcSpearman(map >& axes, ofstream& out) { try { + /*axes.clear(); + axes["1a"].push_back(1); + axes["1b"].push_back(1); + axes["1c"].push_back(1); + axes["2a"].push_back(2); + axes["2b"].push_back(2); + axes["2c"].push_back(2); + axes["3a"].push_back(3); + axes["3b"].push_back(3); + axes["3c"].push_back(3); + axes["4a"].push_back(4); + axes["4b"].push_back(4); + axes["4c"].push_back(4); + axes["5a"].push_back(5); + axes["5b"].push_back(5); + axes["5c"].push_back(5);*/ + //format data + vector< map > tableX; tableX.resize(numaxes); + map::iterator itTable; vector< vector > scores; scores.resize(numaxes); for (map >::iterator it = axes.begin(); it != axes.end(); it++) { vector temp = it->second; for (int i = 0; i < temp.size(); i++) { spearmanRank member(it->first, temp[i]); scores[i].push_back(member); + + //count number of repeats + itTable = tableX[i].find(temp[i]); + if (itTable == tableX[i].end()) { + tableX[i][temp[i]] = 1; + }else { + tableX[i][temp[i]]++; + } + } + } + + //calc LX + //for each axis + vector Lx; Lx.resize(numaxes, 0.0); + for (int i = 0; i < numaxes; i++) { + for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) { + double tx = (double) itTable->second; + Lx[i] += ((pow(tx, 3.0) - tx) / 12.0); } } @@ -396,11 +439,12 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o vector ties; int rankTotal = 0; for (int j = 0; j < scores[i].size(); j++) { - rankTotal += j; + rankTotal += (j+1); ties.push_back(scores[i][j]); - if (j != scores.size()-1) { // you are not the last so you can look ahead + if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue + for (int k = 0; k < ties.size(); k++) { float thisrank = rankTotal / (float) ties.size(); rankAxes[ties[k].name].push_back(thisrank); @@ -409,27 +453,67 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o rankTotal = 0; } }else { // you are the last one + for (int k = 0; k < ties.size(); k++) { float thisrank = rankTotal / (float) ties.size(); rankAxes[ties[k].name].push_back(thisrank); + } } } } - + /*cout << endl; + lookupFloat.clear(); + lookupFloat.resize(15); + for (int i = 0; i < lookupFloat.size(); i++) { + lookupFloat[i] = new SharedRAbundFloatVector(); + } + lookupFloat[0]->push_back(0.2288227, "1a"); lookupFloat[0]->setGroup("1a"); + lookupFloat[1]->push_back(0.7394062, "1b"); lookupFloat[1]->setGroup("1b"); + lookupFloat[2]->push_back(0.4521187, "1c"); lookupFloat[2]->setGroup("1c"); + lookupFloat[3]->push_back(0.1598630, "2a"); lookupFloat[3]->setGroup("2a"); + lookupFloat[4]->push_back(0.09588156, "2b"); lookupFloat[4]->setGroup("2b"); + + lookupFloat[5]->push_back(0.933174, "2c"); lookupFloat[5]->setGroup("2c"); + lookupFloat[6]->push_back(0.3958304, "3a"); lookupFloat[6]->setGroup("3a");; + lookupFloat[7]->push_back(0.2364419, "3b"); lookupFloat[7]->setGroup("3b"); + lookupFloat[8]->push_back(0.1697712, "3c"); lookupFloat[8]->setGroup("3c"); + lookupFloat[9]->push_back(0.4077173, "4a"); lookupFloat[9]->setGroup("4a"); + + lookupFloat[10]->push_back(0.6116547, "4b"); lookupFloat[10]->setGroup("4b"); + lookupFloat[11]->push_back(0.9374322, "4c"); lookupFloat[11]->setGroup("4c"); + lookupFloat[12]->push_back(0.852184, "5a"); lookupFloat[12]->setGroup("5a"); + lookupFloat[13]->push_back(0.845094, "5b"); lookupFloat[13]->setGroup("5b"); + lookupFloat[14]->push_back(0.5795778, "5c"); lookupFloat[14]->setGroup("5c");*/ //for each otu for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) { - if (metadatafile == "") { out << i+1 << '\t'; } - else { out << metadataLabels[i] << '\t'; } + if (metadatafile == "") { out << i+1; } + else { out << metadataLabels[i]; } //find the ranks of this otu - Y vector otuScores; + map tableY; for (int j = 0; j < lookupFloat.size(); j++) { spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i)); otuScores.push_back(member); + + itTable = tableY.find(member.score); + if (itTable == tableY.end()) { + tableY[member.score] = 1; + }else { + tableY[member.score]++; + } + + } + + //calc Ly + double Ly = 0.0; + for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) { + double ty = (double) itTable->second; + Ly += ((pow(ty, 3.0) - ty) / 12.0); } sort(otuScores.begin(), otuScores.end(), compareSpearman); @@ -438,11 +522,12 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o vector ties; int rankTotal = 0; for (int j = 0; j < otuScores.size(); j++) { - rankTotal += j; + rankTotal += (j+1); ties.push_back(otuScores[j]); - if (j != scores.size()-1) { // you are not the last so you can look ahead + if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue + for (int k = 0; k < ties.size(); k++) { float thisrank = rankTotal / (float) ties.size(); rankOtus[ties[k].name] = thisrank; @@ -451,13 +536,15 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o rankTotal = 0; } }else { // you are the last one + for (int k = 0; k < ties.size(); k++) { float thisrank = rankTotal / (float) ties.size(); rankOtus[ties[k].name] = thisrank; } } } - + vector pValues(numaxes); + //calc spearman ranks for each axis for this otu for (int j = 0; j < numaxes; j++) { @@ -470,14 +557,25 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o di += ((xi - yi) * (xi - yi)); } - int n = lookupFloat.size(); - double p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1))); + double p = 0.0; - out << p << '\t'; + double n = (double) lookupFloat.size(); + + double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j]; + double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly; + + p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2))); + + out << '\t' << p; + + pValues[j] = p; } - - - out << endl; + + double sum = 0; + for(int k=0;k >& axes, ofstream& ou //for each otu for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) { - if (metadatafile == "") { out << i+1 << '\t'; } - else { out << metadataLabels[i] << '\t'; } + if (metadatafile == "") { out << i+1; } + else { out << metadataLabels[i]; } //find the ranks of this otu - Y vector otuScores; @@ -569,6 +667,7 @@ int CorrAxesCommand::calcKendall(map >& axes, ofstream& ou } } } + vector pValues(numaxes); //calc spearman ranks for each axis for this otu for (int j = 0; j < numaxes; j++) { @@ -597,10 +696,16 @@ int CorrAxesCommand::calcKendall(map >& axes, ofstream& ou double p = ( (4 * P) / (float) (n * (n - 1)) ) - 1.0; - out << p << '\t'; + out << '\t' << p; + pValues[j] = p; + } - out << endl; + double sum = 0; + for(int k=0;k