*/
#include "phylosummary.h"
-
+#include "referencedb.h"
/**************************************************************************************************/
-PhyloSummary::PhyloSummary(string refTfile, string groupFile){
+PhyloSummary::PhyloSummary(string refTfile, CountTable* c){
try {
m = MothurOut::getInstance();
maxLevel = 0;
ignore = false;
+ numSeqs = 0;
+
+ ct = c;
+ groupmap = NULL;
+
+ //check for necessary files
+ if (refTfile == "saved") { ReferenceDB* rdb = ReferenceDB::getInstance(); refTfile = rdb->getSavedTaxonomy(); }
+ string taxFileNameTest = m->getFullPathName((refTfile.substr(0,refTfile.find_last_of(".")+1) + "tree.sum"));
+ ifstream FileTest(taxFileNameTest.c_str());
- if (groupFile != "") {
- groupmap = new GroupMap(groupFile);
- groupmap->readMap();
+ if (!FileTest) {
+ m->mothurOut("Error: can't find " + taxFileNameTest + "."); m->mothurOutEndLine(); exit(1);
}else{
- groupmap = NULL;
+ readTreeStruct(FileTest);
}
+
+ tree[0].rank = "0";
+ assignRank(0);
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "PhyloSummary", "PhyloSummary");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+PhyloSummary::PhyloSummary(CountTable* c){
+ try {
+ m = MothurOut::getInstance();
+ maxLevel = 0;
+ ignore = true;
+ numSeqs = 0;
+
+ ct = c;
+ groupmap = NULL;
+
+ tree.push_back(rawTaxNode("Root"));
+ tree[0].rank = "0";
+ }
+ catch(exception& e) {
+ m->errorOut(e, "PhyloSummary", "PhyloSummary");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+PhyloSummary::PhyloSummary(string refTfile, GroupMap* g){
+ try {
+ m = MothurOut::getInstance();
+ maxLevel = 0;
+ ignore = false;
+ numSeqs = 0;
+
+ groupmap = g;
+ ct = NULL;
//check for necessary files
+ if (refTfile == "saved") { ReferenceDB* rdb = ReferenceDB::getInstance(); refTfile = rdb->getSavedTaxonomy(); }
string taxFileNameTest = m->getFullPathName((refTfile.substr(0,refTfile.find_last_of(".")+1) + "tree.sum"));
ifstream FileTest(taxFileNameTest.c_str());
/**************************************************************************************************/
-PhyloSummary::PhyloSummary(string groupFile){
+PhyloSummary::PhyloSummary(GroupMap* g){
try {
m = MothurOut::getInstance();
maxLevel = 0;
ignore = true;
+ numSeqs = 0;
- if (groupFile != "") {
- groupmap = new GroupMap(groupFile);
- groupmap->readMap();
- }else{
- groupmap = NULL;
- }
+ groupmap = g;
+ ct = NULL;
tree.push_back(rawTaxNode("Root"));
tree[0].rank = "0";
-
-
}
catch(exception& e) {
m->errorOut(e, "PhyloSummary", "PhyloSummary");
int PhyloSummary::summarize(string userTfile){
try {
-
- ifstream in;
- m->openInputFile(userTfile, in);
-
- //read in users taxonomy file and add sequences to tree
- string name, tax;
- int numSeqs = 0;
- while(!in.eof()){
- in >> name >> tax; m->gobble(in);
-
- addSeqToTree(name, tax);
- numSeqs++;
-
- if (m->control_pressed) { break; }
- }
- in.close();
-
- return numSeqs;
+ map<string, string> temp;
+ m->readTax(userTfile, temp);
+
+ for (map<string, string>::iterator itTemp = temp.begin(); itTemp != temp.end();) {
+ addSeqToTree(itTemp->first, itTemp->second);
+ temp.erase(itTemp++);
+ }
+
+ return numSeqs;
}
catch(exception& e) {
m->errorOut(e, "PhyloSummary", "summarize");
childPointer = tree[currentNode].children.find(taxon);
if(childPointer != tree[currentNode].children.end()){ //if the node already exists, update count and move on
- if (groupmap != NULL) {
+ int thisCount = 1;
+
+ if (groupmap != NULL) {
//find out the sequences group
string group = groupmap->getGroup(seqName);
if (itGroup != tree[childPointer->second].groupCount.end()) {
tree[childPointer->second].groupCount[group]++;
}
- }
+ }else if (ct != NULL) {
+ if (ct->hasGroupInfo()) {
+ vector<int> groupCounts = ct->getGroupCounts(seqName);
+ vector<string> groups = ct->getNamesOfGroups();
+ for (int i = 0; i < groups.size(); i++) {
+
+ if (groupCounts[i] != 0) {
+ //do you have a count for this group?
+ map<string, int>::iterator itGroup = tree[childPointer->second].groupCount.find(groups[i]);
+
+ //if yes, increment it - there should not be a case where we can't find it since we load group in read
+ if (itGroup != tree[childPointer->second].groupCount.end()) {
+ tree[childPointer->second].groupCount[groups[i]] += groupCounts[i];
+ }
+ }
+ }
+ }
+ thisCount = ct->getNumSeqs(seqName);
+ }
- tree[childPointer->second].total++;
+ tree[childPointer->second].total += thisCount;
currentNode = childPointer->second;
}else{
tree[index].parent = currentNode;
tree[index].level = (level+1);
- tree[index].total = 1;
tree[currentNode].children[taxon] = index;
+ int thisCount = 1;
//initialize groupcounts
if (groupmap != NULL) {
//if yes, increment it - there should not be a case where we can't find it since we load group in read
if (itGroup != tree[index].groupCount.end()) {
tree[index].groupCount[group]++;
- }
- }
+ }
+ }else if (ct != NULL) {
+ if (ct->hasGroupInfo()) {
+ vector<string> mGroups = ct->getNamesOfGroups();
+ for (int j = 0; j < mGroups.size(); j++) {
+ tree[index].groupCount[mGroups[j]] = 0;
+ }
+ vector<int> groupCounts = ct->getGroupCounts(seqName);
+ vector<string> groups = ct->getNamesOfGroups();
+
+ for (int i = 0; i < groups.size(); i++) {
+ if (groupCounts[i] != 0) {
+
+ //do you have a count for this group?
+ map<string, int>::iterator itGroup = tree[index].groupCount.find(groups[i]);
+
+ //if yes, increment it - there should not be a case where we can't find it since we load group in read
+ if (itGroup != tree[index].groupCount.end()) {
+ tree[index].groupCount[groups[i]]+=groupCounts[i];
+ }
+ }
+ }
+ }
+ thisCount = ct->getNumSeqs(seqName);
+ }
+ tree[index].total = thisCount;
currentNode = index;
}else{ //otherwise, error
}
/**************************************************************************************************/
-int PhyloSummary::addSeqToTree(string seqTaxonomy, vector<string> names){
+int PhyloSummary::addSeqToTree(string seqTaxonomy, map<string, bool> containsGroup){
try {
numSeqs++;
childPointer = tree[currentNode].children.find(taxon);
if(childPointer != tree[currentNode].children.end()){ //if the node already exists, update count and move on
- if (groupmap != NULL) {
-
- map<string, bool> containsGroup;
- vector<string> mGroups = groupmap->getNamesOfGroups();
- for (int j = 0; j < mGroups.size(); j++) {
- containsGroup[mGroups[j]] = false;
- }
-
- for (int k = 0; k < names.size(); k++) {
- //find out the sequences group
- string group = groupmap->getGroup(names[k]);
+ for (map<string, bool>::iterator itGroup = containsGroup.begin(); itGroup != containsGroup.end(); itGroup++) {
+ if (itGroup->second == true) {
+ tree[childPointer->second].groupCount[itGroup->first]++;
+ }
+ }
- if (group == "not found") { m->mothurOut("[WARNING]: " + names[k] + " is not in your groupfile, and will be included in the overall total, but not any group total."); m->mothurOutEndLine(); }
- else {
- containsGroup[group] = true;
- }
- }
-
- for (map<string, bool>::iterator itGroup = containsGroup.begin(); itGroup != containsGroup.end(); itGroup++) {
- if (itGroup->second == true) {
- tree[childPointer->second].groupCount[itGroup->first]++;
- }
- }
-
- }
-
tree[childPointer->second].total++;
currentNode = childPointer->second;
tree[index].level = (level+1);
tree[index].total = 1;
tree[currentNode].children[taxon] = index;
-
- //initialize groupcounts
- if (groupmap != NULL) {
- map<string, bool> containsGroup;
- vector<string> mGroups = groupmap->getNamesOfGroups();
- for (int j = 0; j < mGroups.size(); j++) {
- tree[index].groupCount[mGroups[j]] = 0;
- containsGroup[mGroups[j]] = false;
- }
-
- for (int k = 0; k < names.size(); k++) {
- //find out the sequences group
- string group = groupmap->getGroup(names[k]);
-
- if (group == "not found") { m->mothurOut("[WARNING]: " + names[k] + " is not in your groupfile, and will be included in the overall total, but not any group total."); m->mothurOutEndLine(); }
- else {
- containsGroup[group] = true;
- }
- }
-
- for (map<string, bool>::iterator itGroup = containsGroup.begin(); itGroup != containsGroup.end(); itGroup++) {
- if (itGroup->second == true) {
- tree[index].groupCount[itGroup->first]++;
- }
- }
- }
+ for (map<string, bool>::iterator itGroup = containsGroup.begin(); itGroup != containsGroup.end(); itGroup++) {
+ if (itGroup->second == true) {
+ tree[index].groupCount[itGroup->first]++;
+ }
+ }
currentNode = index;
try {
if (ignore) { assignRank(0); }
-
+ vector<string> mGroups;
//print labels
out << "taxlevel\t rankID\t taxon\t daughterlevels\t total\t";
if (groupmap != NULL) {
//so the labels match the counts below, since the map sorts them automatically...
//sort(groupmap->namesOfGroups.begin(), groupmap->namesOfGroups.end());
- vector<string> mGroups = groupmap->getNamesOfGroups();
+ mGroups = groupmap->getNamesOfGroups();
for (int i = 0; i < mGroups.size(); i++) {
out << mGroups[i] << '\t';
}
- }
+ }else if (ct != NULL) {
+ if (ct->hasGroupInfo()) {
+ mGroups = ct->getNamesOfGroups();
+ for (int i = 0; i < mGroups.size(); i++) {
+ out << mGroups[i] << '\t';
+ }
+ }
+ }
out << endl;
tree[0].total += tree[it->second].total;
if (groupmap != NULL) {
- vector<string> mGroups = groupmap->getNamesOfGroups();
for (int i = 0; i < mGroups.size(); i++) { tree[0].groupCount[mGroups[i]] += tree[it->second].groupCount[mGroups[i]]; }
- }
+ }else if ( ct != NULL) {
+ if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) { tree[0].groupCount[mGroups[i]] += tree[it->second].groupCount[mGroups[i]]; } }
+ }
}
}
if (groupmap != NULL) {
- //for (itGroup = tree[0].groupCount.begin(); itGroup != tree[0].groupCount.end(); itGroup++) {
- // out << itGroup->second << '\t';
- //}
- vector<string> mGroups = groupmap->getNamesOfGroups();
- for (int i = 0; i < mGroups.size(); i++) { out << tree[0].groupCount[mGroups[i]] << '\t'; }
- }
+ for (int i = 0; i < mGroups.size(); i++) { out << tree[0].groupCount[mGroups[i]] << '\t'; }
+ }else if ( ct != NULL) {
+ if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) { out << tree[0].groupCount[mGroups[i]] << '\t'; } }
+ }
out << endl;
//print rest
exit(1);
}
}
+/**************************************************************************************************/
+void PhyloSummary::print(ofstream& out, bool relabund){
+ try {
+
+ if (ignore) { assignRank(0); }
+
+ int totalChildrenInTree = 0;
+ map<string, int>::iterator itGroup;
+
+ map<string,int>::iterator it;
+ for(it=tree[0].children.begin();it!=tree[0].children.end();it++){
+ if (tree[it->second].total != 0) {
+ totalChildrenInTree++;
+ tree[0].total += tree[it->second].total;
+
+ if (groupmap != NULL) {
+ vector<string> mGroups = groupmap->getNamesOfGroups();
+ for (int i = 0; i < mGroups.size(); i++) { tree[0].groupCount[mGroups[i]] += tree[it->second].groupCount[mGroups[i]]; }
+ }else if ( ct != NULL) {
+ vector<string> mGroups = ct->getNamesOfGroups();
+ if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) { tree[0].groupCount[mGroups[i]] += tree[it->second].groupCount[mGroups[i]]; } }
+ }
+ }
+ }
+
+ //print root
+ out << tree[0].name << "\t" << "1.0000" << "\t"; //root relative abundance is 1, everyone classifies to root
+
+ /*
+ if (groupmap != NULL) {
+ for (int i = 0; i < mGroups.size(); i++) { out << tree[0].groupCount[mGroups[i]] << '\t'; }
+ }else if ( ct != NULL) {
+ if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) { out << tree[0].groupCount[mGroups[i]] << '\t'; } }
+ }*/
+
+ if (groupmap != NULL) {
+ vector<string> mGroups = groupmap->getNamesOfGroups();
+ for (int i = 0; i < mGroups.size(); i++) { out << "1.0000" << '\t'; }
+ }else if ( ct != NULL) {
+ vector<string> mGroups = ct->getNamesOfGroups();
+ if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) { out << "1.0000" << '\t'; } }
+ }
+
+ out << endl;
+
+ //print rest
+ print(0, out, relabund);
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "PhyloSummary", "print");
+ exit(1);
+ }
+}
/**************************************************************************************************/
void PhyloSummary::print(int i, ofstream& out){
for(it=tree[i].children.begin();it!=tree[i].children.end();it++){
if (tree[it->second].total != 0) {
-
+
int totalChildrenInTree = 0;
-
+
map<string,int>::iterator it2;
- for(it2=tree[it->second].children.begin();it2!=tree[it->second].children.end();it2++){
+ for(it2=tree[it->second].children.begin();it2!=tree[it->second].children.end();it2++){
if (tree[it2->second].total != 0) { totalChildrenInTree++; }
}
-
+
out << tree[it->second].level << "\t" << tree[it->second].rank << "\t" << tree[it->second].name << "\t" << totalChildrenInTree << "\t" << tree[it->second].total << "\t";
map<string, int>::iterator itGroup;
// out << itGroup->second << '\t';
//}
vector<string> mGroups = groupmap->getNamesOfGroups();
- for (int i = 0; i < mGroups.size(); i++) { out << tree[it->second].groupCount[mGroups[i]] << '\t'; }
- }
+ for (int i = 0; i < mGroups.size(); i++) { out << tree[it->second].groupCount[mGroups[i]] << '\t'; }
+ }else if (ct != NULL) {
+ if (ct->hasGroupInfo()) {
+ vector<string> mGroups = ct->getNamesOfGroups();
+ for (int i = 0; i < mGroups.size(); i++) { out << tree[it->second].groupCount[mGroups[i]] << '\t'; }
+ }
+ }
out << endl;
}
exit(1);
}
}
+
+/**************************************************************************************************/
+
+void PhyloSummary::print(int i, ofstream& out, bool relabund){
+ try {
+ map<string,int>::iterator it;
+ for(it=tree[i].children.begin();it!=tree[i].children.end();it++){
+
+ if (tree[it->second].total != 0) {
+
+ int totalChildrenInTree = 0;
+
+ map<string,int>::iterator it2;
+ for(it2=tree[it->second].children.begin();it2!=tree[it->second].children.end();it2++){
+ if (tree[it2->second].total != 0) { totalChildrenInTree++; }
+ }
+
+ string nodeName = "";
+ int thisNode = it->second;
+ while (tree[thisNode].rank != "0") { //while you are not at top
+ if (m->control_pressed) { break; }
+ nodeName = tree[thisNode].name + "|" + nodeName;
+ thisNode = tree[thisNode].parent;
+ }
+ if (nodeName != "") { nodeName = nodeName.substr(0, nodeName.length()-1); }
+
+ out << nodeName << "\t" << (tree[it->second].total / (float)tree[i].total) << "\t";
+
+ map<string, int>::iterator itGroup;
+ if (groupmap != NULL) {
+ vector<string> mGroups = groupmap->getNamesOfGroups();
+ for (int j = 0; j < mGroups.size(); j++) {
+ if (tree[i].groupCount[mGroups[j]] == 0) {
+ out << 0 << '\t';
+ }else { out << (tree[it->second].groupCount[mGroups[j]] / (float)tree[i].groupCount[mGroups[j]]) << '\t'; }
+ }
+ }else if (ct != NULL) {
+ if (ct->hasGroupInfo()) {
+ vector<string> mGroups = ct->getNamesOfGroups();
+ for (int j = 0; j < mGroups.size(); j++) {
+ if (tree[i].groupCount[mGroups[j]] == 0) {
+ out << 0 << '\t';
+ }else { out << (tree[it->second].groupCount[mGroups[j]] / (float)tree[i].groupCount[mGroups[j]]) << '\t'; }
+ }
+ }
+ }
+ out << endl;
+
+ }
+
+ print(it->second, out, relabund);
+ }
+ }
+ catch(exception& e) {
+ m->errorOut(e, "PhyloSummary", "print");
+ exit(1);
+ }
+}
/**************************************************************************************************/
void PhyloSummary::readTreeStruct(ifstream& in){
try {
for (int j = 0; j < (groupmap->getNamesOfGroups()).size(); j++) {
tree[i].groupCount[(groupmap->getNamesOfGroups())[j]] = 0;
}
- }
+ }else if (ct != NULL) {
+ if (ct->hasGroupInfo()) {
+ for (int j = 0; j < (ct->getNamesOfGroups()).size(); j++) {
+ tree[i].groupCount[(ct->getNamesOfGroups())[j]] = 0;
+ }
+ }
+ }
tree[i].total = 0;