doTaxonomy.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 "doTaxonomy.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
  29 string PhyloTree::getNextTaxon(string& heirarchy){
  30         try {
  31                 string currentLevel = "";
  32                 if(heirarchy != ""){
  33                         currentLevel=heirarchy.substr(0,heirarchy.find_first_of(';'));
  34                         heirarchy=heirarchy.substr(heirarchy.find_first_of(';')+1);
  35                 }
  36                 return currentLevel;
  37         }
  38         catch(exception& e) {
  39                 errorOut(e, "PhyloTree", "getNextTaxon");
  40                 exit(1);
  41         }
  42 }
  43
  44 /**************************************************************************************************/
  45
  46 void PhyloTree::addSeqToTree(string seqName, string seqTaxonomy){
  47         try {
  48                 numSeqs++;
  49
  50                 map<string, int>::iterator childPointer;
  51
  52                 int currentNode = 0;
  53                 int level = 1;
  54
  55                 tree[0].accessions.push_back(seqName);
  56                 string taxon;// = getNextTaxon(seqTaxonomy);
  57
  58                 while(seqTaxonomy != ""){
  59
  60                         level++;
  61
  62                         //somehow the parent is getting one too many accnos
  63                         //use print to reassign the taxa id
  64                         taxon = getNextTaxon(seqTaxonomy);
  65
  66                         childPointer = tree[currentNode].children.find(taxon);
  67
  68                         if(childPointer != tree[currentNode].children.end()){   //if the node already exists, move on
  69                                 currentNode = childPointer->second;
  70                                 tree[currentNode].accessions.push_back(seqName);
  71                                 name2Taxonomy[seqName] = currentNode;
  72                         }
  73                         else{                                                                                   //otherwise, create it
  74                                 tree.push_back(TaxNode(taxon));
  75                                 numNodes++;
  76                                 tree[currentNode].children[taxon] = numNodes-1;
  77                                 tree[numNodes-1].parent = currentNode;
  78
  79                                 //                      int numChildren = tree[currentNode].children.size();
  80                                 //                      string heirarchyID = tree[currentNode].heirarchyID;
  81                                 //                      tree[currentNode].accessions.push_back(seqName);
  82
  83                                 currentNode = tree[currentNode].children[taxon];
  84                                 tree[currentNode].accessions.push_back(seqName);
  85                                 name2Taxonomy[seqName] = currentNode;
  86                                 //                      tree[currentNode].level = level;
  87                                 //                      tree[currentNode].childNumber = numChildren;
  88                                 //                      tree[currentNode].heirarchyID = heirarchyID + '.' + toString(tree[currentNode].childNumber);
  89                         }
  90
  91                         if (seqTaxonomy == "") {   uniqueTaxonomies[currentNode] = currentNode; }
  92                 }
  93
  94         }
  95         catch(exception& e) {
  96                 errorOut(e, "PhyloTree", "addSeqToTree");
  97                 exit(1);
  98         }
  99 }
 100 /**************************************************************************************************/
 101 vector<int> PhyloTree::getGenusNodes()  {
 102         try {
 103                 genusIndex.clear();
 104                 //generate genusIndexes
 105                 map<int, int>::iterator it2;
 106                 for (it2=uniqueTaxonomies.begin(); it2!=uniqueTaxonomies.end(); it2++) {  genusIndex.push_back(it2->first);     }
 107
 108                 return genusIndex;
 109         }
 110         catch(exception& e) {
 111                 errorOut(e, "PhyloTree", "getGenusNodes");
 112                 exit(1);
 113         }
 114 }
 115 /**************************************************************************************************/
 116
 117 void PhyloTree::assignHeirarchyIDs(int index){
 118         try {
 119                 map<string,int>::iterator it;
 120                 int counter = 1;
 121
 122                 for(it=tree[index].children.begin();it!=tree[index].children.end();it++){
 123                         tree[it->second].heirarchyID = tree[index].heirarchyID + '.' + toString(counter);
 124                         counter++;
 125                         tree[it->second].level = tree[index].level + 1;
 126                         assignHeirarchyIDs(it->second);
 127                 }
 128         }
 129         catch(exception& e) {
 130                 errorOut(e, "PhyloTree", "assignHeirarchyIDs");
 131                 exit(1);
 132         }
 133 }
 134
 135 /**************************************************************************************************/
 136
 137 void PhyloTree::print(ofstream& out){
 138         try {
 139                 out << tree[0].level << '\t'<< tree[0].heirarchyID << '\t' << tree[0].name << '\t' << tree[0].children.size() << '\t' << tree[0].accessions.size() << endl;
 140                 print(0, out);
 141         }
 142         catch(exception& e) {
 143                 errorOut(e, "PhyloTree", "print");
 144                 exit(1);
 145         }
 146 }
 147
 148 /**************************************************************************************************/
 149
 150 void PhyloTree::print(int i, ofstream& out){
 151         try {
 152                 map<string,int>::iterator it;
 153                 for(it=tree[i].children.begin();it!=tree[i].children.end();it++){
 154                         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;
 155                         print(it->second, out);
 156                 }
 157         }
 158         catch(exception& e) {
 159                 errorOut(e, "PhyloTree", "print");
 160                 exit(1);
 161         }
 162 }
 163
 164 /**************************************************************************************************/
 165
 166
 167