]> git.donarmstrong.com Git - mothur.git/blobdiff - corraxescommand.cpp
working on current change
[mothur.git] / corraxescommand.cpp
index 429ef141b54b621bb9bff43904e78038c31df438..cd188b8a921a33aaac255852b5ab22d6008f605e 100644 (file)
@@ -8,41 +8,57 @@
  */
 
 #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(){  
+vector<string> CorrAxesCommand::setParameters(){       
        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",false,true); parameters.push_back(paxes);
+               CommandParameter pshared("shared", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pshared);
+               CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); 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) {
-               m->errorOut(e, "CorrAxesCommand", "getValidParameters");
+               m->errorOut(e, "CorrAxesCommand", "setParameters");
                exit(1);
        }
 }
 //**********************************************************************************************************************
-vector<string> CorrAxesCommand::getRequiredParameters(){       
+string CorrAxesCommand::getHelpString(){       
        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) {
-               m->errorOut(e, "CorrAxesCommand", "getRequiredParameters");
+               m->errorOut(e, "CorrAxesCommand", "getHelpString");
                exit(1);
        }
 }
 //**********************************************************************************************************************
 CorrAxesCommand::CorrAxesCommand(){    
        try {
-               abort = true;
-               //initialize outputTypes
+               abort = true; calledHelp = true; 
+               setParameters();
                vector<string> tempOutNames;
                outputTypes["corr.axes"] = tempOutNames;
        }
@@ -51,31 +67,17 @@ CorrAxesCommand::CorrAxesCommand(){
                exit(1);
        }
 }
-
-//**********************************************************************************************************************
-vector<string> CorrAxesCommand::getRequiredFiles(){    
-       try {
-               vector<string> myArray;
-               return myArray;
-       }
-       catch(exception& e) {
-               m->errorOut(e, "CorrAxesCommand", "getRequiredFiles");
-               exit(1);
-       }
-}
 //**********************************************************************************************************************
 CorrAxesCommand::CorrAxesCommand(string option)  {
        try {
-               abort = false;
-               globaldata = GlobalData::getInstance();
+               abort = false; calledHelp = false;   
                
                //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 {
-                       //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();
@@ -138,16 +140,17 @@ CorrAxesCommand::CorrAxesCommand(string option)  {
                        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 = ""; }
-                       else { inputFileName = relabundfile; }
+                       else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
                        
                        metadatafile = validParameter.validFile(parameters, "metadata", true);
                        if (metadatafile == "not open") { abort = true; }
                        else if (metadatafile == "not found") { metadatafile = ""; }
+                       else { inputFileName = metadatafile; }
                        
                        groups = validParameter.validFile(parameters, "groups", false);                 
                        if (groups == "not found") { groups = "";  pickedGroups = false;  }
@@ -155,17 +158,33 @@ CorrAxesCommand::CorrAxesCommand(string option)  {
                                pickedGroups = true;
                                m->splitAtDash(groups, Groups); 
                        }                       
-                       globaldata->Groups = Groups;
+                       m->Groups = 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=""; }  
                        
-                       if ((relabundfile == "") && (sharedfile == "")) { m->mothurOut("You must provide either a shared or relabund file."); m->mothurOutEndLine(); abort = true;  }
-                       
-                       if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may not use both a shared and relabund 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;  }
+                       }else {
+                               if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared, relabund or metadata file."); m->mothurOutEndLine(); abort = true;  }
+                       }
                        string temp;
                        temp = validParameter.validFile(parameters, "numaxes", false);          if (temp == "not found"){       temp = "3";                             }
                        convert(temp, numaxes); 
@@ -183,56 +202,38 @@ CorrAxesCommand::CorrAxesCommand(string option)  {
 }
 //**********************************************************************************************************************
 
