*/
#include "corraxescommand.h"
+#include "sharedutilities.h"
//********************************************************************************************************************
//sorts highes to lowest
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; }
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 == "")) { 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);
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, 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");
//fills lookupFloat with relative abundance values from lookup
convertToRelabund();
- for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
- }else {
+ for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
+
+ }else if (relabundfile != "") {
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; }
if (pickedGroups) { eliminateZeroOTUS(lookupFloat); }
- }
+
+ }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; }
// calc the r values //
/************************************************************************************/
- string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + "corr.axes";
+ string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + method + ".corr.axes";
outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName);
ofstream out;
m->openOutputFile(outputFileName, out);
out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
//output headings
- out << "OTU\t";
+ if (metadatafile == "") { out << "OTU\t"; }
+ else { out << "Feature\t"; }
+
for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; }
out << endl;
//for each otu
for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
- out << i+1 << '\t';
+ if (metadatafile == "") { out << i+1 << '\t'; }
+ else { out << metadataLabels[i] << '\t'; }
//find the averages this otu - Y
float sumOtu = 0.0;
//for each otu
for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
- out << i+1 << '\t';
+ if (metadatafile == "") { out << i+1 << '\t'; }
+ else { out << metadataLabels[i] << '\t'; }
//find the ranks of this otu - Y
vector<spearmanRank> otuScores;
//sort each axis
for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); }
- //find ranks of xi in each axis
- map<string, vector<float> > rankAxes;
+ //convert scores to ranks of xi in each axis
for (int i = 0; i < numaxes; i++) {
- vector<spearmanRank> ties;
+ vector<spearmanRank*> ties;
int rankTotal = 0;
for (int j = 0; j < scores[i].size(); j++) {
rankTotal += j;
- 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 (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);
+ (*ties[k]).score = thisrank;
}
ties.clear();
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);
+ (*ties[k]).score = thisrank;
}
}
}
}
-
-
//for each otu
for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
-
- out << i+1 << '\t';
+
+ if (metadatafile == "") { out << i+1 << '\t'; }
+ else { out << metadataLabels[i] << '\t'; }
//find the ranks of this otu - Y
vector<spearmanRank> otuScores;
spearmanRank member(lookupFloat[j]->getGroup(), lookupFloat[j]->getAbundance(i));
otuScores.push_back(member);
}
-
+
sort(otuScores.begin(), otuScores.end(), compareSpearman);
map<string, float> rankOtus;
}
}
- //calc kendall coeffient for each axis for this otu
+ //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
+ vector<spearmanRank> otus;
+ 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 numConcordant = 0;
- int numDiscordant = 0;
-
- for (int f = 0; f < j; f++) {
+ for (int l = 0; l < scores[j].size(); l++) {
- for (int k = 0; k < lookupFloat.size(); k++) {
-
- float xi = rankAxes[lookupFloat[k]->getGroup()][j];
- float yi = rankOtus[lookupFloat[k]->getGroup()];
-
- for (int h = 0; h < k; h++) {
-
- float xj = rankAxes[lookupFloat[k]->getGroup()][f];
- float yj = rankOtus[lookupFloat[h]->getGroup()];
-
- if ( ((xi > xj) && (yi < yj)) || ((xi < xj) && (yi > yj)) ){ numDiscordant++; }
- if ( ((xi > xj) && (yi > yj)) || ((xi < xj) && (yi < yj)) ){ numConcordant++; }
- }
+ int numWithHigherRank = 0;
+ float thisrank = otus[l].score;
+
+ for (int u = l; u < scores[j].size(); u++) {
+ if (otus[u].score > thisrank) { numWithHigherRank++; }
}
+
+ P += numWithHigherRank;
}
int n = lookupFloat.size();
- double p = (numConcordant - numDiscordant) / (float) (0.5 * n * (n - 1));
+
+ double p = ( (4 * P) / (float) (n * (n - 1)) ) - 1.0;
out << p << '\t';
}
-
out << endl;
}
}else { done = true; }
}
+ if (numaxes > count) { m->mothurOut("You requested " + toString(numaxes) + " axes, but your file only includes " + toString(count) + ". Using " + toString(count) + "."); m->mothurOutEndLine(); numaxes = count; }
+
while (!in.eof()) {
if (m->control_pressed) { in.close(); return axes; }
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(axesfile, 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(globaldata->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) {
+ 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);
}
}
projects / mothur.git / blobdiff