fixes while testing 1.33.0

[mothur.git] / corraxescommand.cpp

diff --git a/corraxescommand.cpp b/corraxescommand.cpp
index eac1f3505b1bca098ab01d100a513d24bd5c9d29..72fa03b75edc196aafcaaa4b6ade6f95d1106d48 100644 (file)
--- a/corraxescommand.cpp
+++ b/corraxescommand.cpp
@@ -9,74 +9,92 @@
 
 #include "corraxescommand.h"
 #include "sharedutilities.h"
 
 #include "corraxescommand.h"
 #include "sharedutilities.h"

+#include "linearalgebra.h"

 
 

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

 //**********************************************************************************************************************
 //**********************************************************************************************************************

-vector<string> CorrAxesCommand::getValidParameters(){  
+vector<string> CorrAxesCommand::setParameters(){       

        try {
        try {

-               string Array[] =  {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
-               vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+               CommandParameter paxes("axes", "InputTypes", "", "", "none", "none", "none","corraxes",false,true,true); parameters.push_back(paxes);
+               CommandParameter pshared("shared", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none","",false,false,true); parameters.push_back(pshared);
+               CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none","",false,false,true); parameters.push_back(prelabund);
+               CommandParameter pmetadata("metadata", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none","",false,false); parameters.push_back(pmetadata);
+               CommandParameter pnumaxes("numaxes", "Number", "", "3", "", "", "","",false,false); parameters.push_back(pnumaxes);
+               CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
+               CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
+               CommandParameter pmethod("method", "Multiple", "pearson-spearman-kendall", "pearson", "", "", "","",false,false); parameters.push_back(pmethod);
+               CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
+               CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
+               
+               vector<string> myArray;
+               for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }

                return myArray;
        }
        catch(exception& e) {
                return myArray;
        }
        catch(exception& e) {

-               m->errorOut(e, "CorrAxesCommand", "getValidParameters");
+               m->errorOut(e, "CorrAxesCommand", "setParameters");

                exit(1);
        }
 }
 //**********************************************************************************************************************
                exit(1);
        }
 }
 //**********************************************************************************************************************

-vector<string> CorrAxesCommand::getRequiredParameters(){       
+string CorrAxesCommand::getHelpString(){       

        try {
        try {

-               string Array[] =  {"axes"};
-               vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
-               return myArray;
+               string helpString = "";
+               helpString += "The corr.axes command reads a shared, relabund or metadata file as well as an axes file and calculates the correlation coefficient.\n";
+               helpString += "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";
+               helpString += "The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes.\n";
+               helpString += "The label parameter allows you to select what distance level you would like used, if none is given the first distance is used.\n";
+               helpString += "The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n";
+               helpString += "The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n";
+               helpString += "The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n";
+               helpString += "Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n";
+               helpString += "The corr.axes command outputs a .corr.axes file.\n";
+               helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
+               return helpString;

        }
        catch(exception& e) {
        }
        catch(exception& e) {

-               m->errorOut(e, "CorrAxesCommand", "getRequiredParameters");
+               m->errorOut(e, "CorrAxesCommand", "getHelpString");

                exit(1);
        }
 }
 //**********************************************************************************************************************
                exit(1);
        }
 }
 //**********************************************************************************************************************

-CorrAxesCommand::CorrAxesCommand(){    
-       try {
-               abort = true;
-               //initialize outputTypes
-               vector<string> tempOutNames;
-               outputTypes["corr.axes"] = tempOutNames;
-       }
-       catch(exception& e) {
-               m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");
-               exit(1);
-       }
+string CorrAxesCommand::getOutputPattern(string type) {
+    try {
+        string pattern = "";
+        
+        if (type == "corraxes") {  pattern = "[filename],[tag],corr.axes"; }
+        else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
+        
+        return pattern;
+    }
+    catch(exception& e) {
+        m->errorOut(e, "CorrAxesCommand", "getOutputPattern");
+        exit(1);
+    }

 }
 
 //**********************************************************************************************************************
 }
 
 //**********************************************************************************************************************

-vector<string> CorrAxesCommand::getRequiredFiles(){    
+CorrAxesCommand::CorrAxesCommand(){    

        try {
        try {

-               vector<string> myArray;
-               return myArray;
+               abort = true; calledHelp = true; 
+               setParameters();
+               vector<string> tempOutNames;
+               outputTypes["corraxes"] = tempOutNames;

        }
        catch(exception& e) {
        }
        catch(exception& e) {

-               m->errorOut(e, "CorrAxesCommand", "getRequiredFiles");
+               m->errorOut(e, "CorrAxesCommand", "CorrAxesCommand");

                exit(1);
        }
 }
 //**********************************************************************************************************************
 CorrAxesCommand::CorrAxesCommand(string option)  {
        try {
                exit(1);
        }
 }
 //**********************************************************************************************************************
 CorrAxesCommand::CorrAxesCommand(string option)  {
        try {

-               abort = false;
-               globaldata = GlobalData::getInstance();
+               abort = false; calledHelp = false;   

                
                //allow user to run help
                
                //allow user to run help

-               if(option == "help") { help(); abort = true; }
+               if(option == "help") { help(); abort = true; calledHelp = true; }
+               else if(option == "citation") { citation(); abort = true; calledHelp = true;}

                
                else {
                
                else {

-                       //valid paramters for this command
-                       string Array[] =  {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
-                       vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+                       vector<string> myArray = setParameters();

                        
                        OptionParser parser(option);
                        map<string, string> parameters = parser.getParameters();
                        
                        OptionParser parser(option);
                        map<string, string> parameters = parser.getParameters();


@@ -90,7 +108,7 @@ CorrAxesCommand::CorrAxesCommand(string option)  {
                        }
                        
                        vector<string> tempOutNames;
                        }
                        
                        vector<string> tempOutNames;

-                       outputTypes["corr.axes"] = tempOutNames;
+                       outputTypes["corraxes"] = tempOutNames;

                        
                        //if the user changes the input directory command factory will send this info to us in the output parameter 
                        string inputDir = validParameter.validFile(parameters, "inputdir", false);              
                        
                        //if the user changes the input directory command factory will send this info to us in the output parameter 
                        string inputDir = validParameter.validFile(parameters, "inputdir", false);              


@@ -139,12 +157,12 @@ CorrAxesCommand::CorrAxesCommand(string option)  {
                        sharedfile = validParameter.validFile(parameters, "shared", true);
                        if (sharedfile == "not open") { abort = true; }
                        else if (sharedfile == "not found") { sharedfile = ""; }
                        sharedfile = validParameter.validFile(parameters, "shared", true);
                        if (sharedfile == "not open") { abort = true; }
                        else if (sharedfile == "not found") { sharedfile = ""; }

-                       else { inputFileName = sharedfile; }
+                       else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }

                        
                        relabundfile = validParameter.validFile(parameters, "relabund", true);
                        if (relabundfile == "not open") { abort = true; }
                        else if (relabundfile == "not found") { relabundfile = ""; }
                        
                        relabundfile = validParameter.validFile(parameters, "relabund", true);
                        if (relabundfile == "not open") { abort = true; }
                        else if (relabundfile == "not found") { relabundfile = ""; }

-                       else { inputFileName = relabundfile; }
+                       else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }

                        
                        metadatafile = validParameter.validFile(parameters, "metadata", true);
                        if (metadatafile == "not open") { abort = true; }
                        
                        metadatafile = validParameter.validFile(parameters, "metadata", true);
                        if (metadatafile == "not open") { abort = true; }


@@ -157,14 +175,27 @@ CorrAxesCommand::CorrAxesCommand(string option)  {
                                pickedGroups = true;
                                m->splitAtDash(groups, Groups); 
                        }                       
                                pickedGroups = true;
                                m->splitAtDash(groups, Groups); 
                        }                       

-                       globaldata->Groups = Groups;
+                       m->setGroups(Groups);

                        
                        outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(inputFileName);  }
                        
                        label = validParameter.validFile(parameters, "label", false);                   
                        if (label == "not found") { label = ""; m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); label=""; }  
                        
                        
                        outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(inputFileName);  }
                        
                        label = validParameter.validFile(parameters, "label", false);                   
                        if (label == "not found") { label = ""; m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); label=""; }  
                        

-                       if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) { m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true;  }
+                       if ((relabundfile == "") && (sharedfile == "") && (metadatafile == "")) { 
+                               //is there are current file available for any of these?
+                               //give priority to shared, then relabund
+                               //if there is a current shared file, use it
+                               sharedfile = m->getSharedFile(); 
+                               if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
+                               else { 
+                                       relabundfile = m->getRelAbundFile(); 
+                                       if (relabundfile != "") { inputFileName = relabundfile;  m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
+                                       else { 
+                                               m->mothurOut("You must provide either a shared, relabund, or metadata file."); m->mothurOutEndLine(); abort = true; 
+                                       }
+                               }
+                       }       

                        
                        if (metadatafile != "") {
                                if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true;  }
                        
                        if (metadatafile != "") {
                                if ((relabundfile != "") || (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true;  }


@@ -173,7 +204,7 @@ CorrAxesCommand::CorrAxesCommand(string option)  {
                        }
                        string temp;
                        temp = validParameter.validFile(parameters, "numaxes", false);          if (temp == "not found"){       temp = "3";                             }
                        }
                        string temp;
                        temp = validParameter.validFile(parameters, "numaxes", false);          if (temp == "not found"){       temp = "3";                             }

-                       convert(temp, numaxes); 
+                       m->mothurConvert(temp, numaxes); 

                        
                        method = validParameter.validFile(parameters, "method", false);         if (method == "not found"){     method = "pearson";             }
                        
                        
                        method = validParameter.validFile(parameters, "method", false);         if (method == "not found"){     method = "pearson";             }
                        


@@ -188,35 +219,10 @@ CorrAxesCommand::CorrAxesCommand(string option)  {
 }
 //**********************************************************************************************************************
 
 }
 //**********************************************************************************************************************
 

-void CorrAxesCommand::help(){
-       try {
-               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 method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n");
-               m->mothurOut("The numaxes parameter allows you to select the number of axes you would like to use. Default=3.\n");
-               m->mothurOut("The corr.axes command should be in the following format: corr.axes(axes=yourPcoaFile, shared=yourSharedFile, method=yourMethod).\n");
-               m->mothurOut("Example corr.axes(axes=genus.pool.thetayc.genus.lt.pcoa, shared=genus.pool.shared, method=kendall).\n");
-               m->mothurOut("The corr.axes command outputs a .corr.axes file.\n");
-               m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
-       }
-       catch(exception& e) {
-               m->errorOut(e, "CorrAxesCommand", "help");      
-               exit(1);
-       }
-}
-
-//**********************************************************************************************************************
-
-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  //


@@ -275,9 +281,11 @@ int CorrAxesCommand::execute(){
                /*************************************************************************************/
                // calc the r values                                                                                                                            //
                /************************************************************************************/
                /*************************************************************************************/
                // calc the r values                                                                                                                            //
                /************************************************************************************/

-               
-               string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
-               outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);      
+        map<string, string> variables; 
+        variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(inputFileName));
+        variables["[tag]"] = method;
+               string outputFileName = getOutputFileName("corraxes", variables);
+               outputNames.push_back(outputFileName); outputTypes["corraxes"].push_back(outputFileName);       

                ofstream out;
                m->openOutputFile(outputFileName, out);
                out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
                ofstream out;
                m->openOutputFile(outputFileName, out);
                out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);