-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 parameters are shared, relabund, axes, metadata, groups, method, numaxes and label.  The shared or relabund and axes parameters are required.  If shared is given the relative abundance is calculated.\n");
-               m->mothurOut("The metadata parameter.....\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 {
                
-               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  //
                /************************************************************************************/
                if (sharedfile != "") {  
-                       getShared(); 
-                       if (m->control_pressed) {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } return 0; }
-                       if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
+                       InputData* input = new InputData(sharedfile, "sharedfile");
+                       getSharedFloat(input); 
+                       delete input;
+                       
+                       if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
+                       if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
                        
-                       //fills lookupFloat with relative abundance values from lookup
-                       convertToRelabund();
+               }else if (relabundfile != "") { 
+                       InputData* input = new InputData(relabundfile, "relabund");
+                       getSharedFloat(input); 
+                       delete input;
                        
-                       for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
-               }else { 
-                       getSharedFloat(); 
                        if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
                        if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
                        
+               }else if (metadatafile != "") { 
+                       getMetadata();  //reads metadata file and store in lookupFloat, saves column headings in metadataLabels for later
+                       if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
+                       if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading metadata file."); m->mothurOutEndLine(); return 0; }
+                       
                        if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
-               }
+                       
+               }else { m->mothurOut("[ERROR]: no file given."); m->mothurOutEndLine(); return 0; }
                
                if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
                
@@ -271,9 +272,11 @@ int CorrAxesCommand::execute(){
                out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
                
                //output headings
-               out << "OTU\t";
-               for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; }
-               out << endl;
+               if (metadatafile == "") {  out << "OTU";        }
+               else {  out << "Feature";                                               }
+
+               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); }
@@ -315,8 +318,9 @@ int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& ou
           //for each otu
           for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
                   
-                  out << i+1 << '\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++) {
@@ -324,6 +328,8 @@ int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& ou
                   }
                   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++) {
                           
@@ -343,11 +349,15 @@ int CorrAxesCommand::calcPearson(map<string, vector<float> >& axes, ofstream& ou
                           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;
@@ -362,12 +372,32 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
        try {
                
                //format data
+               vector< map<float, int> > tableX; tableX.resize(numaxes);
+               map<float, int>::iterator itTable;
                vector< vector<spearmanRank> > scores; scores.resize(numaxes);
                for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
                        vector<float> temp = it->second;
                        for (int i = 0; i < temp.size(); i++) {
                                spearmanRank member(it->first, temp[i]);
                                scores[i].push_back(member);  
+                               
+                               //count number of repeats
+                               itTable = tableX[i].find(temp[i]);
+                               if (itTable == tableX[i].end()) { 
+                                       tableX[i][temp[i]] = 1;
+                               }else {
+                                       tableX[i][temp[i]]++;
+                               }
+                       }
+               }
+               
+               //calc LX
+               //for each axis
+               vector<double> Lx; Lx.resize(numaxes, 0.0);
+               for (int i = 0; i < numaxes; i++) {
+                       for (itTable = tableX[i].begin(); itTable != tableX[i].end(); itTable++) {
+                               double tx = (double) itTable->second;
+                               Lx[i] += ((pow(tx, 3.0) - tx) / 12.0);
                        }
                }
                
@@ -381,11 +411,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++) {
-                               rankTotal += j;
+                               rankTotal += (j+1);
                                ties.push_back(scores[i][j]);
                                
-                               if (j != scores.size()-1) { // you are not the last so you can look ahead
+                               if (j != (scores[i].size()-1)) { // you are not the last so you can look ahead
                                        if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
+
                                                for (int k = 0; k < ties.size(); k++) {
                                                        float thisrank = rankTotal / (float) ties.size();
                                                        rankAxes[ties[k].name].push_back(thisrank);
@@ -394,26 +425,44 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                                                rankTotal = 0;
                                        }
                                }else { // you are the last one
+                                       
                                        for (int k = 0; k < ties.size(); k++) {
                                                float thisrank = rankTotal / (float) ties.size();
                                                rankAxes[ties[k].name].push_back(thisrank);
+                                               
                                        }
                                }
                        }
                }
                
                                
-               
                //for each otu
                for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
                        
-                       out << i+1 << '\t';
+                       if (metadatafile == "") {  out << i+1;  }
+                       else {  out << metadataLabels[i];               }
                        
                        //find the ranks of this otu - Y
                        vector<spearmanRank> otuScores;
+                       map<float, int> tableY;
                        for (int j = 0; j < lookupFloat.size(); j++) {
                                spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
                                otuScores.push_back(member);
+                               
+                               itTable = tableY.find(member.score);
+                               if (itTable == tableY.end()) { 
+                                       tableY[member.score] = 1;
+                               }else {
+                                       tableY[member.score]++;
+                               }
+                               
+                       }
+                       
+                       //calc Ly
+                       double Ly = 0.0;
+                       for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
+                               double ty = (double) itTable->second;
+                               Ly += ((pow(ty, 3.0) - ty) / 12.0);
                        }
                        
                        sort(otuScores.begin(), otuScores.end(), compareSpearman);
