phylotree.cpp

   1 /*
   2  *  doTaxonomy.cpp
   3  *
   4  *
   5  *  Created by Pat Schloss on 6/17/09.
   6  *  Copyright 2009 Patrick D. Schloss. All rights reserved.
   7  *
   8  */
   9
  10 #include "phylotree.h"
  11
  12 /**************************************************************************************************/
  13
  14 PhyloTree::PhyloTree(){
  15         try {
  16                 numNodes = 1;
  17                 numSeqs = 0;
  18                 tree.push_back(TaxNode("Root"));
  19                 tree[0].heirarchyID = "0";
  20         }
  21         catch(exception& e) {
  22                 errorOut(e, "PhyloTree", "PhyloTree");
  23                 exit(1);
  24         }
  25 }
  26 /**************************************************************************************************/
  27
  28 PhyloTree::PhyloTree(string tfile){
  29         try {
  30                 numNodes = 1;
  31                 numSeqs = 0;
  32                 tree.push_back(TaxNode("Root"));
  33                 tree[0].heirarchyID = "0";
  34
  35                 ifstream in;
  36                 openInputFile(tfile, in);
  37
  38                 //read in users taxonomy file and add sequences to tree
  39                 string name, tax;
  40                 while(!in.eof()){
  41                         in >> name >> tax; gobble(in);
  42
  43                         addSeqToTree(name, tax);
  44                 }
  45                 in.close();
  46
  47                 assignHeirarchyIDs(0);
  48         }
  49         catch(exception& e) {
  50                 errorOut(e, "PhyloTree", "PhyloTree");
  51                 exit(1);
  52         }
  53 }
  54
  55 /**************************************************************************************************/
  56
  57 string PhyloTree::getNextTaxon(string& heirarchy){
  58         try {
  59                 string currentLevel = "";
  60                 if(heirarchy != ""){
  61                         currentLevel=heirarchy.substr(0,heirarchy.find_first_of(';'));
  62                         heirarchy=heirarchy.substr(heirarchy.find_first_of(';')+1);
  63                 }
  64                 return currentLevel;
  65         }
  66         catch(exception& e) {
  67                 errorOut(e, "PhyloTree", "getNextTaxon");
  68                 exit(1);
  69         }
  70 }
  71
  72 /**************************************************************************************************/
  73
  74 void PhyloTree::addSeqToTree(string seqName, string seqTaxonomy){
  75         try {
  76                 numSeqs++;
  77
  78                 map<string, int>::iterator childPointer;
  79
  80                 int currentNode = 0;
  81                 int level = 1;
  82
  83                 tree[0].accessions.push_back(seqName);
  84                 string taxon;// = getNextTaxon(seqTaxonomy);
  85
  86                 while(seqTaxonomy != ""){
  87
  88                         level++;
  89
  90                         //somehow the parent is getting one too many accnos
  91                         //use print to reassign the taxa id
  92                         taxon = getNextTaxon(seqTaxonomy);
  93
  94                         childPointer = tree[currentNode].children.find(taxon);
  95
  96                         if(childPointer != tree[currentNode].children.end()){   //if the node already exists, move on
  97                                 currentNode = childPointer->second;
  98                                 tree[currentNode].accessions.push_back(seqName);
  99                                 name2Taxonomy[seqName] = currentNode;
 100                         }
 101                         else{                                                                                   //otherwise, create it
 102                                 tree.push_back(TaxNode(taxon));
 103                                 numNodes++;
 104                                 tree[currentNode].children[taxon] = numNodes-1;
 105                                 tree[numNodes-1].parent = currentNode;
 106
 107                                 //                      int numChildren = tree[currentNode].children.size();
 108                                 //                      string heirarchyID = tree[currentNode].heirarchyID;
 109                                 //                      tree[currentNode].accessions.push_back(seqName);
 110
 111                                 currentNode = tree[currentNode].children[taxon];
 112                                 tree[currentNode].accessions.push_back(seqName);
 113                                 name2Taxonomy[seqName] = currentNode;
 114                                 //                      tree[currentNode].level = level;
 115                                 //                      tree[currentNode].childNumber = numChildren;
 116                                 //                      tree[currentNode].heirarchyID = heirarchyID + '.' + toString(tree[currentNode].childNumber);
 117                         }
 118
 119                         if (seqTaxonomy == "") {   uniqueTaxonomies[currentNode] = currentNode; }
 120                 }
 121
 122         }
 123         catch(exception& e) {
 124                 errorOut(e, "PhyloTree", "addSeqToTree");
 125                 exit(1);
 126         }
 127 }
 128 /**************************************************************************************************/
 129 vector<int> PhyloTree::getGenusNodes()  {
 130         try {
 131                 genusIndex.clear();
 132                 //generate genusIndexes
 133                 map<int, int>::iterator it2;
 134                 for (it2=uniqueTaxonomies.begin(); it2!=uniqueTaxonomies.end(); it2++) {  genusIndex.push_back(it2->first);     }
 135
 136                 return genusIndex;
 137         }
 138         catch(exception& e) {
 139                 errorOut(e, "PhyloTree", "getGenusNodes");
 140                 exit(1);
 141         }
 142 }
 143 /**************************************************************************************************/
 144
 145 void PhyloTree::assignHeirarchyIDs(int index){
 146         try {
 147                 map<string,int>::iterator it;
 148                 int counter = 1;
 149
 150                 for(it=tree[index].children.begin();it!=tree[index].children.end();it++){
 151                         tree[it->second].heirarchyID = tree[index].heirarchyID + '.' + toString(counter);
 152                         counter++;
 153                         tree[it->second].level = tree[index].level + 1;
 154                         assignHeirarchyIDs(it->second);
 155                 }
 156         }
 157         catch(exception& e) {
 158                 errorOut(e, "PhyloTree", "assignHeirarchyIDs");
 159                 exit(1);
 160         }
 161 }
 162
 163 /**************************************************************************************************/
 164
 165 void PhyloTree::print(ofstream& out){
 166         try {
 167                 out << tree[0].level << '\t'<< tree[0].heirarchyID << '\t' << tree[0].name << '\t' << tree[0].children.size() << '\t' << tree[0].accessions.size() << endl;
 168                 print(0, out);
 169         }
 170         catch(exception& e) {
 171                 errorOut(e, "PhyloTree", "print");
 172                 exit(1);
 173         }
 174 }
 175
 176 /**************************************************************************************************/
 177
 178 void PhyloTree::print(int i, ofstream& out){
 179         try {
 180                 map<string,int>::iterator it;
 181                 for(it=tree[i].children.begin();it!=tree[i].children.end();it++){
 182                         out <<tree[it->second].level << '\t' << tree[it->second].heirarchyID << '\t' << tree[it->second].name << '\t' << tree[it->second].children.size() << '\t' << tree[it->second].accessions.size() << endl;
 183                         print(it->second, out);
 184                 }
 185         }
 186         catch(exception& e) {
 187                 errorOut(e, "PhyloTree", "print");
 188                 exit(1);
 189         }
 190 }
 191
 192 /**************************************************************************************************/
 193
 194
 195