#include "corraxescommand.h"
#include "sharedutilities.h"
-//********************************************************************************************************************
-//sorts highest to lowest
-inline bool compareSpearman(spearmanRank left, spearmanRank right){
- return (left.score > right.score);
-}
-//********************************************************************************************************************
-//sorts lowest to highest
-inline bool compareSpearmanReverse(spearmanRank left, spearmanRank right){
- return (left.score < right.score);
-}
//**********************************************************************************************************************
-vector<string> CorrAxesCommand::getValidParameters(){
+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;
}
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 {
- //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();
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 == "") && (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; }
}
//**********************************************************************************************************************
-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 {
- 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 //
int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
try {
- /*axes.clear();
- axes["1a"].push_back(1);
- axes["1b"].push_back(1);
- axes["1c"].push_back(1);
- axes["2a"].push_back(2);
- axes["2b"].push_back(2);
- axes["2c"].push_back(2);
- axes["3a"].push_back(3);
- axes["3b"].push_back(3);
- axes["3c"].push_back(3);
- axes["4a"].push_back(4);
- axes["4b"].push_back(4);
- axes["4c"].push_back(4);
- axes["5a"].push_back(5);
- axes["5b"].push_back(5);
- axes["5c"].push_back(5);*/
-
//format data
vector< map<float, int> > tableX; tableX.resize(numaxes);
map<float, int>::iterator itTable;
int CorrAxesCommand::calcKendall(map<string, vector<float> >& axes, ofstream& out) {
try {
- numaxes = 1;
- axes.clear();
- axes["1a"].push_back(1);
- axes["1b"].push_back(1);
- axes["1c"].push_back(1);
- axes["2a"].push_back(2);
- axes["2b"].push_back(2);
- axes["2c"].push_back(2);
- axes["3a"].push_back(3);
- axes["3b"].push_back(3);
- axes["3c"].push_back(3);
- axes["4a"].push_back(4);
- axes["4b"].push_back(4);
- axes["4c"].push_back(4);
- axes["5a"].push_back(5);
- axes["5b"].push_back(5);
- axes["5c"].push_back(5);
-
-
//format data
vector< vector<spearmanRank> > scores; scores.resize(numaxes);
for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
}
}
}
- cout << "axes ranks = ";
- for(int i = 0; i < scores[0].size(); i++) { cout << scores[0][i].score << '\t'; }
- cout << endl;
- lookupFloat.clear();
- lookupFloat.resize(15);
- for (int i = 0; i < lookupFloat.size(); i++) {
- lookupFloat[i] = new SharedRAbundFloatVector();
- }
- lookupFloat[0]->push_back(0.2288227, "1a"); lookupFloat[0]->setGroup("1a");
- lookupFloat[1]->push_back(0.7394062, "1b"); lookupFloat[1]->setGroup("1b");
- lookupFloat[2]->push_back(0.4521187, "1c"); lookupFloat[2]->setGroup("1c");
- lookupFloat[3]->push_back(0.1598630, "2a"); lookupFloat[3]->setGroup("2a");
- lookupFloat[4]->push_back(0.09588156, "2b"); lookupFloat[4]->setGroup("2b");
-
- lookupFloat[5]->push_back(0.933174, "2c"); lookupFloat[5]->setGroup("2c");
- lookupFloat[6]->push_back(0.3958304, "3a"); lookupFloat[6]->setGroup("3a");;
- lookupFloat[7]->push_back(0.2364419, "3b"); lookupFloat[7]->setGroup("3b");
- lookupFloat[8]->push_back(0.1697712, "3c"); lookupFloat[8]->setGroup("3c");
- lookupFloat[9]->push_back(0.4077173, "4a"); lookupFloat[9]->setGroup("4a");
-
- lookupFloat[10]->push_back(0.6116547, "4b"); lookupFloat[10]->setGroup("4b");
- lookupFloat[11]->push_back(0.9374322, "4c"); lookupFloat[11]->setGroup("4c");
- lookupFloat[12]->push_back(0.852184, "5a"); lookupFloat[12]->setGroup("5a");
- lookupFloat[13]->push_back(0.845094, "5b"); lookupFloat[13]->setGroup("5b");
- lookupFloat[14]->push_back(0.5795778, "5c"); lookupFloat[14]->setGroup("5c");
-
//for each otu
for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
int numCoor = 0;
int numDisCoor = 0;
- //assemble otus ranks in same order as axis ranks, if ties, put in the best order possible
- //NOTE: after this ordering the scores[j] indexes may not match the otus indexes within tied sections
- //since we do not use the axes ranks except to order the otu ranks, I did not take the time to reorder the axes ranks
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]);
- otusTemp.push_back(member);
-
- if (l != (scores[j].size()-1)) { // you are not the last so you can look ahead
- if (scores[j][l].score != scores[j][l+1].score) { // you are done with ties, order them and continue
- //order otus within tied section in the best way possible to make coor pairs
- sort(otusTemp.begin(), otusTemp.end(), compareSpearmanReverse);
-
- //save order in otus
- for (int h = 0; h < otusTemp.size(); h++) { ; otus.push_back(otusTemp[h]); }
-
- otusTemp.clear();
- }
- }else { // you are the last one
- //order otus within tied section in the best way possible to make coor pairs
- sort(otusTemp.begin(), otusTemp.end(), compareSpearmanReverse);
-
- //save order in otus
- for (int h = 0; h < otusTemp.size(); h++) { otus.push_back(otusTemp[h]); }
- }
+ otus.push_back(member);
}
- cout << "otu ranks = ";
- for (int h = 0; h < otus.size(); h++ ) { cout << otus[h].score << '\t'; }
- cout << endl;
-
int count = 0;
for (int l = 0; l < scores[j].size(); l++) {
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++;
numDisCoor += numWithLowerRank;
}
- //sample size
- //int n = lookupFloat.size();
- //comparing to yourself
- count -= lookupFloat.size();
-
- //double p = ( (4 * P) / (float) ((n * (n - 1)) - numTies) ) - 1.0;
double p = (numCoor - numDisCoor) / (float) count;
- cout << "numCoor = " << numCoor << " numDisCoor = " << numDisCoor << " p = " << p << " count = " << count << endl;
+
out << '\t' << p;
pValues[j] = p;
vector<string> groupNames;
ifstream in;
- m->openInputFile(axesfile, 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);
SharedUtil* util;
util = new SharedUtil();
- util->setGroups(globaldata->Groups, groupNames);
+ 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(), globaldata->Groups) == false) {
+ if (util->isValidGroup(lookupFloat[i]->getGroup(), m->Groups) == false) {
delete lookupFloat[i]; lookupFloat[i] = NULL;
lookupFloat.erase(lookupFloat.begin()+i);
i--;
projects / mothur.git / blobdiff