@@ -422,11 +471,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++) {
-                               rankTotal += j;
+                               rankTotal += (j+1);
                                ties.push_back(otuScores[j]);
                                
-                               if (j != scores.size()-1) { // you are not the last so you can look ahead
+                               if (j != (otuScores.size()-1)) { // you are not the last so you can look ahead
                                        if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
+                                               
                                                for (int k = 0; k < ties.size(); k++) {
                                                        float thisrank = rankTotal / (float) ties.size();
                                                        rankOtus[ties[k].name] = thisrank;
@@ -435,13 +485,15 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                                                rankTotal = 0;
                                        }
                                }else { // you are the last one
+                                       
                                        for (int k = 0; k < ties.size(); k++) {
                                                float thisrank = rankTotal / (float) ties.size();
                                                rankOtus[ties[k].name] = thisrank;
                                        }
                                }
                        }
-                       
+                       vector<double> pValues(numaxes);        
+
                        //calc spearman ranks for each axis for this otu
                        for (int j = 0; j < numaxes; j++) {
                                
@@ -454,14 +506,25 @@ int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& o
                                        di += ((xi - yi) * (xi - yi));
                                }
                                
-                               int n = lookupFloat.size();
-                               double p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1)));
+                               double p = 0.0;
+                               
+                               double n = (double) lookupFloat.size();
+                               
+                               double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx[j];
+                               double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
+                               
+                               p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
+                               
+                               out  << '\t' << p;
                                
