class Ace : public Calculator {
public:
- Ace(int n) : abund(n), Calculator("ace", 3) {};
+ Ace(int n) : abund(n), Calculator("ace", 3, false) {};
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
private:
int abund;
};
class BergerParker : public Calculator {
public:
- BergerParker() : Calculator("bergerparker", 1) {};
+ BergerParker() : Calculator("bergerparker", 1, false) {};
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
private:
};
class Bootstrap : public Calculator {
public:
- Bootstrap() : Calculator("Bootstrap", 1) {};
+ Bootstrap() : Calculator("Bootstrap", 1, false) {};
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
};
try {
int count = 1;
EstOutput data;
+ vector<SharedRAbundVector*> subset;
//if the users entered no valid calculators don't execute command
if (treeCalculators.size() == 0) { return 0; }
for (int k = 0; k < lookup.size(); k++) { // pass cdd each set of groups to commpare
for (int l = k; l < lookup.size(); l++) {
if (k != l) { //we dont need to similiarity of a groups to itself
+ subset.clear(); //clear out old pair of sharedrabunds
+ //add new pair of sharedrabunds
+ subset.push_back(lookup[k]); subset.push_back(lookup[l]);
+
//get estimated similarity between 2 groups
- data = treeCalculators[i]->getValues(lookup[k], lookup[l]); //saves the calculator outputs
+ data = treeCalculators[i]->getValues(subset); //saves the calculator outputs
//save values in similarity matrix
simMatrix[k][l] = data[0];
simMatrix[l][k] = data[0];
}
}
}
-
+
//creates tree from similarity matrix and write out file
createTree(out[i]);
}
//do merges and create tree structure by setting parents and children
//there are numGroups - 1 merges to do
for (int i = 0; i < (numGroups - 1); i++) {
-
- float largest = 0.0;
+
+ float largest = -1.0;
int row, column;
//find largest value in sims matrix by searching lower triangle
for (int j = 1; j < simMatrix.size(); j++) {
//set non-leaf node info and update leaves to know their parents
//non-leaf
t->tree[numGroups + i].setChildren(index[row], index[column]);
-
+
//parents
t->tree[index[row]].setParent(numGroups + i);
t->tree[index[column]].setParent(numGroups + i);
//branchlengths
t->tree[index[row]].setBranchLength(blength - t->tree[index[row]].getLengthToLeaves());
t->tree[index[column]].setBranchLength(blength - t->tree[index[column]].getLengthToLeaves());
-
+
//set your length to leaves to your childs length plus branchlength
t->tree[numGroups + i].setLengthToLeaves(t->tree[index[row]].getLengthToLeaves() + t->tree[index[row]].getBranchLength());
-
+
//update index
index[row] = numGroups+i;
index[column] = numGroups+i;
//zero out highest value that caused the merge.
- simMatrix[row][column] = 0.0;
- simMatrix[column][row] = 0.0;
-
+ simMatrix[row][column] = -1.0;
+ simMatrix[column][row] = -1.0;
+
//merge values in simsMatrix
for (int n = 0; n < simMatrix.size(); n++) {
//row becomes merge of 2 groups
simMatrix[row][n] = (simMatrix[row][n] + simMatrix[column][n]) / 2;
simMatrix[n][row] = simMatrix[row][n];
//delete column
- simMatrix[column][n] = 0.0;
- simMatrix[n][column] = 0.0;
+ simMatrix[column][n] = -1.0;
+ simMatrix[n][column] = -1.0;
}
}
-
+
//assemble tree
t->assembleTree();
-
+
//print newick file
t->print(*out);
-
+
//delete tree
delete t;
class BStick : public Calculator {
public:
- BStick() : Calculator("bstick", 3) {};
+ BStick() : Calculator("bstick", 3, false) {};
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
private:
double invSum(int, double);
public:
Calculator(){};
- Calculator(string n, int c) : name(n), cols(c) {};
+ Calculator(string n, int c, bool m) : name(n), cols(c), multiple(m) {};
virtual EstOutput getValues(SAbundVector*) = 0;
- virtual EstOutput getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) = 0;
+ virtual EstOutput getValues(vector<SharedRAbundVector*>) = 0;
virtual void print(ostream& f) { f.setf(ios::fixed, ios::floatfield); f.setf(ios::showpoint);
f << data[0]; for(int i=1;i<data.size();i++){ f << '\t' << data[i]; }}
virtual string getName() { return name; }
virtual int getCols() { return cols; }
+ virtual bool getMultiple() { return multiple; }
protected:
EstOutput data;
string name;
int cols;
+ bool multiple;
};
class Chao1 : public Calculator {
public:
- Chao1() : Calculator("Chao", 3) {};
+ Chao1() : Calculator("Chao", 3, false) {};
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
};
try {
globaldata = GlobalData::getInstance();
vector<SharedRAbundVector*> lookup;
+ vector<SharedRAbundVector*> subset;
//create and initialize vector of sharedvectors, one for each group
for (int i = 0; i < globaldata->Groups.size(); i++) {
//initialize labels for output
//makes 'uniqueAB uniqueAC uniqueBC' if your groups are A, B, C
getGroupComb();
- groupLabel = "";
+ groupLabel = "";
for (int s = 0; s < groupComb.size(); s++) {
groupLabel = groupLabel + label + groupComb[s] + "\t";
}
-
+
+ //for multi displays
+ string groupLabelAll = groupLabel + label + "all\t";
+
for(int i=0;i<displays.size();i++){
ccd->registerDisplay(displays[i]); //adds a display[i] to cdd
- displays[i]->init(groupLabel);
+ if (displays[i]->isCalcMultiple() == true) { displays[i]->init(groupLabelAll); }
+ else { displays[i]->init(groupLabel); }
}
//sample all the members
//calculate at 0 and the given increment
if((i == 0) || (i+1) % increment == 0){
- //randomize group order
- if (globaldata->getJumble() == "1") { random_shuffle(lookup.begin(), lookup.end()); }
- //how many comparisons to make i.e. for group a, b, c = ab, ac, bc.
+ //how many comparisons to make i.e. for group a, b, c = ab, ac, bc.
int n = 1;
for (int k = 0; k < (lookup.size() - 1); k++) { // pass cdd each set of groups to commpare
for (int l = n; l < lookup.size(); l++) {
- ccd->updateSharedData(lookup[k], lookup[l], i+1, globaldata->Groups.size());
+ subset.clear(); //clear out old pair of sharedrabunds
+ //add new pair of sharedrabund vectors
+ subset.push_back(lookup[k]); subset.push_back(lookup[l]);
+ ccd->updateSharedData(subset, i+1, globaldata->Groups.size());
}
n++;
}
+ //if this is a calculator that can do multiples then do them
+ ccd->updateSharedData(lookup, i+1, globaldata->Groups.size());
}
totalNumSeq = i+1;
}
int n = 1;
for (int k = 0; k < (lookup.size() - 1); k++) { // pass cdd each set of groups to commpare
for (int l = n; l < lookup.size(); l++) {
- ccd->updateSharedData(lookup[k], lookup[l], totalNumSeq, globaldata->Groups.size());
+ subset.clear(); //clear out old pair of sharedrabunds
+ //add new pair of sharedrabund vectors
+ subset.push_back(lookup[k]); subset.push_back(lookup[l]);
+ ccd->updateSharedData(subset, totalNumSeq, globaldata->Groups.size());
}
n++;
}
+ //if this is a calculator that can do multiples then do them
+ ccd->updateSharedData(lookup, totalNumSeq, globaldata->Groups.size());
}
//resets output files
output->output(nSeqs, data);
};
- void update(SharedRAbundVector* shared1, SharedRAbundVector* shared2, int numSeqs, int numGroups){
+ void update(vector<SharedRAbundVector*> shared, int numSeqs, int numGroups){
timesCalled++;
- data = estimate->getValues(shared1, shared2); //passes estimators a shared vector from each group to be compared
+ data = estimate->getValues(shared); //passes estimators a shared vector from each group to be compared
//figure out what groups are being compared in getValues
//because the jumble parameter randomizes the order we need to put the results in the correct column in the output file
int group1Index, group2Index, pos;
- group1Index = shared1->getGroupIndex();
- group2Index = shared2->getGroupIndex();
+ group1Index = shared[0]->getGroupIndex();
+ group2Index = shared[1]->getGroupIndex();
numGroupComb = 0;
int n = 1;
}
n++;
}
-
- groupData.resize((numGroupComb*data.size()), 0);
- //fills groupdata with datas info
- for (int i = 0; i < data.size(); i++) {
- groupData[pos+i] = data[i];
+ if (estimate->getMultiple() == true) {
+ numGroupComb++;
+ groupData.resize((numGroupComb*data.size()), 0);
+ //is this the time its called with all values
+ if ((timesCalled % numGroupComb) == 0) {
+ //last spot
+ pos = ((groupData.size()-1) * data.size());
+ }
+ //fills groupdata with datas info
+ for (int i = 0; i < data.size(); i++) {
+ groupData[pos+i] = data[i];
+ }
+ }else {
+ groupData.resize((numGroupComb*data.size()), 0);
+ //fills groupdata with datas info
+ for (int i = 0; i < data.size(); i++) {
+ groupData[pos+i] = data[i];
+ }
}
//when you get all your groups info then output
void init(string s) { output->initFile(s); };
void reset() { output->resetFile(); };
void close() { output->resetFile(); };
+ bool isCalcMultiple() { return estimate->getMultiple(); }
private:
Calculator* estimate;
class SharedCollectorsCurveData : public Observable {
public:
- SharedCollectorsCurveData() : shared1(0), shared2(0) {};
+ SharedCollectorsCurveData() { }; //: shared1(0), shared2(0)
void registerDisplay(Display* o) { displays.insert(o); };
void removeDisplay(Display* o) { displays.erase(o); delete o; };
- void SharedDataChanged() { notifyDisplays(); };
- void updateSharedData(SharedRAbundVector* rv, SharedRAbundVector* rv2, int numSeqs, int numGroupComb) { shared1 = rv; shared2 = rv2; NumSeqs = numSeqs; NumGroupComb = numGroupComb; SharedDataChanged(); };
-
+ void SharedDataChanged() { notifyDisplays(); };
+ void updateSharedData(vector<SharedRAbundVector*> s, int numSeqs, int numGroupComb) { shared = s; NumSeqs = numSeqs; NumGroupComb = numGroupComb; SharedDataChanged(); };
+
void notifyDisplays(){
for(set<Display*>::iterator pos=displays.begin();pos!=displays.end();pos++){
- (*pos)->update(shared1, shared2, NumSeqs, NumGroupComb);
+ if ( ((*pos)->isCalcMultiple() == true) || ( ((*pos)->isCalcMultiple() == false) && (shared.size() == 2) ) ) {
+ (*pos)->update(shared, NumSeqs, NumGroupComb);
+ }
}
};
private:
set<Display*> displays;
- SharedRAbundVector* shared1;
- SharedRAbundVector* shared2;
+ vector<Display*> multiDisplays;
+ vector<SharedRAbundVector*> shared;
int NumSeqs, NumGroupComb;
};
#include "getoturepcommand.h"
#include "treegroupscommand.h"
#include "bootstrapsharedcommand.h"
+#include "concensuscommand.h"
/***********************************************************/
else if(commandName == "get.oturep") { command = new GetOTURepCommand(); }
else if(commandName == "tree.shared") { command = new TreeGroupCommand(); }
else if(commandName == "bootstrap.shared") { command = new BootSharedCommand(); }
+ else if(commandName == "concensus") { command = new ConcensusCommand(); }
else { command = new NoCommand(); }
return command;
class Coverage : public Calculator {
public:
- Coverage() : Calculator("coverage", 1) {};
+ Coverage() : Calculator("coverage", 1, false) {};
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
};
public:
virtual void update(SAbundVector* rank) = 0;
- virtual void update(SharedRAbundVector* shared1, SharedRAbundVector* shared2, int numSeqs, int numGroupComb) = 0;
+ virtual void update(vector<SharedRAbundVector*> shared, int numSeqs, int numGroupComb) = 0;
virtual void init(string) = 0;
virtual void reset() = 0;
virtual void close() = 0;
+ virtual bool isCalcMultiple() = 0;
private:
}
}
- if ((commandName == "unifrac.weighted") || (commandName == "unifrac.unweighted")) {
+ if ((commandName == "unifrac.weighted") || (commandName == "unifrac.unweighted") || (commandName == "concensus")) {
if (globaldata->gTree.size() == 0) {//no trees were read
- cout << "You must execute the read.tree command, before you may execute the unifrac.weighted or unifrac.unweighted command." << endl; return false; }
+ cout << "You must execute the read.tree command, before you may execute the unifrac.weighted, unifrac.unweighted or concensus command." << endl; return false; }
}
//check for valid method
renameOk = rename(outName.c_str(), inName.c_str());
//checks to make sure user was able to rename and remove successfully
- if ((renameOk != 0)) { cout << "Unable to rename necessary files." << endl; cout << outName << " g " << inName << endl;}
+ if ((renameOk != 0)) { cout << "Unable to rename necessary files." << endl; }
}
catch(exception& e) {
class ThreeColumnFile : public FileOutput {
public:
- ThreeColumnFile(string n) : FileOutput(), inName(n), counter(0), outName(getPathName(n) + ".temp." + getSimpleName(n)) { };
+ ThreeColumnFile(string n) : FileOutput(), inName(n), counter(0), outName(getPathName(n) + ".temp") { };
~ThreeColumnFile();
void initFile(string);
void output(int, vector<double>);
public:
- OneColumnFile(string n) : inName(n), counter(0), outName(getPathName(n) + ".temp." + getSimpleName(n)) {};
+ OneColumnFile(string n) : inName(n), counter(0), outName(getPathName(n) + ".temp") {};
~OneColumnFile();
void output(int, vector<double>);
void initFile(string);
public:
- SharedOneColumnFile(string n) : inName(n), counter(0), outName(getPathName(n) + ".temp." + getSimpleName(n)) {};
+ SharedOneColumnFile(string n) : inName(n), counter(0), outName(getPathName(n) + ".temp") {};
~SharedOneColumnFile();
void output(int, vector<double>);
void initFile(string);
class SharedThreeColumnFile : public FileOutput {
public:
- SharedThreeColumnFile(string n, string groups) : FileOutput(), groupLabel(groups), inName(n), counter(0), numGroup(1), outName(getPathName(n) + ".temp." + getSimpleName(n)) { };
+ SharedThreeColumnFile(string n, string groups) : FileOutput(), groupLabel(groups), inName(n), counter(0), numGroup(1), outName(getPathName(n) + ".temp") { };
~SharedThreeColumnFile();
void initFile(string);
void output(int, vector<double>);
};
/***********************************************************************/
-
+//used by parsimony, unifrac.weighted and unifrac.unweighted
class ColumnFile : public FileOutput {
public:
- ColumnFile(string n) : FileOutput(), inName(n), counter(0), outName(getPathName(n) + ".temp." + getSimpleName(n)) { globaldata = GlobalData::getInstance(); };
+ ColumnFile(string n) : FileOutput(), inName(n), counter(0), outName(getPathName(n) + ".temp") { globaldata = GlobalData::getInstance(); };
~ColumnFile();
//to make compatible with parent class
class Geom : public Calculator {
public:
- Geom() : Calculator("geometric", 3) {};
+
+ Geom() : Calculator("geometric", 3, false) {};
+
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
private:
double kEq(double, double);
cout << "Note: No spaces between parameter labels (i.e. freq), '=' and parameters (i.e.yourFreq)." << "\n" << "\n";
}else if (globaldata->helpRequest == "collect.shared") {
cout << "The collect.shared command can only be executed after a successful read.otu command." << "\n";
- cout << "The collect.shared command parameters are label, line, freq, jumble, calc and groups. No parameters are required, but you may not use " << "\n";
+ cout << "The collect.shared command parameters are label, line, freq, calc and groups. No parameters are required, but you may not use " << "\n";
cout << "both the line and label parameters at the same time. The collect.shared command should be in the following format: " << "\n";
- cout << "collect.shared(label=yourLabel, line=yourLines, freq=yourFreq, jumble=yourJumble, calc=yourEstimators, groups=yourGroups)." << "\n";
- cout << "Example collect.shared(label=unique-.01-.03, line=0-5-10, freq=10, jumble=1, groups=B-C, calc=sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan)." << "\n";
- cout << "The default values for jumble is 1 (meaning jumble, if it’s set to 0 then it will not jumble), freq is 100 and calc are sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan." << "\n";
+ cout << "collect.shared(label=yourLabel, line=yourLines, freq=yourFreq, calc=yourEstimators, groups=yourGroups)." << "\n";
+ cout << "Example collect.shared(label=unique-.01-.03, line=0-5-10, freq=10, groups=B-C, calc=sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan)." << "\n";
+ cout << "The default values for freq is 100 and calc are sharedsobs-sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan." << "\n";
cout << "The default value for groups is all the groups in your groupfile." << "\n";
validCalcs->printCalc("shared", cout);
cout << "The label and line parameters are used to analyze specific lines in your input." << "\n";
cout << "Note: No spaces between parameter labels (i.e. line), '=' and parameters (i.e.yourLines)." << "\n" << "\n";
}else if (globaldata->helpRequest == "summary.shared") {
cout << "The summary.shared command can only be executed after a successful read.otu command." << "\n";
- cout << "The summary.shared command parameters are label, line, jumble and calc. No parameters are required, but you may not use " << "\n";
+ cout << "The summary.shared command parameters are label, line and calc. No parameters are required, but you may not use " << "\n";
cout << "both the line and label parameters at the same time. The summary.shared command should be in the following format: " << "\n";
- cout << "summary.shared(label=yourLabel, line=yourLines, jumble=yourJumble, calc=yourEstimators, groups=yourGroups)." << "\n";
- cout << "Example summary.shared(label=unique-.01-.03, line=0,5,10, jumble=1, groups=B-C, calc=sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan)." << "\n";
+ cout << "summary.shared(label=yourLabel, line=yourLines, calc=yourEstimators, groups=yourGroups)." << "\n";
+ cout << "Example summary.shared(label=unique-.01-.03, line=0,5,10, groups=B-C, calc=sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan)." << "\n";
validCalcs->printCalc("sharedsummary", cout);
- cout << "The default value for jumble is 1 (meaning jumble, if it’s set to 0 then it will not jumble) and calc is sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan" << "\n";
+ cout << "The default value for calc is sharedsobs-sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan" << "\n";
cout << "The default value for groups is all the groups in your groupfile." << "\n";
cout << "The label and line parameters are used to analyze specific lines in your input." << "\n";
cout << "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed. You must enter at least 2 valid groups." << "\n";
validCalcs->printCalc("boot", cout);
cout << "The bootstrap.shared command outputs a .tre file for each calculator you specify at each distance you choose containing iters number of trees." << "\n";
cout << "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups)." << "\n" << "\n";
+ }else if (globaldata->helpRequest == "concensus") {
+ cout << "The concensus command can only be executed after a successful read.tree command." << "\n";
+ cout << "The concensus command has no parameters." << "\n";
+ cout << "The concensus command should be in the following format: concensus()." << "\n";
+ cout << "The concensus command output two files: .concensus.tre and .concensuspairs." << "\n";
+ cout << "The .concensus.tre file contains the concensus tree of the trees in your input file." << "\n";
+ cout << "The branch lengths are the percentage of trees in your input file that had the given pair." << "\n";
+ cout << "The .concensuspairs file contains a list of the internal nodes in your tree. For each node, the pair that was used in the concensus tree " << "\n";
+ cout << "is reported with its percentage, as well as the other pairs that were seen for that node but not used and their percentages." << "\n" << "\n";
}else if (globaldata->helpRequest == "bin.seqs") {
cout << "The bin.seqs command can only be executed after a successful read.otu command of a list file." << "\n";
cout << "The bin.seqs command parameters are fasta, name, line and label. The fasta parameter is required, and you may not use line and label at the same time." << "\n";
class Jackknife : public Calculator {
public:
- Jackknife() : Calculator("Jackknife", 3) { getAMatrix(); };
+ Jackknife() : Calculator("Jackknife", 3, false) { getAMatrix(); };
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
private:
static const int maxOrder = 30;
class LogSD : public Calculator {
public:
- LogSD() : Calculator("logseries", 3) {};
+
+ LogSD() : Calculator("logseries", 3, false) {};
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) { return data; };
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
private:
double logS(double);
typedef unsigned long long ull;
-
+struct IntNode {
+ int lvalue;
+ int rvalue;
+ int lcoef;
+ int rcoef;
+ IntNode* left;
+ IntNode* right;
+};
+
/***********************************************************************/
// snagged from http://www.parashift.com/c++-faq-lite/misc-technical-issues.html#faq-39.2
return output.str();
}
-
//**********************************************************************************************************************
template<typename T>
return simpleName;
}
+/***********************************************************************/
+
+inline int factorial(int num){
+ int total = 1;
+
+ for (int i = 1; i <= num; i++) {
+ total *= i;
+ }
+
+ return total;
+}
/***********************************************************************/
class NPShannon : public Calculator {
public:
- NPShannon() : Calculator("NPShannon", 1) {};
+ NPShannon() : Calculator("NPShannon", 1, false) {};
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
private:
};
class NSeqs : public Calculator {
public:
- NSeqs() : Calculator("NSeqs", 1) {};
+ NSeqs() : Calculator("NSeqs", 1, false) {};
EstOutput getValues(SAbundVector* rank){
data.resize(1,0);
data[0] = (double)rank->getNumSeqs();
return data;
}
- EstOutput getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
};
/***********************************************************************/
class QStat : public Calculator {
public:
- QStat() : Calculator("qstat", 1) {};
+ QStat() : Calculator("qstat", 1, false) {};
+
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
private:
RAbundVector rdata;
};
/***********************************************************************/
-void RareDisplay::update(SharedRAbundVector* shared1, SharedRAbundVector* shared2, int numSeqs, int numGroupComb) {
+void RareDisplay::update(vector<SharedRAbundVector*> shared, int numSeqs, int numGroupComb) {
try {
- vector<double> data = estimate->getValues(shared1, shared2);
+ vector<double> data = estimate->getValues(shared);
double newNSeqs = data[0];
if(nIters != 1){
public:
RareDisplay(Calculator* calc, FileOutput* file) : estimate(calc), output(file), nIters(1),
- tempInName(getPathName(output->getFileName()) + ".tempin."+ getSimpleName(output->getFileName())), tempOutName(getPathName(output->getFileName()) + ".tempout."+ getSimpleName(output->getFileName())) {};
+ tempInName(getPathName(output->getFileName()) + ".tempin"), tempOutName(getPathName(output->getFileName()) + ".tempout") {};
~RareDisplay() { delete estimate; delete output; };
void init(string);
void reset();
void update(SAbundVector*);
- void update(SharedRAbundVector* shared1, SharedRAbundVector* shared2, int numSeqs, int numGroupComb);
+ void update(vector<SharedRAbundVector*> shared, int numSeqs, int numGroupComb);
void close();
+ bool isCalcMultiple() { return estimate->getMultiple(); }
private:
Calculator* estimate;
//randomize the groups
random_shuffle(lookup.begin(), lookup.end());
-
+
+ vector<SharedRAbundVector*> subset;
//send each group one at a time
for (int k = 0; k < lookup.size(); k++) {
- rcd->updateSharedData(lookup[0], lookup[k], k+1, numGroupComb);
+ subset.clear(); //clears out old pair of sharedrabunds
+ //add in new pair of sharedrabunds
+ subset.push_back(lookup[0]); subset.push_back(lookup[k]);
+
+ rcd->updateSharedData(subset, k+1, numGroupComb);
mergeVectors(lookup[0], lookup[k]);
}
class SharedRarefactionCurveData : public Observable {
public:
- SharedRarefactionCurveData() : shared1(0), shared2(0) {};
+ SharedRarefactionCurveData() {}; //: shared1(0), shared2(0)
void registerDisplay(Display* o) { displays.insert(o); };
void removeDisplay(Display* o) { displays.erase(o); delete o; };
void SharedDataChanged() { notifyDisplays(); };
- void updateSharedData(SharedRAbundVector* rv, SharedRAbundVector* rv2, int numSeqs, int numGroupComb) { shared1 = rv; shared2 = rv2; NumSeqs = numSeqs; NumGroupComb = numGroupComb; SharedDataChanged(); };
+ void updateSharedData(vector<SharedRAbundVector*> r, int numSeqs, int numGroupComb) { shared = r; NumSeqs = numSeqs; NumGroupComb = numGroupComb; SharedDataChanged(); };
void notifyDisplays(){
for(set<Display*>::iterator pos=displays.begin();pos!=displays.end();pos++){
- (*pos)->update(shared1, shared2, NumSeqs, NumGroupComb);
+ (*pos)->update(shared, NumSeqs, NumGroupComb);
}
};
private:
set<Display*> displays;
- SharedRAbundVector* shared1;
- SharedRAbundVector* shared2;
+ vector<SharedRAbundVector*> shared;
int NumSeqs, NumGroupComb;
};
class Shannon : public Calculator {
public:
- Shannon() : Calculator("Shannon", 3) {};
+ Shannon() : Calculator("Shannon", 3, false) {};
EstOutput getValues(SAbundVector* rank);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
private:
};
/***********************************************************************/
-EstOutput SharedAce::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput SharedAce::getValues(vector<SharedRAbundVector*> shared) {
try {
data.resize(1,0);
string label;
- label = shared1->getLabel();
+ label = shared[0]->getLabel();
int fARare, fBRare, S12Rare, S12Abund, S12, f11, tempA, tempB, t10, t01, t11, t21, t12, t22, C12Numerator;
fARare = 0; fBRare = 0; S12Rare = 0; S12Abund = 0; S12 = 0; f11 = 0; t10 = 0; t01 = 0; t11= 0; t21= 0; t12= 0; t22= 0; C12Numerator = 0;
- float SharedAce, C12, part1, part2, part3, part4, part5, Gamma1, Gamma2, Gamma3;
+ float Sharedace, C12, part1, part2, part3, part4, part5, Gamma1, Gamma2, Gamma3;
/*fARare = number of OTUs with one individual found in A and less than or equal to 10 in B.
fBRare = number of OTUs with one individual found in B and less than or equal to 10 in A.
S12 = number of shared OTUs in A and B
This estimator was changed to reflect Caldwell's changes, eliminating the nrare / nrare - 1 */
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
if ((tempA != 0) && (tempB != 0)) {//they are shared
S12++;
//do both A and B have one
part4 = fBRare * Gamma2;
part5 = f11 * Gamma3;
- SharedAce = S12Abund + part1 + (part2 * (part3 + part4 + part5));
- data[0] = SharedAce;
+ Sharedace = S12Abund + part1 + (part2 * (part3 + part4 + part5));
+ data[0] = Sharedace;
return data;
}
class SharedAce : public Calculator {
public:
- SharedAce(int n=10) : abund(n), Calculator("sharedace", 3) {};
+ SharedAce(int n=10) : abund(n), Calculator("sharedace", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
int abund;
};
/***********************************************************************/
-EstOutput Anderberg::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput Anderberg::getValues(vector<SharedRAbundVector*> shared) {
try {
int S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
data.resize(1,0);
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
if (tempA != 0) { S1++; }
if (tempB != 0) { S2++; }
class Anderberg : public Calculator {
public:
- Anderberg() : Calculator("Anderberg", 1) {};
+ Anderberg() : Calculator("Anderberg", 1, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
#include "sharedbraycurtis.h"
/***********************************************************************/
-
-EstOutput BrayCurtis::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+//This is used by SharedJAbund and SharedSorAbund
+EstOutput BrayCurtis::getValues(vector<SharedRAbundVector*> shared) {
try {
data.resize(1,0);
sumSharedAB = the sum of the minimum otus int all shared otus in AB.
*/
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
sumSharedA += tempA;
sumSharedB += tempB;
class BrayCurtis : public Calculator {
public:
- BrayCurtis() : Calculator("BrayCurtis", 1) {};
+ BrayCurtis() : Calculator("BrayCurtis", 1, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
#include "sharedchao1.h"
/***********************************************************************/
-EstOutput SharedChao1::getValues(SharedRAbundVector* sharedA, SharedRAbundVector* sharedB){
+EstOutput SharedChao1::getValues(vector<SharedRAbundVector*> shared){
try {
data.resize(1,0);
- int f11, f1A, f2A, f1B, f2B, S12, tempA, tempB;
- f11 = 0; f1A = 0; f2A = 0; f1B = 0; f2B = 0; S12 = 0;
- float ChaoAB, part1, part2, part3;
-
- /* f11 = number of shared OTUs with one observed individual in A and B
- f1A, f2A = number of shared OTUs with one or two individuals observed in A
- f1B, f2B = number of shared OTUs with one or two individuals observed in B
- S12 = number of shared OTUs in A and B */
-
- //loop through vectors calculating the f11, f1A, f2A, f1B, f2B, S12 values
- for (int i = 0; i< sharedA->size(); i++) {
- tempA = sharedA->getAbundance(i); //store in temps to avoid calling getAbundance multiple times
- tempB = sharedB->getAbundance(i);
- if ((tempA != 0) && (tempB != 0)) {//they are shared
- S12++; //they are shared
- //do both A and B have one
- if ((tempA == 1) && (tempB == 1)) { f11++; }
-
- //does A have one or two
- if (tempA == 1) { f1A++; }
- else if (tempA == 2) { f2A++; }
+ vector<int> temp;
+ int numGroups = shared.size();
+ float Chao = 0.0; float leftvalue, rightvalue;
- //does B have one or two
- if (tempB == 1) { f1B++; }
- else if (tempB == 2) { f2B++; }
-
+ // IntNode is defined in mothur.h
+ // The tree used here is a binary tree used to represent the f1+++, f+1++, f++1+, f+++1, f11++, f1+1+...
+ // combinations required to solve the chao estimator equation for any number of groups. Conceptually, think
+ // of each node as having a 1 and a + value, or for f2 values a 2 and a + value, and 2 pointers to intnodes, and 2 coeffient values.
+ // The coeffient value is how many times you chose branch 1 to get to that fvalue.
+ // If you choose left you are selecting the 1 or 2 value and right means the + value. For instance, to find
+ // the number of bins that have f1+1+ you would start at the root, go left, right, left, and select the rightvalue.
+ // the coeffient is 2. Note: we only set the coeffient in f2 values.
+
+ //create and initialize trees to 0.
+ initialTree(numGroups);
+
+ //loop through vectors calculating the f11, f1A, f2A, f1B, f2B, S12 values
+ for (int i = 0; i < shared[0]->size(); i++) {
+ //get bin values and calc shared
+ bool sharedByAll = true;
+ temp.clear();
+ for (int j = 0; j < numGroups; j++) {
+ temp.push_back(shared[j]->getAbundance(i));
+ if (temp[j] == 0) { sharedByAll = false; }
+ }
+
+ //they are shared
+ if (sharedByAll == true) {
+
+ // cout << "temp = ";
+ // for (int h = 0; h < temp.size(); h++) { cout << temp[h] << " "; } cout << endl;
+ //find f1 and f2values
+ updateTree(temp);
}
}
-
- //checks for divide by zero error
- if ((f2A == 0) || (f2B == 0)) {
- part1 = ((float)(f1A*f1B)/(float)(4*(f2A+1)*(f2B+1)));
- part2 = ((float)(f1A*(f1A-1))/(float)(2*f2A+2));
- part3 = ((float)(f1B*(f1B-1))/(float)(2*f2B+2));
- }else {
- part1 = ((float)(f1A*f1B)/(float)(4*f2A*f2B));
- part2 = ((float)(f1A*f1A)/(float)(2*f2A));
- part3 = ((float)(f1B*f1B)/(float)(2*f2B));
+
+
+ //cout << "Entering " << endl;
+ //calculate chao1, (numleaves-1) because numleaves contains the ++ values.
+ for (int i = 0; i < numLeaves; i++) {
+ //divide by zero error
+ if (f2leaves[i]->lvalue == 0) { f2leaves[i]->lvalue++; }
+ if (f2leaves[i]->rvalue == 0) { f2leaves[i]->rvalue++; }
+
+ //cout << "f2 leaves ";
+ //cout << f2leaves[i]->lvalue << f2leaves[i]->rvalue << endl;
+
+ // cout << "f2 leaf coef ";
+ //cout << f2leaves[i]->lcoef << '\t' << f2leaves[i]->rcoef << endl;
+
+ // cout << "f1 leaves ";
+ // cout << f1leaves[i]->lvalue << f1leaves[i]->rvalue << endl;
+
+
+ leftvalue = (float)(f1leaves[i]->lvalue * f1leaves[i]->lvalue) / (float)((pow(2, (float)f2leaves[i]->lcoef)) * f2leaves[i]->lvalue);
+ if (i != (numLeaves-1)) {
+ rightvalue = (float)(f1leaves[i]->rvalue * f1leaves[i]->rvalue) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * f2leaves[i]->rvalue);
+ }else{
+ rightvalue = (float)(f1leaves[i]->rvalue);
+ }
+ //cout << "left = " << leftvalue << " right = " << rightvalue << endl;
+ Chao += leftvalue + rightvalue;
}
- ChaoAB = (float)S12 + (float)(f11*part1) + part2 + part3;
- data[0] = ChaoAB;
-
+ // cout << "exiting " << endl;
+ data[0] = Chao;
return data;
}
catch(exception& e) {
cout << "An unknown error has occurred in the SharedChao1 class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
+}
+
+/***********************************************************************/
+//builds trees structure with n leaf nodes initialized to 0.
+void SharedChao1::initialTree(int n) {
+ try {
+ // (2^n) / 2. Divide by 2 because each leaf node contains 2 values. One for + and one for 1 or 2.
+ numLeaves = pow(2, (float)n) / 2;
+ numNodes = 2*numLeaves - 1;
+ int countleft = 0;
+ int countright = 1;
+
+ f1leaves.resize(numNodes);
+ f2leaves.resize(numNodes);
+
+ //initialize leaf values
+ for (int i = 0; i < numLeaves; i++) {
+ f1leaves[i] = new IntNode;
+ f1leaves[i]->lvalue = 0;
+ f1leaves[i]->rvalue = 0;
+ f1leaves[i]->left = NULL;
+ f1leaves[i]->right = NULL;
+
+ f2leaves[i] = new IntNode;
+ f2leaves[i]->lvalue = 0;
+ f2leaves[i]->rvalue = 0;
+ f2leaves[i]->left = NULL;
+ f2leaves[i]->right = NULL;
+ }
+
+ //set pointers to children
+ for (int j = numLeaves; j < numNodes; j++) {
+ f1leaves[j] = new IntNode;
+ f1leaves[j]->left = f1leaves[countleft];
+ f1leaves[j]->right = f1leaves[countright];
+
+ f2leaves[j] = new IntNode;
+ f2leaves[j]->left = f2leaves[countleft];
+ f2leaves[j]->right =f2leaves[countright];
+
+ countleft = countleft + 2;
+ countright = countright + 2;
+ }
+
+ //point to root
+ f1root = f1leaves[numNodes-1];
+
+ //point to root
+ f2root = f2leaves[numNodes-1];
+
+ //set coeffients
+ setCoef(f2root, 0);
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function initialTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the SharedChao1 class function initialTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+/***********************************************************************/
+//take vector containing the abundance info. for a bin and updates trees.
+void SharedChao1::updateTree(vector<int> bin) {
+ try {
+ updateBranchf1(f1root, bin, 0);
+ updateBranchf2(f2root, bin, 0);
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the SharedChao1 class function updateTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
}
/***********************************************************************/
+void SharedChao1::updateBranchf1(IntNode* node, vector<int> bin, int index) {
+ try {
+ //if you have more than one group
+ if (index == (bin.size()-1)) {
+ if (bin[index] == 1) { node->lvalue++; node->rvalue++; }
+ else { node->rvalue++; }
+ }else {
+ if (bin[index] == 1) {
+ //follow path as if you are 1
+ updateBranchf1(node->left, bin, index+1);
+ }
+ //follow path as if you are +
+ updateBranchf1(node->right, bin, index+1);
+ }
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateBranchf1. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the SharedChao1 class function updateBranchf1. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+
+/***********************************************************************/
+void SharedChao1::updateBranchf2(IntNode* node, vector<int> bin, int index) {
+ try {
+ //if you have more than one group
+ if (index == (bin.size()-1)) {
+ if (bin[index] == 2) { node->lvalue++; node->rvalue++; }
+ else { node->rvalue++; }
+ }else {
+ if (bin[index] == 2) {
+ //follow path as if you are 1
+ updateBranchf2(node->left, bin, index+1);
+ }
+ //follow path as if you are +
+ updateBranchf2(node->right, bin, index+1);
+ }
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateBranchf2. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the SharedChao1 class function updateBranchf2. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+
+/***********************************************************************/
+void SharedChao1::setCoef(IntNode* node, int coef) {
+ try {
+ if (node->left != NULL) {
+ setCoef(node->left, coef+1);
+ setCoef(node->right, coef);
+ }else {
+ node->lcoef = coef+1;
+ node->rcoef = coef;
+ }
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function getCoef. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the SharedChao1 class function getCoef. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+
+/***********************************************************************/
+//for debugging purposes
+void SharedChao1::printTree() {
+
+ cout << "F1 leaves" << endl;
+ printBranch(f1root);
+
+ cout << "F2 leaves" << endl;
+ printBranch(f2root);
+
+
+}
+/*****************************************************************/
+void SharedChao1::printBranch(IntNode* node) {
+ try {
+
+ // you are not a leaf
+ if (node->left != NULL) {
+ printBranch(node->left);
+ printBranch(node->right);
+ }else { //you are a leaf
+ cout << node->lvalue << endl;
+ cout << node->rvalue << endl;
+ }
+
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function printBranch. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the Tree class function printBranch. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+
+/*****************************************************************/
+
+
+
+
/***********************************************************************/
+
class SharedChao1 : public Calculator {
-public:
- SharedChao1() : Calculator("sharedchao", 3) {};
- EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ public:
+ SharedChao1() : Calculator("sharedchao", 3, true) {};
+ EstOutput getValues(SAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>);
+ private:
+ IntNode* f1root;
+ IntNode* f2root;
+ vector<IntNode*> f1leaves;
+ vector<IntNode*> f2leaves;
+ int numLeaves;
+ int numNodes;
+
+ void initialTree(int); //builds trees structure with n leaf nodes initialized to 0.
+ void setCoef(IntNode*, int);
+ void updateTree(vector<int>); //take vector containing the abundance info. for a bin and updates trees.
+ void updateBranchf1(IntNode*, vector<int>, int); //pointer, vector of abundance values, index into vector
+ void updateBranchf2(IntNode*, vector<int>, int); //pointer, vector of abundance values, index into vector
+
+ //for debugging
+ void printTree();
+ void printBranch(IntNode*);
};
/***********************************************************************/
/***********************************************************************/
-EstOutput JAbund::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput JAbund::getValues(vector<SharedRAbundVector*> shared) {
try {
EstOutput UVest;
UVest.resize(2,0);
data.resize(1,0);
- UVest = uv->getUVest(shared1, shared2);
+ UVest = uv->getUVest(shared);
//UVest[0] is Uest UVest[1] is Vest
data[0] = (UVest[0] * UVest[1]) / ((float)(UVest[0] + UVest[1] - (UVest[0] * UVest[1])));
class JAbund : public Calculator {
public:
- JAbund() : Calculator("JAbund", 3) {};
+ JAbund() : Calculator("JAbund", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
UVEst* uv;
/***********************************************************************/
-EstOutput Jclass::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput Jclass::getValues(vector<SharedRAbundVector*> shared) {
try {
int S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
data.resize(1,0);
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
//find number of bins in shared1 and shared2
if (tempA != 0) { S1++; }
class Jclass : public Calculator {
public:
- Jclass() : Calculator("Jclass", 3) {};
+ Jclass() : Calculator("Jclass", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
/***********************************************************************/
-EstOutput Jest::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput Jest::getValues(vector<SharedRAbundVector*> shared) {
try {
EstOutput S1, S2, S12;
S12.resize(1,0);
SAbundVector* chaoS1Sabund = new SAbundVector();
SAbundVector* chaoS2Sabund = new SAbundVector();
- *chaoS1Sabund = shared1->getSAbundVector();
- *chaoS2Sabund = shared2->getSAbundVector();
+ *chaoS1Sabund = shared[0]->getSAbundVector();
+ *chaoS2Sabund = shared[1]->getSAbundVector();
- S12 = sharedChao->getValues(shared1, shared2);
+ S12 = sharedChao->getValues(shared);
S1 = chaoS1->getValues(chaoS1Sabund);
S2 = chaoS2->getValues(chaoS2Sabund);
class Jest : public Calculator {
public:
- Jest() : Calculator("Jest", 3) {};
+ Jest() : Calculator("Jest", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
/***********************************************************************/
-EstOutput KSTest::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
+EstOutput KSTest::getValues(vector<SharedRAbundVector*> shared){
try {
data.resize(2,0);
//Must return shared1 and shared2 to original order at conclusion of kstest
- vector <individual> initData1 = shared1->getData();
- vector <individual> initData2 = shared2->getData();
- shared1->sortD();
- shared2->sortD();
+ vector <individual> initData1 = shared[0]->getData();
+ vector <individual> initData2 = shared[1]->getData();
+ shared[0]->sortD();
+ shared[1]->sortD();
- int numNZ1 = shared1->numNZ();
- int numNZ2 = shared2->numNZ();
- double numInd1 = (double)shared1->getNumSeqs();
- double numInd2 = (double)shared2->getNumSeqs();
+ int numNZ1 = shared[0]->numNZ();
+ int numNZ2 = shared[1]->numNZ();
+ double numInd1 = (double)shared[0]->getNumSeqs();
+ double numInd2 = (double)shared[1]->getNumSeqs();
double maxDiff = -1;
double sum1 = 0;
double sum2 = 0;
- for(int i = 1; i < shared1->getNumBins(); i++)
+ for(int i = 1; i < shared[0]->getNumBins(); i++)
{
- sum1 += shared1->get(i).abundance;
- sum2 += shared2->get(i).abundance;
+ sum1 += shared[0]->get(i).abundance;
+ sum2 += shared[1]->get(i).abundance;
double diff = fabs((double)sum1/numInd1 - (double)sum2/numInd2);
if(diff > maxDiff)
maxDiff = diff;
cout << "If D-Statistic is greater than the critical value then the data sets are significantly different at the 95% confidence level.\n\n";
}*/
- shared1->setData(initData1);
- shared2->setData(initData2);
+ shared[0]->setData(initData1);
+ shared[1]->setData(initData2);
data[0] = DStatistic;
data[1] = critVal;
class KSTest : public Calculator {
public:
- KSTest() : Calculator("kstest", 3) {};
+ KSTest() : Calculator("kstest", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
/***********************************************************************/
-EstOutput Kulczynski::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput Kulczynski::getValues(vector<SharedRAbundVector*> shared) {
try {
int S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
data.resize(1,0);
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
if (tempA != 0) { S1++; }
if (tempB != 0) { S2++; }
class Kulczynski : public Calculator {
public:
- Kulczynski() : Calculator("Kulczynski", 1) {};
+ Kulczynski() : Calculator("Kulczynski", 1, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
/***********************************************************************/
-EstOutput KulczynskiCody::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput KulczynskiCody::getValues(vector<SharedRAbundVector*> shared) {
try {
int S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
data.resize(1,0);
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
if (tempA != 0) { S1++; }
if (tempB != 0) { S2++; }
class KulczynskiCody : public Calculator {
public:
- KulczynskiCody() : Calculator("KulczynskiCody", 1) {};
+ KulczynskiCody() : Calculator("KulczynskiCody", 1, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
/***********************************************************************/
-EstOutput Lennon::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput Lennon::getValues(vector<SharedRAbundVector*> shared) {
try {
int S1, S2, S12, tempA, tempB, min;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0; min = 0;
data.resize(1,0);
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
if (tempA != 0) { S1++; }
if (tempB != 0) { S2++; }
class Lennon : public Calculator {
public:
- Lennon() : Calculator("Lennon", 1) {};
+ Lennon() : Calculator("Lennon", 1, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
#include "sharedmorisitahorn.h"
/***********************************************************************/
-EstOutput MorHorn::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput MorHorn::getValues(vector<SharedRAbundVector*> shared) {
try {
data.resize(1,0);
morhorn = 0.0; sumSharedA = 0.0; sumSharedB = 0.0; a = 0.0; b = 0.0; d = 0.0;
//get the total values we need to calculate the theta denominator sums
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- Atotal += shared1->getAbundance(i);
- Btotal += shared2->getAbundance(i);
+ Atotal += shared[0]->getAbundance(i);
+ Btotal += shared[1]->getAbundance(i);
}
//calculate the theta denominator sums
- for (int j = 0; j < shared1->size(); j++) {
+ for (int j = 0; j < shared[0]->size(); j++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(j);
- tempB = shared2->getAbundance(j);
+ tempA = shared[0]->getAbundance(j);
+ tempB = shared[1]->getAbundance(j);
a += tempA * tempA;
b += tempB * tempB;
class MorHorn : public Calculator {
public:
- MorHorn() : Calculator("MorisitaHorn", 1) {};
+ MorHorn() : Calculator("MorisitaHorn", 1, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
class SharedNSeqs : public Calculator {
public:
- SharedNSeqs() : Calculator("SharedNSeqs", 1) {};
+ SharedNSeqs() : Calculator("SharedNSeqs", 1, false) {};
EstOutput getValues(SAbundVector* rank){ return data; };
- EstOutput getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+ EstOutput getValues(vector<SharedRAbundVector*> shared) {
data.resize(1,0);
- data[0] = (double)shared1->getNumSeqs() + (double)shared2->getNumSeqs();
+ data[0] = (double)shared[0]->getNumSeqs() + (double)shared[1]->getNumSeqs();
return data;
}
};
/***********************************************************************/
-EstOutput Ochiai::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput Ochiai::getValues(vector<SharedRAbundVector*> shared) {
try {
int S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
data.resize(1,0);
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
if (tempA != 0) { S1++; }
if (tempB != 0) { S2++; }
class Ochiai : public Calculator {
public:
- Ochiai() : Calculator("Ochiai", 1) {};
+ Ochiai() : Calculator("Ochiai", 1, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
/***********************************************************************/
SharedOrderVector SharedOrderVector::getSharedOrderVector(){
+ random_shuffle(data.begin(), data.end());
return *this;
}
//This returns the number of unique species observed in several groups.
//The shared vector is each groups sharedrabundvector.
-EstOutput SharedSobs::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
+EstOutput SharedSobs::getValues(vector<SharedRAbundVector*> shared){
try {
data.resize(1,0);
int observed = 0;
//loop through the species in each group
- for (int k = 0; k < shared1->size(); k++) {
+ for (int k = 0; k < shared[0]->size(); k++) {
//if you have found a new species
- if (shared1->getAbundance(k) != 0) { observed++; }
- else if ((shared1->getAbundance(k) == 0) && (shared2->getAbundance(k) != 0)) { observed++; }
+ if (shared[0]->getAbundance(k) != 0) { observed++; }
+ else if ((shared[0]->getAbundance(k) == 0) && (shared[1]->getAbundance(k) != 0)) { observed++; }
}
data[0] = observed;
class SharedSobs : public Calculator {
public:
- SharedSobs() : Calculator("sharedsobs", 1) {};
+ SharedSobs() : Calculator("sharedsobs", 1, false) {};
EstOutput getValues(SAbundVector* rank){ return data; };
- EstOutput getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2);
+ EstOutput getValues(vector<SharedRAbundVector*>);
};
/***********************************************************************/
//This returns the number of shared species observed in several groups.
//The shared vector is each groups sharedrabundvector.
-EstOutput SharedSobsCS::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
+EstOutput SharedSobsCS::getValues(vector<SharedRAbundVector*> shared){
try {
data.resize(1,0);
int observed = 0;
- int tempA, tempB;
+ int numGroups = shared.size();
- //loop through the species in each group
- for (int k = 0; k < shared1->size(); k++) {
- tempA = shared1->getAbundance(k); //store in temps to avoid calling getAbundance multiple times
- tempB = shared2->getAbundance(k);
-
- //if you have found a new species
- if ((tempA != 0) && (tempB != 0)) {//they are shared
- observed++;
+ for (int i = 0; i < shared[0]->size(); i++) {
+ //get bin values and set sharedByAll
+ bool sharedByAll = true;
+ for (int j = 0; j < numGroups; j++) {
+ if (shared[j]->getAbundance(i) == 0) { sharedByAll = false; }
}
+
+ //they are shared
+ if (sharedByAll == true) { observed++; }
}
data[0] = observed;
class SharedSobsCS : public Calculator {
public:
- SharedSobsCS() : Calculator("sharedsobs", 1) {};
+ SharedSobsCS() : Calculator("sharedsobs", 1, true) {};
EstOutput getValues(SAbundVector* rank){ return data; };
- EstOutput getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2);
+ EstOutput getValues(vector<SharedRAbundVector*>);
};
/***********************************************************************/
/***********************************************************************/
-EstOutput SorAbund::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput SorAbund::getValues(vector<SharedRAbundVector*> shared) {
try {
EstOutput UVest;
UVest.resize(2,0);
data.resize(1,0);
- UVest = uv->getUVest(shared1, shared2);
+ UVest = uv->getUVest(shared);
//UVest[0] is Uest, UVest[1] is Vest
data[0] = (2 * UVest[0] * UVest[1]) / ((float)(UVest[0] + UVest[1]));
class SorAbund : public Calculator {
public:
- SorAbund() : Calculator("SorAbund", 3) {};
+ SorAbund() : Calculator("SorAbund", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
UVEst* uv;
/***********************************************************************/
-EstOutput SorClass::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput SorClass::getValues(vector<SharedRAbundVector*> shared) {
try {
int S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
data.resize(1,0);
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
//find number of bins in shared1 and shared2
if (tempA != 0) { S1++; }
class SorClass : public Calculator {
public:
- SorClass() : Calculator("SorClass", 3) {};
+ SorClass() : Calculator("SorClass", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
/***********************************************************************/
-EstOutput SorEst::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput SorEst::getValues(vector<SharedRAbundVector*> shared) {
try {
EstOutput S1, S2, S12;
S12.resize(1,0);
SAbundVector* chaoS1Sabund = new SAbundVector();
SAbundVector* chaoS2Sabund = new SAbundVector();
- *chaoS1Sabund = shared1->getSAbundVector();
- *chaoS2Sabund = shared2->getSAbundVector();
+ *chaoS1Sabund = shared[0]->getSAbundVector();
+ *chaoS2Sabund = shared[1]->getSAbundVector();
- S12 = sharedChao->getValues(shared1, shared2);
+ S12 = sharedChao->getValues(shared);
S1 = chaoS1->getValues(chaoS1Sabund);
S2 = chaoS2->getValues(chaoS2Sabund);
class SorEst : public Calculator {
public:
- SorEst() : Calculator("SorEst", 3) {};
+ SorEst() : Calculator("SorEst", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
#include "sharedthetan.h"
/***********************************************************************/
-EstOutput ThetaN::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput ThetaN::getValues(vector<SharedRAbundVector*> shared) {
try {
data.resize(1,0);
numerator = 0.0; denominator = 0.0; thetaN = 0.0; sumSharedA = 0.0; sumSharedB = 0.0; a = 0.0; b = 0.0; d = 0.0;
//get the total values we need to calculate the theta denominator sums
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- Atotal += shared1->getAbundance(i);
- Btotal += shared2->getAbundance(i);
+ Atotal += shared[0]->getAbundance(i);
+ Btotal += shared[1]->getAbundance(i);
}
//calculate the theta denominator sums
- for (int j = 0; j < shared1->size(); j++) {
+ for (int j = 0; j < shared[0]->size(); j++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(j);
- tempB = shared2->getAbundance(j);
+ tempA = shared[0]->getAbundance(j);
+ tempB = shared[1]->getAbundance(j);
//they are shared
if ((tempA != 0) && (tempB != 0)) {
class ThetaN : public Calculator {
public:
- ThetaN() : Calculator("ThetaN", 3) {};
+ ThetaN() : Calculator("ThetaN", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
#include "sharedthetayc.h"
/***********************************************************************/
-EstOutput ThetaYC::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput ThetaYC::getValues(vector<SharedRAbundVector*> shared) {
try {
data.resize(1,0);
float b = 0;
//get the total values we need to calculate the theta denominator sums
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- Atotal += shared1->getAbundance(i);
- Btotal += shared2->getAbundance(i);
+ Atotal += shared[0]->getAbundance(i);
+ Btotal += shared[1]->getAbundance(i);
}
//calculate the theta denominator sums
- for (int j = 0; j < shared1->size(); j++) {
+ for (int j = 0; j < shared[0]->size(); j++) {
//store in temps to avoid multiple repetitive function calls
- relA = shared1->getAbundance(j) / (float)Atotal;
- relB = shared2->getAbundance(j) / (float)Btotal;
+ relA = shared[0]->getAbundance(j) / (float)Atotal;
+ relB = shared[1]->getAbundance(j) / (float)Btotal;
a += relA * relB;
b += pow((relA-relB),2);
class ThetaYC : public Calculator {
public:
- ThetaYC() : Calculator("ThetaYC", 3) {};
+ ThetaYC() : Calculator("ThetaYC", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
private:
};
class Simpson : public Calculator {
public:
- Simpson() : Calculator("Simpson", 3) {};
+ Simpson() : Calculator("Simpson", 3, false) {};
EstOutput getValues(SAbundVector*);
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
};
/***********************************************************************/
class Sobs : public Calculator {
public:
- Sobs() : Calculator("sobs", 1) {};
+ Sobs() : Calculator("sobs", 1, false) {};
EstOutput getValues(SAbundVector* rank){
data.resize(1,0);
data[0] = (double)rank->getNumBins();
return data;
}
- EstOutput getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
};
/***********************************************************************/
format = globaldata->getFormat();
validCalculator = new ValidCalculators();
util = new SharedUtil();
+ mult = false;
int i;
for (i=0; i<globaldata->Estimators.size(); i++) {
}else if (globaldata->Estimators[i] == "whittaker") {
sumCalculators.push_back(new Whittaker());
}
-
}
}
//reset calc for next command
//if the users entered no valid calculators don't execute command
if (sumCalculators.size() == 0) { return 0; }
-
+ //check if any calcs can do multiples
+ else{
+ for (int i = 0; i < sumCalculators.size(); i++) {
+ if (sumCalculators[i]->getMultiple() == true) { mult = true; }
+ }
+ }
+
if (format == "sharedfile") {
read = new ReadOTUFile(globaldata->inputFileName);
read->read(&*globaldata);
}
outputFileHandle << endl;
+ //create file and put column headers for multiple groups file
+ if (mult = true) {
+ outAllFileName = ((getRootName(globaldata->inputFileName)) + "sharedmultiple.summary");
+ openOutputFile(outAllFileName, outAll);
+
+ outAll << "label" <<'\t' << "comparison" << '\t';
+ for(int i=0;i<sumCalculators.size();i++){
+ if (sumCalculators[i]->getMultiple() == true) {
+ outAll << '\t' << sumCalculators[i]->getName();
+ }
+ }
+ outAll << endl;
+ }
+
while(order != NULL){
if(globaldata->allLines == 1 || globaldata->lines.count(count) == 1 || globaldata->labels.count(order->getLabel()) == 1){
cout << order->getLabel() << '\t' << count << endl;
util->getSharedVectors(globaldata->Groups, lookup, order); //fills group vectors from order vector. //fills group vectors from order vector.
- //randomize group order
- if (globaldata->getJumble() == "1") { random_shuffle(lookup.begin(), lookup.end()); }
+ //loop through calculators and add to file all for all calcs that can do mutiple groups
+ if (mult = true) {
+ //output label
+ outAll << order->getLabel() << '\t';
+
+ //output groups names
+ string outNames = "";
+ for (int j = 0; j < lookup.size(); j++) {
+ outNames += lookup[j]->getGroup() + "-";
+ }
+ outNames = outNames.substr(0, outNames.length()-1); //rip off extra '-';
+ outAll << outNames << '\t';
+
+ for(int i=0;i<sumCalculators.size();i++){
+ if (sumCalculators[i]->getMultiple() == true) {
+ sumCalculators[i]->getValues(lookup);
+ outAll << '\t';
+ sumCalculators[i]->print(outAll);
+ }
+ }
+ outAll << endl;
+ }
int n = 1;
+ vector<SharedRAbundVector*> subset;
for (int k = 0; k < (lookup.size() - 1); k++) { // pass cdd each set of groups to commpare
for (int l = n; l < lookup.size(); l++) {
outputFileHandle << order->getLabel() << '\t';
+ subset.clear(); //clear out old pair of sharedrabunds
+ //add new pair of sharedrabunds
+ subset.push_back(lookup[k]); subset.push_back(lookup[l]);
+
//sort groups to be alphanumeric
if (lookup[k]->getGroup() > lookup[l]->getGroup()) {
outputFileHandle << (lookup[l]->getGroup() +'\t' + lookup[k]->getGroup()) << '\t'; //print out groups
outputFileHandle << (lookup[k]->getGroup() +'\t' + lookup[l]->getGroup()) << '\t'; //print out groups
}
- for(int i=0;i<sumCalculators.size();i++){
- sumCalculators[i]->getValues(lookup[k], lookup[l]); //saves the calculator outputs
+ for(int i=0;i<sumCalculators.size();i++) {
+ sumCalculators[i]->getValues(subset); //saves the calculator outputs
outputFileHandle << '\t';
sumCalculators[i]->print(outputFileHandle);
}
SharedListVector* SharedList;
SharedOrderVector* order;
vector<SharedRAbundVector*> lookup;
- string outputFileName, format;
- ofstream outputFileHandle;
+ string outputFileName, format, outAllFileName;
+ ofstream outputFileHandle, outAll;
+ bool mult;
};
try {
int count = 1;
EstOutput data;
-
+ vector<SharedRAbundVector*> subset;
+
//if the users entered no valid calculators don't execute command
if (treeCalculators.size() == 0) { return 0; }
for (int l = k; l < lookup.size(); l++) {
if (k != l) { //we dont need to similiarity of a groups to itself
//get estimated similarity between 2 groups
- data = treeCalculators[i]->getValues(lookup[k], lookup[l]); //saves the calculator outputs
+
+ subset.clear(); //clear out old pair of sharedrabunds
+ //add new pair of sharedrabunds
+ subset.push_back(lookup[k]); subset.push_back(lookup[l]);
+
+ data = treeCalculators[i]->getValues(subset); //saves the calculator outputs
//save values in similarity matrix
simMatrix[k][l] = data[0];
simMatrix[l][k] = data[0];
/***********************************************************************/
//This is used by SharedJAbund and SharedSorAbund
-EstOutput UVEst::getUVest(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
+EstOutput UVEst::getUVest(vector<SharedRAbundVector*> shared) {
try {
EstOutput results;
results.resize(2,0);
sumSharedA1 = the sum of all shared otus in A where B = 1
sumSharedB1 = the sum of all shared otus in B where A = 1 */
- for (int i = 0; i < shared1->size(); i++) {
+ for (int i = 0; i < shared[0]->size(); i++) {
//store in temps to avoid multiple repetitive function calls
- tempA = shared1->getAbundance(i);
- tempB = shared2->getAbundance(i);
+ tempA = shared[0]->getAbundance(i);
+ tempB = shared[1]->getAbundance(i);
Atotal += tempA;
Btotal += tempB;
/***********************************************************************/
class UVEst {
public:
- EstOutput getUVest(SharedRAbundVector* shared1, SharedRAbundVector* shared2);
+ EstOutput getUVest(vector<SharedRAbundVector*>);
};
/***********************************************************************/
commands["get.label"] = "get.label";
commands["get.line"] = "get.line";
commands["bootstrap.shared"] = "bootstrap.shared";
+ commands["concensus"] = "concensus";
commands["help"] = "help";
commands["quit"] = "quit";
string collectsingleArray[] = {"freq","line","label","calc","abund"};
commandParameters["collect.single"] = addParameters(collectsingleArray, sizeof(collectsingleArray)/sizeof(string));
- string collectsharedArray[] = {"jumble","freq","line","label","calc","groups"};
+ string collectsharedArray[] = {"freq","line","label","calc","groups"};
commandParameters["collect.shared"] = addParameters(collectsharedArray, sizeof(collectsharedArray)/sizeof(string));
string getgroupArray[] = {};
string summarysingleArray[] = {"line","label","calc","abund"};
commandParameters["summary.single"] = addParameters(summarysingleArray, sizeof(summarysingleArray)/sizeof(string));
- string summarysharedArray[] = {"jumble","line","label","calc","groups"};
+ string summarysharedArray[] = {"line","label","calc","groups"};
commandParameters["summary.shared"] = addParameters(summarysharedArray, sizeof(summarysharedArray)/sizeof(string));
string parsimonyArray[] = {"random","groups","iters"};
string bootstrapArray[] = {"line","label","calc","groups","iters"};
commandParameters["bootstrap.shared"] = addParameters(bootstrapArray, sizeof(bootstrapArray)/sizeof(string));
+ string concensusArray[] = {};
+ commandParameters["concensus"] = addParameters(concensusArray, sizeof(concensusArray)/sizeof(string));
+
string quitArray[] = {};
commandParameters["quit"] = addParameters(quitArray, sizeof(quitArray)/sizeof(string));
void Venn::getPic(vector<SharedRAbundVector*> lookup, vector<Calculator*> vCalcs) {
try {
- //fills vector of sharedsabunds - lookup
- //util->getSharedVectors(globaldata->Groups, lookup, sharedorder); //fills group vectors from order vector.
+ vector<SharedRAbundVector*> subset;
/******************* 1 Group **************************/
if (lookup.size() == 1) {
sabundA = new SAbundVector(lookup[0]->getSAbundVector());// sabundA = &sA;
sabundB = new SAbundVector(lookup[1]->getSAbundVector());// sabundB = &sB;
+ subset.clear();
+ subset.push_back(lookup[0]); subset.push_back(lookup[1]);
+
//make a file for each calculator
for(int i=0;i<vCalcs.size();i++){
string filenamesvg = getRootName(globaldata->inputFileName) + lookup[0]->getLabel() + ".venn." + vCalcs[i]->getName() + ".svg";
openOutputFile(filenamesvg, outsvg);
//get estimates for sharedAB
- vector<double> shared = vCalcs[i]->getValues(lookup[0], lookup[1]);
+ vector<double> shared = vCalcs[i]->getValues(subset);
//in essence you want to run it like a single
if (vCalcs[i]->getName() == "sharedsobs") {
for(int i=0;i<vCalcs.size();i++){
string filenamesvg = getRootName(globaldata->inputFileName) + lookup[0]->getLabel() + ".venn." + vCalcs[i]->getName() + ".svg";
openOutputFile(filenamesvg, outsvg);
-
+
//get estimates for sharedAB, sharedAC and sharedBC
- vector<double> sharedAB = vCalcs[i]->getValues(lookup[0], lookup[1]);
- vector<double> sharedAC = vCalcs[i]->getValues(lookup[0], lookup[2]);
- vector<double> sharedBC = vCalcs[i]->getValues(lookup[1], lookup[2]);
+ subset.clear();
+ subset.push_back(lookup[0]); subset.push_back(lookup[1]);
+ vector<double> sharedAB = vCalcs[i]->getValues(subset);
+
+ subset.clear();
+ subset.push_back(lookup[0]); subset.push_back(lookup[2]);
+ vector<double> sharedAC = vCalcs[i]->getValues(subset);
+
+ subset.clear();
+ subset.push_back(lookup[1]); subset.push_back(lookup[2]);
+ vector<double> sharedBC = vCalcs[i]->getValues(subset);
//merge BC and estimate with shared with A
SharedRAbundVector* merge = new SharedRAbundVector();
merge->push_back((lookup[1]->getAbundance(j) + lookup[2]->getAbundance(j)), j, "");
}
- vector<double> sharedAwithBC = vCalcs[i]->getValues(lookup[0], merge);
+ subset.clear();
+ subset.push_back(lookup[0]); subset.push_back(merge);
+ vector<double> sharedAwithBC = vCalcs[i]->getValues(subset);
delete merge;
//merge AC and estimate with shared with B
for (int j = 0; j < lookup[0]->size(); j++) {
merge->push_back((lookup[0]->getAbundance(j) + lookup[2]->getAbundance(j)), j, "");
}
-
- vector<double> sharedBwithAC = vCalcs[i]->getValues(lookup[1], merge);
+
+ subset.clear();
+ subset.push_back(merge); subset.push_back(lookup[1]);
+ vector<double> sharedBwithAC = vCalcs[i]->getValues(subset);
delete merge;
//merge AB and estimate with shared with C
merge->push_back((lookup[0]->getAbundance(j) + lookup[1]->getAbundance(j)), j, "");
}
- vector<double> sharedCwithAB = vCalcs[i]->getValues(lookup[2], merge);
+ subset.clear();
+ subset.push_back(lookup[2]); subset.push_back(merge);
+ vector<double> sharedCwithAB = vCalcs[i]->getValues(subset);
delete merge;
//in essence you want to run it like a single
/***********************************************************************/
-EstOutput Whittaker::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
+EstOutput Whittaker::getValues(vector<SharedRAbundVector*> shared){
try{
data.resize(1);
int countA = 0;
int countB = 0;
- int sTotal = shared1->getNumBins();
+ int sTotal = shared[0]->getNumBins();
for(int i=0;i<sTotal;i++){
- if(shared1->getAbundance(i) != 0){ countA++; }
- if(shared2->getAbundance(i) != 0){ countB++; }
+ if(shared[0]->getAbundance(i) != 0){ countA++; }
+ if(shared[1]->getAbundance(i) != 0){ countB++; }
}
data[0] = 2*sTotal/(float)(countA+countB)-1;
class Whittaker : public Calculator {
public:
- Whittaker() : Calculator("whittaker", 3) {};
+ Whittaker() : Calculator("whittaker", 3, false) {};
EstOutput getValues(SAbundVector*) {return data;};
- EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>);
};
projects / mothur.git / commitdiff