#include "removeotuscommand.h"
#include "indicatorcommand.h"
#include "consensusseqscommand.h"
+#include "corraxescommand.h"
/*******************************************************/
commands["remove.otus"] = "remove.otus";
commands["indicator"] = "indicator";
commands["consensus.seqs"] = "consensus.seqs";
+ commands["corr.axes"] = "corr.axes";
commands["pairwise.seqs"] = "MPIEnabled";
commands["pipeline.pds"] = "MPIEnabled";
commands["classify.seqs"] = "MPIEnabled";
else if(commandName == "sub.sample") { command = new SubSampleCommand(optionString); }
else if(commandName == "indicator") { command = new IndicatorCommand(optionString); }
else if(commandName == "consensus.seqs") { command = new ConsensusSeqsCommand(optionString); }
+ else if(commandName == "corr.axes") { command = new CorrAxesCommand(optionString); }
else { command = new NoCommand(optionString); }
return command;
else if(commandName == "sub.sample") { pipecommand = new SubSampleCommand(optionString); }
else if(commandName == "indicator") { pipecommand = new IndicatorCommand(optionString); }
else if(commandName == "consensus.seqs") { pipecommand = new ConsensusSeqsCommand(optionString); }
+ else if(commandName == "corr.axes") { pipecommand = new CorrAxesCommand(optionString); }
else { pipecommand = new NoCommand(optionString); }
return pipecommand;
else if(commandName == "sub.sample") { shellcommand = new SubSampleCommand(); }
else if(commandName == "indicator") { shellcommand = new IndicatorCommand(); }
else if(commandName == "consensus.seqs") { shellcommand = new ConsensusSeqsCommand(); }
+ else if(commandName == "corr.axes") { shellcommand = new CorrAxesCommand(); }
else { shellcommand = new NoCommand(); }
return shellcommand;
//**********************************************************************************************************************
vector<string> CorrAxesCommand::getValidParameters(){
try {
- string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metastats","outputdir","inputdir"};
+ string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
return myArray;
}
else {
//valid paramters for this command
- string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metastats","outputdir","inputdir"};
+ string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"};
vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
OptionParser parser(option);
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["relabund"] = inputDir + it->second; }
}
+
+ it = parameters.find("metadata");
+ //user has given a template file
+ if(it != parameters.end()){
+ path = m->hasPath(it->second);
+ //if the user has not given a path then, add inputdir. else leave path alone.
+ if (path == "") { parameters["metadata"] = inputDir + it->second; }
+ }
}
else if (relabundfile == "not found") { relabundfile = ""; }
else { inputFileName = relabundfile; }
+ metadatafile = validParameter.validFile(parameters, "metadata", true);
+ if (metadatafile == "not open") { abort = true; }
+ else if (metadatafile == "not found") { metadatafile = ""; }
+
groups = validParameter.validFile(parameters, "groups", false);
if (groups == "not found") { groups = ""; pickedGroups = false; }
else {
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";
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";
+ for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; }
+ out << endl;
if (method == "pearson") { calcPearson(axes, out); }
- //else if (method == "spearman") { calcSpearman(axes, out); }
+ else if (method == "spearman") { calcSpearman(axes, out); }
//else if (method == "kendall") { calcKendal(axes, out); }
//else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); }
//for each otu
for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
+ out << i+1 << '\t';
+
//find the averages this otu - Y
float sumOtu = 0.0;
for (int j = 0; j < lookupFloat.size(); j++) {
}
float Ybar = sumOtu / (float) lookupFloat.size();
- //find the average
+ //find r value for each axis
+ for (int k = 0; k < averageAxes.size(); k++) {
+
+ double r = 0.0;
+ double numerator = 0.0;
+ double denomTerm1 = 0.0;
+ double denomTerm2 = 0.0;
+ for (int j = 0; j < lookupFloat.size(); j++) {
+ float Yi = lookupFloat[j]->getAbundance(i);
+ float Xi = axes[lookupFloat[j]->getGroup()][k];
+
+ numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
+ denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
+ denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
+ }
+
+ double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
+
+ r = numerator / denom;
+
+ out << r << '\t';
+ }
+
+ out << endl;
}
return 0;
}
}
//**********************************************************************************************************************
+int CorrAxesCommand::calcSpearman(map<string, vector<float> >& axes, ofstream& out) {
+ try {
+
+ //find average of each axis - X
+ vector<float> averageAxes; averageAxes.resize(numaxes, 0.0);
+ for (map<string, vector<float> >::iterator it = axes.begin(); it != axes.end(); it++) {
+ vector<float> temp = it->second;
+ for (int i = 0; i < temp.size(); i++) {
+ averageAxes[i] += temp[i];
+ }
+ }
+
+ for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); }
+
+ //for each otu
+ for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) {
+
+ out << i+1 << '\t';
+
+ //find the averages this otu - Y
+ float sumOtu = 0.0;
+ for (int j = 0; j < lookupFloat.size(); j++) {
+ sumOtu += lookupFloat[j]->getAbundance(i);
+ }
+ float Ybar = sumOtu / (float) lookupFloat.size();
+
+ //find r value for each axis
+ for (int k = 0; k < averageAxes.size(); k++) {
+
+ double r = 0.0;
+ double numerator = 0.0;
+ double denomTerm1 = 0.0;
+ double denomTerm2 = 0.0;
+ for (int j = 0; j < lookupFloat.size(); j++) {
+ float Yi = lookupFloat[j]->getAbundance(i);
+ float Xi = axes[lookupFloat[j]->getGroup()][k];
+
+ numerator += ((Xi - averageAxes[k]) * (Yi - Ybar));
+ denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k]));
+ denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
+ }
+
+ double denom = (sqrt(denomTerm1 * denomTerm2));
+
+ r = numerator / denom;
+
+ out << r << '\t';
+ }
+
+ out << endl;
+ }
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "CorrAxesCommand", "calcSpearman");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
int CorrAxesCommand::getShared(){
try {
InputData* input = new InputData(sharedfile, "sharedfile");
headerLine = headerLine.substr(pos+4);
}else { done = true; }
}
- cout << "count = " << count << endl;
while (!in.eof()) {
indexes.clear();
vector<Sequence*> seqsMatches;
+
vector<SeqDist> distsLeft;
vector<SeqDist> distsRight;
float distRight = distcalculator->getDist();
SeqDist subjectLeft;
- subjectLeft.seq = db[j];
+ subjectLeft.seq = NULL;
subjectLeft.dist = distLeft;
subjectLeft.index = j;
distsLeft.push_back(subjectLeft);
SeqDist subjectRight;
- subjectRight.seq = db[j];
+ subjectRight.seq = NULL;
subjectRight.dist = distRight;
subjectRight.index = j;
int lasti = 0;
for (int i = 0; i < distsLeft.size(); i++) {
//add left if you havent already
- it = seen.find(distsLeft[i].seq->getName());
+ it = seen.find(db[distsLeft[i].index]->getName());
if (it == seen.end()) {
dists.push_back(distsLeft[i]);
- seen[distsLeft[i].seq->getName()] = distsLeft[i].seq->getName();
+ seen[db[distsLeft[i].index]->getName()] = db[distsLeft[i].index]->getName();
lastLeft = distsLeft[i].dist;
}
//add right if you havent already
- it = seen.find(distsRight[i].seq->getName());
+ it = seen.find(db[distsRight[i].index]->getName());
if (it == seen.end()) {
dists.push_back(distsRight[i]);
- seen[distsRight[i].seq->getName()] = distsRight[i].seq->getName();
+ seen[db[distsRight[i].index]->getName()] = db[distsRight[i].index]->getName();
lastRight = distsRight[i].dist;
}
//cout << numWanted << endl;
for (int i = 0; i < numWanted; i++) {
//cout << dists[i].seq->getName() << '\t' << dists[i].dist << endl;
- Sequence* temp = new Sequence(dists[i].seq->getName(), dists[i].seq->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother.
+ Sequence* temp = new Sequence(db[dists[i].index]->getName(), db[dists[i].index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother.
seqsMatches.push_back(temp);
indexes.push_back(dists[i].index);
}
+ dists.clear(); distsLeft.clear(); distsRight.clear();
+
return seqsMatches;
}
catch(exception& e) {
map<int, int> trimmedPos;
//check to make sure that is not whole seq
if ((rearPos - frontPos - 1) <= 0) {
- m->mothurOut("[ERROR]: when I trim your sequences, the entire sequence is trimmed."); m->mothurOutEndLine();
+ m->mothurOut("[ERROR]: when I trim " + query->getName() + ", the entire sequence is trimmed. Skipping."); m->mothurOutEndLine();
query->setAligned("");
//trim topMatches
for (int i = 0; i < topMatches.size(); i++) {
projects / mothur.git / commitdiff