-                               out << p << '\t';
+                               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;
@@ -494,10 +557,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++) {
-                               rankTotal += j;
+                               rankTotal += (j+1);
                                ties.push_back(&(scores[i][j]));
                                
-                               if (j != scores.size()-1) { // you are not the last so you can look ahead
+                               if (j != scores[i].size()-1) { // you are not the last so you can look ahead
                                        if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue
                                                for (int k = 0; k < ties.size(); k++) {
                                                        float thisrank = rankTotal / (float) ties.size();
@@ -518,7 +581,8 @@ int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& ou
                //for each otu
                for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
                
-                       out << i+1 << '\t';
+                       if (metadatafile == "") {  out << i+1;  }
+                       else {  out << metadataLabels[i];               }
                        
                        //find the ranks of this otu - Y
                        vector<spearmanRank> otuScores;
@@ -533,10 +597,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++) {
-                               rankTotal += j;
+                               rankTotal += (j+1);
                                ties.push_back(otuScores[j]);
                                
-                               if (j != scores.size()-1) { // you are not the last so you can look ahead
+                               if (j != otuScores.size()-1) { // you are not the last so you can look ahead
                                        if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue
                                                for (int k = 0; k < ties.size(); k++) {
                                                        float thisrank = rankTotal / (float) ties.size();
@@ -553,130 +617,67 @@ 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;
-                               
-                               out << p << '\t';
-                       }
-                       
-                       out << endl;
-               }
-               
-               return 0;
-       }
-       catch(exception& e) {
-               m->errorOut(e, "CorrAxesCommand", "calcKendall");
-               exit(1);
-       }
-}
-//**********************************************************************************************************************
-int CorrAxesCommand::getShared(){
-       try {
-               InputData* input = new InputData(sharedfile, "sharedfile");
-               lookup = input->getSharedRAbundVectors();
-               string lastLabel = lookup[0]->getLabel();
-               
-               if (label == "") { label = lastLabel; delete input; return 0; }
-               
-               //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
-               set<string> labels; labels.insert(label);
-               set<string> processedLabels;
-               set<string> userLabels = labels;
-               
-               //as long as you are not at the end of the file or done wih the lines you want
-               while((lookup[0] != NULL) && (userLabels.size() != 0)) {
-                       if (m->control_pressed) {  delete input; return 0;  }
-                       
-                       if(labels.count(lookup[0]->getLabel()) == 1){
-                               processedLabels.insert(lookup[0]->getLabel());
-                               userLabels.erase(lookup[0]->getLabel());
-                               break;
-                       }
-                       
-                       if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
-                               string saveLabel = lookup[0]->getLabel();
-                               
-                               for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
-                               lookup = input->getSharedRAbundVectors(lastLabel);
-                               
-                               processedLabels.insert(lookup[0]->getLabel());
-                               userLabels.erase(lookup[0]->getLabel());
-                               
-                               //restore real lastlabel to save below
-                               lookup[0]->setLabel(saveLabel);
-                               break;
+                               double p = (numCoor - numDisCoor) / (float) count;
+
+                               out << '\t' << p;
+                               pValues[j] = p;
+
                        }
                        
-                       lastLabel = lookup[0]->getLabel();                      
-                       
-                       //get next line to process
-                       //prevent memory leak
-                       for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
-                       lookup = input->getSharedRAbundVectors();
-               }
-               
-               
-               if (m->control_pressed) { delete input; return 0;  }
-               
-               //output error messages about any remaining user labels
-               set<string>::iterator it;
-               bool needToRun = false;
-               for (it = userLabels.begin(); it != userLabels.end(); it++) {  
-                       m->mothurOut("Your file does not include the label " + *it); 
-                       if (processedLabels.count(lastLabel) != 1) {
-                               m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
-                               needToRun = true;
-                       }else {
-                               m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
+                       double sum = 0;
+                       for(int k=0;k<numaxes;k++){
+                               sum += pValues[k] * pValues[k];
                        }
+                       out << '\t' << sqrt(sum) << endl;
                }
                
-               //run last label if you need to
-               if (needToRun == true)  {
-                       for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } } 
-                       lookup = input->getSharedRAbundVectors(lastLabel);
-               }       
-               
-               delete input;
                return 0;
        }
        catch(exception& e) {
-               m->errorOut(e, "CorrAxesCommand", "getShared"); 
+               m->errorOut(e, "CorrAxesCommand", "calcKendall");
                exit(1);
        }
 }
 //**********************************************************************************************************************
-int CorrAxesCommand::getSharedFloat(){
+int CorrAxesCommand::getSharedFloat(InputData* input){
        try {
-               InputData* input = new InputData(relabundfile, "relabund");
                lookupFloat = input->getSharedRAbundFloatVectors();
                string lastLabel = lookupFloat[0]->getLabel();
                
-               if (label == "") { label = lastLabel; delete input; return 0; }
+               if (label == "") { label = lastLabel;  return 0; }
                
                //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
                set<string> labels; labels.insert(label);
@@ -686,7 +687,7 @@ int CorrAxesCommand::getSharedFloat(){
                //as long as you are not at the end of the file or done wih the lines you want
                while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
                        
-                       if (m->control_pressed) {  delete input; return 0;  }
+                       if (m->control_pressed) {  return 0;  }
                        
                        if(labels.count(lookupFloat[0]->getLabel()) == 1){
                                processedLabels.insert(lookupFloat[0]->getLabel());
@@ -717,7 +718,7 @@ int CorrAxesCommand::getSharedFloat(){
                }
                
                
-               if (m->control_pressed) { delete input; return 0;  }
+               if (m->control_pressed) { return 0;  }
                
                //output error messages about any remaining user labels
                set<string>::iterator it;
@@ -738,7 +739,6 @@ int CorrAxesCommand::getSharedFloat(){
                        lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
                }       
                
-               delete input;
                return 0;
        }
        catch(exception& e) {
@@ -746,43 +746,6 @@ int CorrAxesCommand::getSharedFloat(){
                exit(1);
        }
 }
-/*****************************************************************/
-int CorrAxesCommand::convertToRelabund(){
-       try {
-               
-               vector<SharedRAbundFloatVector*> newLookup;
-               for (int i = 0; i < lookup.size(); i++) {
-                       SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
-                       temp->setLabel(lookup[i]->getLabel());
-                       temp->setGroup(lookup[i]->getGroup());
-                       newLookup.push_back(temp);
-               }
-               
-               for (int i = 0; i < lookup.size(); i++) {
-                       
-                       for (int j = 0; j < lookup[i]->getNumBins(); j++) {
-                               
-                               if (m->control_pressed) { return 0; }
-                               
-                               int abund = lookup[i]->getAbundance(j);
-                               
-                               float relabund = abund / (float) lookup[i]->getNumSeqs();
-                               
-                               newLookup[i]->push_back(relabund, lookup[i]->getGroup());
-                       }
-               }
-               
-               if (pickedGroups) { eliminateZeroOTUS(newLookup); }
-               
-               lookupFloat = newLookup;
-               
-               return 0;
-       }
-       catch(exception& e) {
-               m->errorOut(e, "CorrAxesCommand", "convertToRelabund"); 
-               exit(1);
-       }
-}
 //**********************************************************************************************************************
 int CorrAxesCommand::eliminateZeroOTUS(vector<SharedRAbundFloatVector*>& thislookup) {
        try {
@@ -882,7 +845,79 @@ map<string, vector<float> > CorrAxesCommand::readAxes(){
                return axes;
        }
        catch(exception& e) {
-               m->errorOut(e, "CorrAxesCommand", "convertToRelabund"); 
+               m->errorOut(e, "CorrAxesCommand", "readAxes");  
+               exit(1);
+       }
+}
+/*****************************************************************/
+int CorrAxesCommand::getMetadata(){
+       try {
+               vector<string> groupNames;
+               
+               ifstream in;
+               m->openInputFile(metadatafile, 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); 
+
+               //save names of columns you are reading
+               while (!iss.eof()) {
+                       iss >> columnLabel; m->gobble(iss);
+                       metadataLabels.push_back(columnLabel);
+               }
+               int count = metadataLabels.size();
+               
+               //read rest of file
+               while (!in.eof()) {
+                       
+                       if (m->control_pressed) { in.close(); return 0; }
+                       
+                       string group = "";
+                       in >> group; m->gobble(in);
+                       groupNames.push_back(group);
+                       
+                       SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector();
+                       tempLookup->setGroup(group);
+                       tempLookup->setLabel("1");
+                       
+                       for (int i = 0; i < count; i++) {
+                               float temp = 0.0;
+                               in >> temp; 
+                               
+                               tempLookup->push_back(temp, group);
+                       }
+                       
+                       lookupFloat.push_back(tempLookup);
+                       
+                       m->gobble(in);
+               }
+               in.close();
+               
+               //remove any groups the user does not want, and set globaldata->groups with only valid groups
+               SharedUtil* util;
+               util = new SharedUtil();
+               
+               util->setGroups(m->Groups, groupNames);
+               
+               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(), m->Groups) == false) { 
+                               delete lookupFloat[i]; lookupFloat[i] = NULL;
+                               lookupFloat.erase(lookupFloat.begin()+i); 
+                               i--; 
+                       }
+               }
+               
+               delete util;
+               
+               return 0;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "CorrAxesCommand", "getMetadata");       
                exit(1);
        }
 }