*/
#include "phylodiversitycommand.h"
-#include "phylodiversity.h"
//**********************************************************************************************************************
PhyloDiversityCommand::PhyloDiversityCommand(string option) {
else {
//valid paramters for this command
- string Array[] = {"freq","rarefy","iters","groups","outputdir","inputdir"};
+ string Array[] = {"freq","rarefy","iters","groups","summary","collect","scale","outputdir","inputdir"};
vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
OptionParser parser(option);
}
//if the user changes the output directory command factory will send this info to us in the output parameter
- outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
+ outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = hasPath(globaldata->getTreeFile()); }
if (globaldata->gTree.size() == 0) {//no trees were read
m->mothurOut("You must execute the read.tree command, before you may execute the phylo.diversity command."); m->mothurOutEndLine(); abort = true; }
rarefy = isTrue(temp);
if (!rarefy) { iters = 1; }
+ temp = validParameter.validFile(parameters, "summary", false); if (temp == "not found") { temp = "T"; }
+ summary = isTrue(temp);
+
+ temp = validParameter.validFile(parameters, "scale", false); if (temp == "not found") { temp = "F"; }
+ scale = isTrue(temp);
+
+ temp = validParameter.validFile(parameters, "collect", false); if (temp == "not found") { temp = "F"; }
+ collect = isTrue(temp);
+
groups = validParameter.validFile(parameters, "groups", false);
if (groups == "not found") { groups = ""; Groups = globaldata->gTreemap->namesOfGroups; globaldata->Groups = Groups; }
else {
splitAtDash(groups, Groups);
globaldata->Groups = Groups;
}
+
+ if ((!collect) && (!rarefy) && (!summary)) { m->mothurOut("No outputs selected. You must set either collect, rarefy or summary to true, summary=T by default."); m->mothurOutEndLine(); abort=true; }
}
}
void PhyloDiversityCommand::help(){
try {
-
+ m->mothurOut("The phylo.diversity command can only be executed after a successful read.tree command.\n");
+ m->mothurOut("The phylo.diversity command parameters are groups, iters, freq, scale, rarefy, collect and summary. No parameters are required.\n");
+ m->mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed. The group names are separated by dashes. By default all groups are used.\n");
+ m->mothurOut("The iters parameter allows you to specify the number of randomizations to preform, by default iters=1000, if you set rarefy to true.\n");
+ m->mothurOut("The freq parameter is used indicate when to output your data, by default it is set to 100. But you can set it to a percentage of the number of sequence. For example freq=0.10, means 10%. \n");
+ m->mothurOut("The scale parameter is used indicate that you want your ouptut scaled to the number of sequences sampled, default = false. \n");
+ m->mothurOut("The rarefy parameter allows you to create a rarefaction curve. The default is false.\n");
+ m->mothurOut("The collect parameter allows you to create a collectors curve. The default is false.\n");
+ m->mothurOut("The summary parameter allows you to create a .summary file. The default is true.\n");
+ m->mothurOut("The phylo.diversity command should be in the following format: phylo.diversity(groups=yourGroups, rarefy=yourRarefy, iters=yourIters).\n");
+ m->mothurOut("Example phylo.diversity(groups=A-B-C, rarefy=T, iters=500).\n");
+ m->mothurOut("The phylo.diversity command output two files: .phylo.diversity and if rarefy=T, .rarefaction.\n");
+ m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
}
catch(exception& e) {
for (int i = 0; i < globaldata->Groups.size(); i++) { if (globaldata->Groups[i] == "xxx") { globaldata->Groups.erase(globaldata->Groups.begin()+i); break; } }
vector<string> outputNames;
-
- //diversity calculator
- PhyloDiversity phylo(globaldata->gTreemap);
-
+
vector<Tree*> trees = globaldata->gTree;
//for each of the users trees
if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
- phylo.setTotalGroupBranchLengths(trees[i]);
+ ofstream outSum, outRare, outCollect;
+ string outSumFile = outputDir + getRootName(getSimpleName(globaldata->getTreeFile())) + toString(i+1) + ".phylodiv.summary";
+ string outRareFile = outputDir + getRootName(getSimpleName(globaldata->getTreeFile())) + toString(i+1) + ".phylodiv.rarefaction";
+ string outCollectFile = outputDir + getRootName(getSimpleName(globaldata->getTreeFile())) + toString(i+1) + ".phylodiv";
- string outFile = outputDir + getRootName(getSimpleName(globaldata->getTreeFile())) + toString(i+1) + ".phylo.diversity";
- if (rarefy) { outFile += ".rarefaction"; }
- outputNames.push_back(outFile);
+ if (summary) { openOutputFile(outSumFile, outSum); outputNames.push_back(outSumFile); }
+ if (rarefy) { openOutputFile(outRareFile, outRare); outputNames.push_back(outRareFile); }
+ if (collect) { openOutputFile(outCollectFile, outCollect); outputNames.push_back(outCollectFile); }
int numLeafNodes = trees[i]->getNumLeaves();
numLeafNodes = randomLeaf.size(); //reset the number of leaf nodes you are using
//each group, each sampling, if no rarefy iters = 1;
- vector< vector<float> > diversity;
- diversity.resize(globaldata->Groups.size());
+ map<string, vector<float> > diversity;
+
+ //each group, each sampling, if no rarefy iters = 1;
+ map<string, vector<float> > sumDiversity;
+
+ //find largest group total
+ int largestGroup = 0;
+ for (int j = 0; j < globaldata->Groups.size(); j++) {
+ if (globaldata->gTreemap->seqsPerGroup[globaldata->Groups[j]] > largestGroup) { largestGroup = globaldata->gTreemap->seqsPerGroup[globaldata->Groups[j]]; }
+
+ //initialize diversity
+ diversity[globaldata->Groups[j]].resize(globaldata->gTreemap->seqsPerGroup[globaldata->Groups[j]]+1, 0.0); //numSampled
+ //groupA 0.0 0.0
+
+ //initialize sumDiversity
+ sumDiversity[globaldata->Groups[j]].resize(globaldata->gTreemap->seqsPerGroup[globaldata->Groups[j]]+1, 0.0);
+ }
+
+ //convert freq percentage to number
+ int increment = 100;
+ if (freq < 1.0) { increment = largestGroup * freq;
+ }else { increment = freq; }
//initialize sampling spots
- vector<int> numSampledList;
- for(int k = 0; k < numLeafNodes; k++){ if((k == 0) || (k+1) % freq == 0){ numSampledList.push_back(k); } }
- if(numLeafNodes % freq != 0){ numSampledList.push_back(numLeafNodes); }
-
- //initialize diversity
- for (int j = 0; j < diversity.size(); j++) { diversity[j].resize(numSampledList.size(), 0.0); } // 10sampled 20 sampled ...
- //groupA 0.0 0.0
- //then for each iter you add to score and then when printing divide by iters to get average
+ set<int> numSampledList;
+ for(int k = 1; k <= largestGroup; k++){ if((k == 1) || (k % increment == 0)){ numSampledList.insert(k); } }
+ if(largestGroup % increment != 0){ numSampledList.insert(largestGroup); }
+
+ //add other groups ending points
+ for (int j = 0; j < globaldata->Groups.size(); j++) {
+ if (numSampledList.count(diversity[globaldata->Groups[j]].size()-1) == 0) { numSampledList.insert(diversity[globaldata->Groups[j]].size()-1); }
+ }
+
for (int l = 0; l < iters; l++) {
random_shuffle(randomLeaf.begin(), randomLeaf.end());
- vector<int> leavesSampled;
- EstOutput data;
- int count = 0;
+ //initialize counts
+ map<string, int> counts;
+ for (int j = 0; j < globaldata->Groups.size(); j++) { counts[globaldata->Groups[j]] = 0; }
+
for(int k = 0; k < numLeafNodes; k++){
if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
- leavesSampled.push_back(randomLeaf[k]);
-
- if((k == 0) || (k+1) % freq == 0){ //ready to calc?
-
- data = phylo.getValues(trees[i], leavesSampled);
-
- //datas results are in the same order as globaldatas groups
- for (int h = 0; h < data.size(); h++) { diversity[h][count] += data[h]; }
-
- count++;
+ //calc branch length of randomLeaf k
+ float br = calcBranchLength(trees[i], randomLeaf[k]);
+
+ //for each group in the groups update the total branch length accounting for the names file
+ vector<string> groups = trees[i]->tree[randomLeaf[k]].getGroup();
+ for (int j = 0; j < groups.size(); j++) {
+ int numSeqsInGroupJ = 0;
+ map<string, int>::iterator it;
+ it = trees[i]->tree[randomLeaf[k]].pcount.find(groups[j]);
+ if (it != trees[i]->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
+ numSeqsInGroupJ = it->second;
+ }
+
+ for (int s = (counts[groups[j]]+1); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
+ diversity[groups[j]][s] = diversity[groups[j]][s-1] + ((float) numSeqsInGroupJ * br);
+ }
+ counts[groups[j]] += numSeqsInGroupJ;
}
}
-
- if(numLeafNodes % freq != 0){
-
- data = phylo.getValues(trees[i], leavesSampled);
-
- //datas results are in the same order as globaldatas groups
- for (int h = 0; h < data.size(); h++) { diversity[h][count] += data[h]; }
+
+ if (rarefy) {
+ //add this diversity to the sum
+ for (int j = 0; j < globaldata->Groups.size(); j++) {
+ for (int g = 0; g < diversity[globaldata->Groups[j]].size(); g++) {
+ sumDiversity[globaldata->Groups[j]][g] += diversity[globaldata->Groups[j]][g];
+ }
+ }
}
+
+ if ((collect) && (l == 0)) { printData(numSampledList, diversity, outCollect, 1); }
+ if ((summary) && (l == 0)) { printSumData(diversity, outSum, 1); }
}
- printData(numSampledList, diversity, outFile);
-
+ if (rarefy) { printData(numSampledList, sumDiversity, outRare, iters); }
}
}
//**********************************************************************************************************************
-void PhyloDiversityCommand::printData(vector<int>& num, vector< vector<float> >& div, string file){
+void PhyloDiversityCommand::printSumData(map< string, vector<float> >& div, ofstream& out, int numIters){
try {
- ofstream out;
- openOutputFile(file, out);
out << "numSampled\t";
for (int i = 0; i < globaldata->Groups.size(); i++) { out << globaldata->Groups[i] << '\t'; }
out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
- for (int i = 0; i < num.size(); i++) {
- if (i == (num.size()-1)) { out << num[i] << '\t'; }
- else { out << (num[i]+1) << '\t'; }
+ set<int> num;
+ //find end points to output
+ for (map<string, vector<float> >::iterator itEnds = div.begin(); itEnds != div.end(); itEnds++) { num.insert(itEnds->second.size()-1); }
+
+ for (set<int>::iterator it = num.begin(); it != num.end(); it++) {
+ int numSampled = *it;
+
+ out << numSampled << '\t';
- for (int j = 0; j < div.size(); j++) {
- float score = div[j][i] / (float)iters;
- out << setprecision(6) << score << '\t';
+ for (int j = 0; j < globaldata->Groups.size(); j++) {
+ if (numSampled < div[globaldata->Groups[j]].size()) {
+ float score;
+ if (scale) { score = (div[globaldata->Groups[j]][numSampled] / (float)numIters) / (float)numSampled; }
+ else { score = div[globaldata->Groups[j]][numSampled] / (float)numIters; }
+
+ out << setprecision(4) << score << '\t';
+ }else { out << "NA" << '\t'; }
}
out << endl;
}
}
catch(exception& e) {
- m->errorOut(e, "PhyloDiversityCommand", "printData");
+ m->errorOut(e, "PhyloDiversityCommand", "printSumData");
exit(1);
}
}
+//**********************************************************************************************************************
+
+void PhyloDiversityCommand::printData(set<int>& num, map< string, vector<float> >& div, ofstream& out, int numIters){
+ try {
+
+ out << "numSampled\t";
+ for (int i = 0; i < globaldata->Groups.size(); i++) { out << globaldata->Groups[i] << '\t'; }
+ out << endl;
+
+ out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
+
+ for (set<int>::iterator it = num.begin(); it != num.end(); it++) {
+ int numSampled = *it;
+
+ out << numSampled << '\t';
+
+ for (int j = 0; j < globaldata->Groups.size(); j++) {
+ if (numSampled < div[globaldata->Groups[j]].size()) {
+ float score;
+ if (scale) { score = (div[globaldata->Groups[j]][numSampled] / (float)numIters) / (float)numSampled; }
+ else { score = div[globaldata->Groups[j]][numSampled] / (float)numIters; }
+ out << setprecision(4) << score << '\t';
+ }else { out << "NA" << '\t'; }
+ }
+ out << endl;
+ }
+
+ out.close();
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "PhyloDiversityCommand", "printData");
+ exit(1);
+ }
+}
//**********************************************************************************************************************
+float PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf){
+ try {
+
+ //calc the branch length
+ //while you aren't at root
+ float sum = 0.0;
+ int index = leaf;
+
+ while(t->tree[index].getParent() != -1){
+
+ //if you have a BL
+ if(t->tree[index].getBranchLength() != -1){
+ sum += abs(t->tree[index].getBranchLength());
+ }
+ index = t->tree[index].getParent();
+ }
+
+ //get last breanch length added
+ if(t->tree[index].getBranchLength() != -1){
+ sum += abs(t->tree[index].getBranchLength());
+ }
+
+ return sum;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "PhyloDiversityCommand", "calcBranchLength");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
\ No newline at end of file
projects / mothur.git / blobdiff