X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=corraxescommand.cpp;h=26cdcda0178bcfdbfc8c89fbbbced8866edb796c;hb=5bb0453ee6ab8f7e700b5d0c61a84fdc73db4976;hp=8ebd22c50144bf60ee5c7ff0a6f608745707a045;hpb=e2fa4794a30265d840fd66ed2f24de0c4d96d393;p=mothur.git diff --git a/corraxescommand.cpp b/corraxescommand.cpp index 8ebd22c..26cdcda 100644 --- a/corraxescommand.cpp +++ b/corraxescommand.cpp @@ -11,10 +11,15 @@ #include "sharedutilities.h" //******************************************************************************************************************** -//sorts highes to lowest +//sorts highest to lowest inline bool compareSpearman(spearmanRank left, spearmanRank right){ return (left.score > right.score); } +//******************************************************************************************************************** +//sorts lowest to highest +inline bool compareSpearmanReverse(spearmanRank left, spearmanRank right){ + return (left.score < right.score); +} //********************************************************************************************************************** vector CorrAxesCommand::getValidParameters(){ try { @@ -283,11 +288,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 +334,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 +344,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 +365,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& o try { //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 +427,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 +441,44 @@ 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); + } } } } - //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 +487,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 +501,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 +522,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 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(); @@ -534,8 +597,8 @@ int CorrAxesCommand::calcKendall(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 ranks of this otu - Y vector otuScores; @@ -550,10 +613,10 @@ int CorrAxesCommand::calcKendall(map >& axes, ofstream& ou 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(); @@ -570,37 +633,54 @@ 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++) { - int P = 0; - //assemble otus ranks in same order as axis ranks + int numCoor = 0; + int numDisCoor = 0; + vector otus; + vector otusTemp; for (int l = 0; l < scores[j].size(); l++) { spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]); otus.push_back(member); } + int count = 0; for (int l = 0; l < scores[j].size(); l++) { int numWithHigherRank = 0; + int numWithLowerRank = 0; float thisrank = otus[l].score; for (int u = l; u < scores[j].size(); u++) { if (otus[u].score > thisrank) { numWithHigherRank++; } + else if (otus[u].score < thisrank) { numWithLowerRank++; } + count++; } - P += numWithHigherRank; + numCoor += numWithHigherRank; + numDisCoor += numWithLowerRank; } - int n = lookupFloat.size(); - - double p = ( (4 * P) / (float) (n * (n - 1)) ) - 1.0; + //comparing to yourself + count -= lookupFloat.size(); - out << p << '\t'; + double p = (numCoor - numDisCoor) / (float) count; + + out << '\t' << p; + pValues[j] = p; + } - out << endl; + double sum = 0; + for(int k=0;k