@@ -286,7 +294,7 @@ int CorrAxesCommand::execute(){
                if (metadatafile == "") {  out << "OTU";        }
                else {  out << "Feature";                                               }
 
                if (metadatafile == "") {  out << "OTU";        }
                else {  out << "Feature";                                               }
 

-               for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1); }
+               for (int i = 0; i < numaxes; i++) { out << '\t' << "axis" << (i+1) << "\tp-value"; }

                out << "\tlength" << endl;
                
                if (method == "pearson")                {  calcPearson(axes, out);      }
                out << "\tlength" << endl;
                
                if (method == "pearson")                {  calcPearson(axes, out);      }


@@ -315,6 +323,8 @@ int CorrAxesCommand::execute(){
 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
    try {
           
 int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& out) {
    try {
           

+       LinearAlgebra linear;
+       

           //find average of each axis - X
           vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
           for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
           //find average of each axis - X
           vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
           for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {


@@ -329,7 +339,7 @@ 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;       }
+                  if (metadatafile == "") {  out << m->currentSharedBinLabels[i];      }

                   else {  out << metadataLabels[i];            }
                                   
                   //find the averages this otu - Y
                   else {  out << metadataLabels[i];            }
                                   
                   //find the averages this otu - Y


@@ -360,8 +370,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;

+               
+               if (isnan(r) || isinf(r)) { r = 0.0; }
+               

                           rValues[k] = r;
                           out << '\t' << r; 
                           rValues[k] = r;
                           out << '\t' << r; 

+               
+               double sig = linear.calcPearsonSig(lookupFloat.size(), r);
+               
+               out << '\t' << sig;

                   }
                   
                   double sum = 0;
                   }
                   
                   double sum = 0;


@@ -382,15 +399,36 @@ int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& ou
 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
        try {
                
 int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
        try {
                

-               bool hasTies = false;
-               
+        LinearAlgebra linear;
+        vector<double> sf; 
+        

                //format data
                //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);

                        }
                }
                
                        }
                }
                


