X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=corraxescommand.cpp;h=f47427d62c2cc7f03615149b11cbfdee5135b083;hb=002421a439168e2610a2b62f1318f21c7202fe6d;hp=db2e9891738cdf57b78dffe1a9c232f18e5de31c;hpb=fa08e82ec2dd48f73d051c210dad54a403308949;p=mothur.git diff --git a/corraxescommand.cpp b/corraxescommand.cpp index db2e989..f47427d 100644 --- a/corraxescommand.cpp +++ b/corraxescommand.cpp @@ -10,11 +10,6 @@ #include "corraxescommand.h" #include "sharedutilities.h" -//******************************************************************************************************************** -//sorts highes to lowest -inline bool compareSpearman(spearmanRank left, spearmanRank right){ - return (left.score > right.score); -} //********************************************************************************************************************** vector CorrAxesCommand::getValidParameters(){ try { @@ -42,8 +37,7 @@ vector CorrAxesCommand::getRequiredParameters(){ //********************************************************************************************************************** CorrAxesCommand::CorrAxesCommand(){ try { - abort = true; - //initialize outputTypes + abort = true; calledHelp = true; vector tempOutNames; outputTypes["corr.axes"] = tempOutNames; } @@ -67,11 +61,11 @@ vector CorrAxesCommand::getRequiredFiles(){ //********************************************************************************************************************** CorrAxesCommand::CorrAxesCommand(string option) { try { - abort = false; + abort = false; calledHelp = false; globaldata = GlobalData::getInstance(); //allow user to run help - if(option == "help") { help(); abort = true; } + if(option == "help") { help(); abort = true; calledHelp = true; } else { //valid paramters for this command @@ -216,7 +210,7 @@ CorrAxesCommand::~CorrAxesCommand(){} int CorrAxesCommand::execute(){ try { - if (abort == true) { return 0; } + if (abort == true) { if (calledHelp) { return 0; } return 2; } /*************************************************************************************/ // use smart distancing to get right sharedRabund and convert to relabund if needed // @@ -383,12 +377,32 @@ int CorrAxesCommand::calcSpearman(map >& 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); } } @@ -402,11 +416,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); @@ -415,16 +430,17 @@ 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++) { @@ -433,9 +449,25 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o //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); @@ -444,11 +476,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; @@ -457,14 +490,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); + vector pValues(numaxes); + //calc spearman ranks for each axis for this otu for (int j = 0; j < numaxes; j++) { @@ -477,12 +511,18 @@ 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; + + 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; - } double sum = 0; @@ -522,10 +562,10 @@ int CorrAxesCommand::calcKendall(map >& 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(); @@ -562,10 +602,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(); @@ -581,35 +621,42 @@ 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++) { + for (int u = l+1; 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; - + double p = (numCoor - numDisCoor) / (float) count; + out << '\t' << p; pValues[j] = p;