added mantel command

[mothur.git] / corraxescommand.cpp

diff --git a/corraxescommand.cpp b/corraxescommand.cpp
index 3da6b82d940fc8257a5564fc1c5b9a419343f604..f47427d62c2cc7f03615149b11cbfdee5135b083 100644 (file)
--- 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<string> CorrAxesCommand::getValidParameters(){  
        try {
@@ -42,8 +37,7 @@ vector<string> CorrAxesCommand::getRequiredParameters(){
 //**********************************************************************************************************************
 CorrAxesCommand::CorrAxesCommand(){    
        try {
-               abort = true;
-               //initialize outputTypes
+               abort = true; calledHelp = true; 
                vector<string> tempOutNames;
                outputTypes["corr.axes"] = tempOutNames;
        }
@@ -67,11 +61,11 @@ vector<string> 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  //
@@ -382,23 +376,6 @@ int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& ou
 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& 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<float, int> > tableX; tableX.resize(numaxes);
                map<float, int>::iterator itTable;
@@ -463,30 +440,7 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                        }
                }
                
-               /*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++) {
                        
@@ -667,35 +621,42 @@ 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++) {
                        
-                               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++) {
+                                       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;