@@ -398,124 +436,139 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                for (int i = 0; i < numaxes; i++) {  sort(scores[i].begin(), scores[i].end(), compareSpearman); }
                
                //find ranks of xi in each axis
                for (int i = 0; i < numaxes; i++) {  sort(scores[i].begin(), scores[i].end(), compareSpearman); }
                
                //find ranks of xi in each axis

-               vector<float> averageX; averageX.resize(numaxes, 0.0);

                map<string, vector<float> > rankAxes;
                for (int i = 0; i < numaxes; i++) {
                        
                        vector<spearmanRank> ties;
                        int rankTotal = 0;
                map<string, vector<float> > rankAxes;
                for (int i = 0; i < numaxes; i++) {
                        
                        vector<spearmanRank> ties;
                        int rankTotal = 0;

+            double sfTemp = 0.0;

                        for (int j = 0; j < scores[i].size(); j++) {
                        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

-                                               if (ties.size() > 1) { hasTies = true; }
-                                               averageX[i] += rankTotal;
+

                                                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);
                                                }

+                        int t = ties.size();
+                        sfTemp += (t*t*t-t);

                                                ties.clear();
                                                rankTotal = 0;
                                        }
                                }else { // you are the last one
                                                ties.clear();
                                                rankTotal = 0;
                                        }
                                }else { // you are the last one

-                                       if (ties.size() > 1) { hasTies = true; }
-                                       averageX[i] += rankTotal;
+                                       

                                        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);

