added mantel command

[mothur.git] / corraxescommand.cpp

diff --git a/corraxescommand.cpp b/corraxescommand.cpp
index 0d421b5aba7d6be8e8f9bf13ae00d39eecbac926..f47427d62c2cc7f03615149b11cbfdee5135b083 100644 (file)
--- a/corraxescommand.cpp
+++ b/corraxescommand.cpp
@@ -10,11 +10,6 @@
 #include "corraxescommand.h"
 #include "sharedutilities.h"
 
 #include "corraxescommand.h"
 #include "sharedutilities.h"
 

-//********************************************************************************************************************
-//sorts highes to lowest
-inline bool compareSpearman(spearmanRank left, spearmanRank right){
-       return (left.score > right.score);      
-} 

 //**********************************************************************************************************************
 vector<string> CorrAxesCommand::getValidParameters(){  
        try {
 //**********************************************************************************************************************
 vector<string> CorrAxesCommand::getValidParameters(){  
        try {


@@ -42,8 +37,7 @@ vector<string> CorrAxesCommand::getRequiredParameters(){
 //**********************************************************************************************************************
 CorrAxesCommand::CorrAxesCommand(){    
        try {
 //**********************************************************************************************************************
 CorrAxesCommand::CorrAxesCommand(){    
        try {

-               abort = true;
-               //initialize outputTypes
+               abort = true; calledHelp = true; 

                vector<string> tempOutNames;
                outputTypes["corr.axes"] = tempOutNames;
        }
                vector<string> tempOutNames;
                outputTypes["corr.axes"] = tempOutNames;
        }


@@ -67,11 +61,11 @@ vector<string> CorrAxesCommand::getRequiredFiles(){
 //**********************************************************************************************************************
 CorrAxesCommand::CorrAxesCommand(string option)  {
        try {
 //**********************************************************************************************************************
 CorrAxesCommand::CorrAxesCommand(string option)  {
        try {

-               abort = false;
+               abort = false; calledHelp = false;   

                globaldata = GlobalData::getInstance();
                
                //allow user to run help
                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
                
                else {
                        //valid paramters for this command


@@ -190,7 +184,7 @@ CorrAxesCommand::CorrAxesCommand(string option)  {
 
 void CorrAxesCommand::help(){
        try {
 
 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");
                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");


@@ -216,7 +210,7 @@ CorrAxesCommand::~CorrAxesCommand(){}
 int CorrAxesCommand::execute(){
        try {
                
 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  //
                
                /*************************************************************************************/
                // use smart distancing to get right sharedRabund and convert to relabund if needed  //


@@ -283,11 +277,11 @@ int CorrAxesCommand::execute(){
                out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
                
                //output headings
                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); }
                
                if (method == "pearson")                {  calcPearson(axes, out);      }
                else if (method == "spearman")  {  calcSpearman(axes, out); }


@@ -329,9 +323,9 @@ int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& ou
           //for each otu
           for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
                   
           //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++) {
                   //find the averages this otu - Y
                   float sumOtu = 0.0;
                   for (int j = 0; j < lookupFloat.size(); j++) {


@@ -339,6 +333,8 @@ int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& ou
                   }
                   float Ybar = sumOtu / (float) lookupFloat.size();
                   
                   }
                   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++) {
                           
                   //find r value for each axis
                   for (int k = 0; k < averageAxes.size(); k++) {
                           


@@ -358,11 +354,15 @@ int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& ou
                           double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
                           
                           r = numerator / denom;
                           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;
           }
                   
           return 0;


@@ -377,12 +377,32 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
        try {
                
                //format data
        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);  
                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);

                        }
                }
                
                        }
                }
                


@@ -396,11 +416,12 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                        vector<spearmanRank> ties;
                        int rankTotal = 0;
                        for (int j = 0; j < scores[i].size(); j++) {
                        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]);
                                
                                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
                                        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);
                                                for (int k = 0; k < ties.size(); k++) {
                                                        float thisrank = rankTotal / (float) ties.size();
                                                        rankAxes[ties[k].name].push_back(thisrank);


@@ -409,27 +430,44 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                                                rankTotal = 0;
                                        }
                                }else { // you are the last one
                                                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 (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++) {
                        
                //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;
                        
                        //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);
                        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);
                        }
                        
                        sort(otuScores.begin(), otuScores.end(), compareSpearman);


@@ -438,11 +476,12 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                        vector<spearmanRank> ties;
                        int rankTotal = 0;
                        for (int j = 0; j < otuScores.size(); j++) {
                        vector<spearmanRank> ties;
                        int rankTotal = 0;
                        for (int j = 0; j < otuScores.size(); j++) {

-                               rankTotal += j;
+                               rankTotal += (j+1);

                                ties.push_back(otuScores[j]);
                                
                                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
                                        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;
                                                for (int k = 0; k < ties.size(); k++) {
                                                        float thisrank = rankTotal / (float) ties.size();
                                                        rankOtus[ties[k].name] = thisrank;


@@ -451,13 +490,15 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                                                rankTotal = 0;
                                        }
                                }else { // you are the last one
                                                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;
                                        }
                                }
                        }
                                        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++) {
                                
                        //calc spearman ranks for each axis for this otu
                        for (int j = 0; j < numaxes; j++) {
                                


@@ -470,14 +511,25 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                                        di += ((xi - yi) * (xi - yi));
                                }
                                
                                        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();

                                
                                

-                               out << p << '\t';
+                               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;
                }
                
                return 0;


@@ -510,10 +562,10 @@ int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& ou
                        vector<spearmanRank*> ties;
                        int rankTotal = 0;
                        for (int j = 0; j < scores[i].size(); j++) {
                        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]));
                                
                                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();
                                        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 +586,8 @@ int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& ou
                //for each otu
                for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
                
                //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;
                        
                        //find the ranks of this otu - Y
                        vector<spearmanRank> otuScores;


@@ -550,10 +602,10 @@ int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& ou
                        vector<spearmanRank> ties;
                        int rankTotal = 0;
                        for (int j = 0; j < otuScores.size(); j++) {
                        vector<spearmanRank> ties;
                        int rankTotal = 0;
                        for (int j = 0; j < otuScores.size(); j++) {

-                               rankTotal += j;
+                               rankTotal += (j+1);

                                ties.push_back(otuScores[j]);
                                
                                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();
                                        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 +622,51 @@ int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& ou
                                }
                        }
                        
                                }
                        }
                        

+                       
+                       vector<double> pValues(numaxes);
+                       

                        //calc spearman ranks for each axis for this otu
                        for (int j = 0; j < numaxes; j++) {
                        
                        //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> 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);
                                }
                                
                                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;
                                for (int l = 0; l < scores[j].size(); l++) {
                                        
                                        int numWithHigherRank = 0;

+                                       int numWithLowerRank = 0;

                                        float thisrank = otus[l].score;
                                        
                                        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++; }
                                                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;
-                               
-                               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;
                }
                
                return 0;