#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<string> CorrAxesCommand::getValidParameters(){
try {
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); }
//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++) {
}
float Ybar = sumOtu / (float) lookupFloat.size();
+ vector<float> rValues(averageAxes.size());
+
//find r value for each axis
for (int k = 0; k < averageAxes.size(); k++) {
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<rValues.size();k++){
+ sum += rValues[k] * rValues[k];
+ }
+ out << '\t' << sqrt(sum) << endl;
}
return 0;
try {
//format data
+ vector< map<float, int> > tableX; tableX.resize(numaxes);
+ map<float, int>::iterator itTable;
vector< vector<spearmanRank> > scores; scores.resize(numaxes);
for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
vector<float> 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<double> 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);
}
}
vector<spearmanRank> 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);
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<spearmanRank> otuScores;
+ map<float, int> 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);
vector<spearmanRank> 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;
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<double> pValues(numaxes);
+
//calc spearman ranks for each axis for this otu
for (int j = 0; j < numaxes; j++) {
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<numaxes;k++){
+ sum += pValues[k] * pValues[k];
+ }
+ out << '\t' << sqrt(sum) << endl;
}
return 0;
vector<spearmanRank*> 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();
//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<spearmanRank> otuScores;
vector<spearmanRank> 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();
}
}
+
+ vector<double> 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<spearmanRank> otus;
+ vector<spearmanRank> 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<numaxes;k++){
+ sum += pValues[k] * pValues[k];
+ }
+ out << '\t' << sqrt(sum) << endl;
}
return 0;
projects / mothur.git / blobdiff