+                                               

                                        }
                                }
                        }
                                        }
                                }
                        }

-                       
-                       averageX[i] /= (float) scores[i].size();
+            sf.push_back(sfTemp);

                }
                
                                
                }
                
                                

-               

                //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;  }
+                       if (metadatafile == "") {  out << m->currentSharedBinLabels[i]; }

                        else {  out << metadataLabels[i];               }
                        
                        //find the ranks of this otu - Y
                        vector<spearmanRank> otuScores;
                        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);
                        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);
                        

+            double sg = 0.0;

                        map<string, float> rankOtus;
                        vector<spearmanRank> ties;
                        map<string, float> rankOtus;
                        vector<spearmanRank> ties;

-                       float averageY = 0.0;

                        int rankTotal = 0;
                        for (int j = 0; j < otuScores.size(); j++) {
                        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

-                                               if (ties.size() > 1) { hasTies = true; }
-                                               averageY += rankTotal;
+                                               

                                                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;
                                                }

+                        int t = ties.size();
+                        sg += (t*t*t-t);

                                                ties.clear();
                                                rankTotal = 0;
                                        }
                                }else { // you are the last one
                                                ties.clear();
                                                rankTotal = 0;
                                        }
                                }else { // you are the last one

-                                       if (ties.size() > 1) { hasTies = true; }
-                                       averageY += rankTotal;
+                                       

                                        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;
                                        }
                                }
                        }

