X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=corraxescommand.cpp;h=10669e580998797222fd785b6d80608fd87ba036;hb=fb981b4699571698cbbb09e9ee9c986ca49c3d58;hp=b6c1a705a39a891072f147997f942708fcc52272;hpb=5d6d303e481489e226fdf8d6c5385b99b50718bc;p=mothur.git diff --git a/corraxescommand.cpp b/corraxescommand.cpp index b6c1a70..10669e5 100644 --- a/corraxescommand.cpp +++ b/corraxescommand.cpp @@ -10,40 +10,55 @@ #include "corraxescommand.h" #include "sharedutilities.h" -//******************************************************************************************************************** -//sorts highes to lowest -inline bool compareSpearman(spearmanRank left, spearmanRank right){ - return (left.score > right.score); -} //********************************************************************************************************************** -vector CorrAxesCommand::getValidParameters(){ +vector CorrAxesCommand::setParameters(){ try { - string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"}; - vector 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 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 CorrAxesCommand::getRequiredParameters(){ +string CorrAxesCommand::getHelpString(){ try { - string Array[] = {"axes"}; - vector 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 tempOutNames; outputTypes["corr.axes"] = tempOutNames; } @@ -52,31 +67,17 @@ CorrAxesCommand::CorrAxesCommand(){ exit(1); } } - -//********************************************************************************************************************** -vector CorrAxesCommand::getRequiredFiles(){ - try { - vector 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 myArray (Array, Array+(sizeof(Array)/sizeof(string))); + vector myArray = setParameters(); OptionParser parser(option); map parameters = parser.getParameters(); @@ -139,12 +140,12 @@ 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; } @@ -157,14 +158,27 @@ CorrAxesCommand::CorrAxesCommand(string option) { 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=""; } - 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; } @@ -173,7 +187,7 @@ CorrAxesCommand::CorrAxesCommand(string option) { } 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"; } @@ -188,35 +202,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 { - 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 // @@ -286,7 +275,7 @@ int CorrAxesCommand::execute(){ 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); } @@ -360,8 +349,19 @@ int CorrAxesCommand::calcPearson(map >& axes, ofstream& ou double denom = (sqrt(denomTerm1) * sqrt(denomTerm2)); r = numerator / denom; + + if (isnan(r) || isinf(r)) { r = 0.0; } + rValues[k] = r; out << '\t' << r; + + //signifigance calc - http://faculty.vassar.edu/lowry/ch4apx.html + double temp = (1- (r*r)) / (double) (lookupFloat.size()-2); + temp = sqrt(temp); + double sig = r / temp; + if (isnan(sig) || isinf(sig)) { sig = 0.0; } + + out << '\t' << sig; } double sum = 0; @@ -382,23 +382,6 @@ int CorrAxesCommand::calcPearson(map >& axes, ofstream& ou int CorrAxesCommand::calcSpearman(map >& 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 > tableX; tableX.resize(numaxes); map::iterator itTable; @@ -463,30 +446,7 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o } } - /*cout << endl; - lookupFloat.clear(); - lookupFloat.resize(15); - for (int i = 0; i < lookupFloat.size(); i++) { - lookupFloat[i] = new SharedRAbundFloatVector(); - } - lookupFloat[0]->push_back(0.2288227, "1a"); lookupFloat[0]->setGroup("1a"); - lookupFloat[1]->push_back(0.7394062, "1b"); lookupFloat[1]->setGroup("1b"); - lookupFloat[2]->push_back(0.4521187, "1c"); lookupFloat[2]->setGroup("1c"); - lookupFloat[3]->push_back(0.1598630, "2a"); lookupFloat[3]->setGroup("2a"); - lookupFloat[4]->push_back(0.09588156, "2b"); lookupFloat[4]->setGroup("2b"); - - lookupFloat[5]->push_back(0.933174, "2c"); lookupFloat[5]->setGroup("2c"); - lookupFloat[6]->push_back(0.3958304, "3a"); lookupFloat[6]->setGroup("3a");; - lookupFloat[7]->push_back(0.2364419, "3b"); lookupFloat[7]->setGroup("3b"); - lookupFloat[8]->push_back(0.1697712, "3c"); lookupFloat[8]->setGroup("3c"); - lookupFloat[9]->push_back(0.4077173, "4a"); lookupFloat[9]->setGroup("4a"); - - lookupFloat[10]->push_back(0.6116547, "4b"); lookupFloat[10]->setGroup("4b"); - lookupFloat[11]->push_back(0.9374322, "4c"); lookupFloat[11]->setGroup("4c"); - lookupFloat[12]->push_back(0.852184, "5a"); lookupFloat[12]->setGroup("5a"); - lookupFloat[13]->push_back(0.845094, "5b"); lookupFloat[13]->setGroup("5b"); - lookupFloat[14]->push_back(0.5795778, "5c"); lookupFloat[14]->setGroup("5c");*/ - + //for each otu for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) { @@ -566,9 +526,20 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2))); + if (isnan(p) || isinf(p)) { p = 0.0; } + out << '\t' << p; pValues[j] = p; + + //signifigance calc - http://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient + double temp = (lookupFloat.size()-2) / (double) (1- (p*p)); + temp = sqrt(temp); + double sig = p*temp; + if (isnan(sig) || isinf(sig)) { sig = 0.0; } + + out << '\t' << sig; + } double sum = 0; @@ -608,10 +579,10 @@ int CorrAxesCommand::calcKendall(map >& axes, ofstream& ou vector 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(); @@ -648,10 +619,10 @@ int CorrAxesCommand::calcKendall(map >& axes, ofstream& ou vector 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(); @@ -667,38 +638,56 @@ int CorrAxesCommand::calcKendall(map >& axes, ofstream& ou } } } + + vector 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 otus; + vector otusTemp; for (int l = 0; l < scores[j].size(); l++) { spearmanRank member(scores[j][l].name, rankOtus[scores[j][l].name]); otus.push_back(member); } + int count = 0; for (int l = 0; l < scores[j].size(); l++) { int numWithHigherRank = 0; + int numWithLowerRank = 0; float thisrank = otus[l].score; - for (int u = l; u < scores[j].size(); u++) { + for (int u = l+1; u < scores[j].size(); u++) { if (otus[u].score > thisrank) { numWithHigherRank++; } + else if (otus[u].score < thisrank) { numWithLowerRank++; } + count++; } - P += numWithHigherRank; + numCoor += numWithHigherRank; + numDisCoor += numWithLowerRank; } - int n = lookupFloat.size(); - - double p = ( (4 * P) / (float) (n * (n - 1)) ) - 1.0; - + double p = (numCoor - numDisCoor) / (float) count; + if (isnan(p) || isinf(p)) { p = 0.0; } + out << '\t' << p; pValues[j] = p; - + + //calc signif - zA - http://en.wikipedia.org/wiki/Kendall_tau_rank_correlation_coefficient#Significance_tests + double numer = 3.0 * (numCoor - numDisCoor); + int n = scores[j].size(); + double denom = n * (n-1) * (2*n + 5) / (double) 2.0; + denom = sqrt(denom); + double sig = numer / denom; + + if (isnan(sig) || isinf(sig)) { sig = 0.0; } + + out << '\t' << sig; } double sum = 0; @@ -803,6 +792,8 @@ int CorrAxesCommand::eliminateZeroOTUS(vector& thisloo } //for each bin + vector 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; } @@ -817,12 +808,25 @@ int CorrAxesCommand::eliminateZeroOTUS(vector& thisloo 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->currentBinLabels.size()) { binLabel = m->currentBinLabels[i]; } + + newBinLabels.push_back(binLabel); } } for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; } thislookup = newLookup; + m->currentBinLabels = newBinLabels; return 0; @@ -899,7 +903,7 @@ int CorrAxesCommand::getMetadata(){ vector groupNames; ifstream in; - m->openInputFile(axesfile, in); + m->openInputFile(metadatafile, in); string headerLine = m->getline(in); m->gobble(in); istringstream iss (headerLine,istringstream::in); @@ -914,7 +918,7 @@ int CorrAxesCommand::getMetadata(){ metadataLabels.push_back(columnLabel); } int count = metadataLabels.size(); - + //read rest of file while (!in.eof()) { @@ -923,7 +927,7 @@ int CorrAxesCommand::getMetadata(){ string group = ""; in >> group; m->gobble(in); groupNames.push_back(group); - + SharedRAbundFloatVector* tempLookup = new SharedRAbundFloatVector(); tempLookup->setGroup(group); tempLookup->setLabel("1"); @@ -931,7 +935,6 @@ int CorrAxesCommand::getMetadata(){ for (int i = 0; i < count; i++) { float temp = 0.0; in >> temp; - tempLookup->push_back(temp, group); } @@ -944,12 +947,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(); - - 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. - 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--;