-                       
-
-                       averageY /= (float) otuScores.size();
-                       //hasTies = false;              
+                       vector<double> pValues(numaxes);        

 
                        //calc spearman ranks for each axis for this otu
                        for (int j = 0; j < numaxes; j++) {
                                
                                double di = 0.0;
 
                        //calc spearman ranks for each axis for this otu
                        for (int j = 0; j < numaxes; j++) {
                                
                                double di = 0.0;

-                               double denom1 = 0.0;
-                               double denom2 = 0.0;

                                for (int k = 0; k < lookupFloat.size(); k++) {
                                        
                                        float xi = rankAxes[lookupFloat[k]->getGroup()][j];
                                        float yi = rankOtus[lookupFloat[k]->getGroup()];
                                        
                                for (int k = 0; k < lookupFloat.size(); k++) {
                                        
                                        float xi = rankAxes[lookupFloat[k]->getGroup()][j];
                                        float yi = rankOtus[lookupFloat[k]->getGroup()];
                                        

-                                       if (hasTies) {
-                                               di += ((xi - averageX[j]) * (yi - averageY));
-                                               denom1 += ((xi - averageX[j]) * (xi - averageX[j]));
-                                               denom2 += ((yi - averageY) * (yi - averageY));
-                                       }else {
-                                               di += ((xi - yi) * (xi - yi));
-                                       }
+                                       di += ((xi - yi) * (xi - yi));

                                }
                                
                                double p = 0.0;
                                
                                }
                                
                                double p = 0.0;
                                

-                               if (hasTies) {
-                                       double denom = sqrt((denom1 * denom2));
-                                       p = di / denom;
-                               }else {
-                                       int n = lookupFloat.size();
-                                       p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1)));
-                               }
+                               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)));
+                               
+                if (isnan(p) || isinf(p)) { p = 0.0; }
+                

                                out  << '\t' << p;
                                out  << '\t' << p;

+                               
+                               pValues[j] = p;
+                
+                double sig = linear.calcSpearmanSig(n, sf[j], sg, di);            
+                out  << '\t' << sig;
+                

                        }
 
                        }
 

-                       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;


@@ -529,6 +582,8 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
        try {
                
 int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
        try {
                

+        LinearAlgebra linear;
+        

                //format data
                vector< vector<spearmanRank> > scores; scores.resize(numaxes);
                for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
                //format data
                vector< vector<spearmanRank> > scores; scores.resize(numaxes);
                for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {


@@ -548,10 +603,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();


@@ -572,7 +627,7 @@ 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;  }
+                       if (metadatafile == "") {  out << m->currentSharedBinLabels[i]; }

                        else {  out << metadataLabels[i];               }
                        
                        //find the ranks of this otu - Y
                        else {  out << metadataLabels[i];               }
                        
                        //find the ranks of this otu - Y


@@ -588,10 +643,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();


@@ -607,38 +662,49 @@ 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++) {
                        
                        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> 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;
-                               
+                               double p = (numCoor - numDisCoor) / (float) count;
+                 if (isnan(p) || isinf(p)) { p = 0.0; }
+                

                                out << '\t' << p;
                                pValues[j] = p;
                                out << '\t' << p;
                                pValues[j] = p;

-
+                
+                double sig = linear.calcKendallSig(scores[j].size(), p);
+                
+                out << '\t' << sig;

                        }
                        
                        double sum = 0;
                        }
                        
                        double sum = 0;


@@ -743,6 +809,8 @@ int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thisloo
                }
                
                //for each bin
                }
                
                //for each bin

+               vector<string> newBinLabels;
+               string snumBins = toString(thislookup[0]->getNumBins());

                for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
                        if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
                        
                for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
                        if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
                        


@@ -757,12 +825,25 @@ int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thisloo
                                for (int j = 0; j < thislookup.size(); j++) {
                                        newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
                                }
                                for (int j = 0; j < thislookup.size(); j++) {
                                        newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
                                }

+                               
+                               //if there is a bin label use it otherwise make one
+                               string binLabel = "Otu";
+                               string sbinNumber = toString(i+1);
+                               if (sbinNumber.length() < snumBins.length()) { 
+                                       int diff = snumBins.length() - sbinNumber.length();
+                                       for (int h = 0; h < diff; h++) { binLabel += "0"; }
+                               }
+                               binLabel += sbinNumber; 
+                               if (i < m->currentSharedBinLabels.size()) {  binLabel = m->currentSharedBinLabels[i]; }
+                               
+                               newBinLabels.push_back(binLabel);

                        }
                }
                
                for (int j = 0; j < thislookup.size(); j++) {  delete thislookup[j];  }
                
                thislookup = newLookup;
                        }
                }
                
                for (int j = 0; j < thislookup.size(); j++) {  delete thislookup[j];  }
                
                thislookup = newLookup;

+               m->currentSharedBinLabels = newBinLabels;

                
                return 0;
                
                
                return 0;
                


@@ -839,22 +920,17 @@ int CorrAxesCommand::getMetadata(){
                vector<string> groupNames;
                
                ifstream in;
                vector<string> groupNames;
                
                ifstream in;

-               m->openInputFile(axesfile, in);
+               m->openInputFile(metadatafile, in);

                
                string headerLine = m->getline(in); m->gobble(in);
                
                string headerLine = m->getline(in); m->gobble(in);

-               istringstream iss (headerLine,istringstream::in);
-               
-               //read the first label, because it refers to the groups
-               string columnLabel;
-               iss >> columnLabel; m->gobble(iss); 
+               vector<string> pieces = m->splitWhiteSpace(headerLine);

                
                //save names of columns you are reading
                
                //save names of columns you are reading

-               while (!iss.eof()) {
-                       iss >> columnLabel; m->gobble(iss);
-                       metadataLabels.push_back(columnLabel);
+               for (int i = 1; i < pieces.size(); i++) {
+                       metadataLabels.push_back(pieces[i]);

                }
                int count = metadataLabels.size();
                }
                int count = metadataLabels.size();

-               
+                       

                //read rest of file
                while (!in.eof()) {
                        
                //read rest of file
                while (!in.eof()) {
                        


@@ -863,7 +939,7 @@ int CorrAxesCommand::getMetadata(){
                        string group = "";
                        in >> group; m->gobble(in);
                        groupNames.push_back(group);
                        string group = "";
                        in >> group; m->gobble(in);
                        groupNames.push_back(group);

-                       
+                               

                        SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
                        tempLookup->setGroup(group);
                        tempLookup->setLabel("1");
                        SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
                        tempLookup->setGroup(group);
                        tempLookup->setLabel("1");


@@ -871,7 +947,6 @@ int CorrAxesCommand::getMetadata(){
                        for (int i = 0; i < count; i++) {
                                float temp = 0.0;
                                in >> temp; 
                        for (int i = 0; i < count; i++) {
                                float temp = 0.0;
                                in >> temp; 

-                               

                                tempLookup->push_back(temp, group);
                        }
                        
                                tempLookup->push_back(temp, group);
                        }
                        


@@ -884,12 +959,13 @@ int CorrAxesCommand::getMetadata(){
                //remove any groups the user does not want, and set globaldata->groups with only valid groups
                SharedUtil* util;
                util = new SharedUtil();
                //remove any groups the user does not want, and set globaldata->groups with only valid groups
                SharedUtil* util;
                util = new SharedUtil();

-               
-               util->setGroups(globaldata->Groups, groupNames);
+               Groups = m->getGroups();
+               util->setGroups(Groups, groupNames);
+               m->setGroups(Groups);

                
                for (int i = 0; i < lookupFloat.size(); i++) {
                        //if this sharedrabund is not from a group the user wants then delete it.
                
                for (int i = 0; i < lookupFloat.size(); i++) {
                        //if this sharedrabund is not from a group the user wants then delete it.

-                       if (util->isValidGroup(lookupFloat[i]->getGroup(), globaldata->Groups) == false) { 
+                       if (util->isValidGroup(lookupFloat[i]->getGroup(), m->getGroups()) == false) { 

                                delete lookupFloat[i]; lookupFloat[i] = NULL;
                                lookupFloat.erase(lookupFloat.begin()+i); 
                                i--; 
                                delete lookupFloat[i]; lookupFloat[i] = NULL;
                                lookupFloat.erase(lookupFloat.begin()+i); 
                                i--;