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
  16         numNodes = 1;
  17         numSeqs = 0;
  18         tree.push_back(TaxNode("Root"));
  19         tree[0].heirarchyID = "0";
  20 }
  21
  22 /**************************************************************************************************/
  23
  24 string PhyloTree::getNextTaxon(string& heirarchy){
  25
  26         string currentLevel = "";
  27         if(heirarchy != ""){
  28                 currentLevel=heirarchy.substr(0,heirarchy.find_first_of(';'));
  29                 heirarchy=heirarchy.substr(heirarchy.find_first_of(';')+1);
  30         }
  31         return currentLevel;
  32
  33 }
  34
  35 /**************************************************************************************************/
  36
  37 void PhyloTree::addSeqToTree(string seqName, string seqTaxonomy){
  38
  39         numSeqs++;
  40
  41         map<string, int>::iterator childPointer;
  42
  43         int currentNode = 0;
  44         int level = 1;
  45
  46         tree[0].accessions.push_back(seqName);
  47         string taxon;// = getNextTaxon(seqTaxonomy);
  48
  49         while(seqTaxonomy != ""){
  50
  51                 level++;
  52
  53 //somehow the parent is getting one too many accnos
  54 //use print to reassign the taxa id
  55                 taxon = getNextTaxon(seqTaxonomy);
  56
  57                 childPointer = tree[currentNode].children.find(taxon);
  58
  59                 if(childPointer != tree[currentNode].children.end()){   //if the node already exists, move on
  60                         currentNode = childPointer->second;
  61                         tree[currentNode].accessions.push_back(seqName);
  62                 }
  63                 else{                                                                                   //otherwise, create it
  64                         tree.push_back(TaxNode(taxon));
  65                         numNodes++;
  66                         tree[currentNode].children[taxon] = numNodes-1;
  67
  68 //                      int numChildren = tree[currentNode].children.size();
  69 //                      string heirarchyID = tree[currentNode].heirarchyID;
  70 //                      tree[currentNode].accessions.push_back(seqName);
  71
  72                         currentNode = tree[currentNode].children[taxon];
  73                         tree[currentNode].accessions.push_back(seqName);
  74
  75 //                      tree[currentNode].level = level;
  76 //                      tree[currentNode].childNumber = numChildren;
  77 //                      tree[currentNode].heirarchyID = heirarchyID + '.' + toString(tree[currentNode].childNumber);
  78                 }
  79
  80         }
  81 }
  82
  83 /**************************************************************************************************/
  84
  85 void PhyloTree::assignHeirarchyIDs(int index){
  86
  87         map<string,int>::iterator it;
  88         int counter = 1;
  89
  90         for(it=tree[index].children.begin();it!=tree[index].children.end();it++){
  91                 tree[it->second].heirarchyID = tree[index].heirarchyID + '.' + toString(counter);
  92                 counter++;
  93                 tree[it->second].level = tree[index].level + 1;
  94                 assignHeirarchyIDs(it->second);
  95
  96         }
  97
  98 }
  99
 100 /**************************************************************************************************/
 101
 102 void PhyloTree::print(ofstream& out){
 103
 104         out << tree[0].level << '\t'<< tree[0].heirarchyID << '\t' << tree[0].name << '\t' << tree[0].children.size() << '\t' << tree[0].accessions.size() << endl;
 105         print(0, out);
 106
 107
 108 }
 109
 110 /**************************************************************************************************/
 111
 112 void PhyloTree::print(int i, ofstream& out){
 113
 114         map<string,int>::iterator it;
 115         for(it=tree[i].children.begin();it!=tree[i].children.end();it++){
 116                 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;
 117                 print(it->second, out);
 118         }
 119
 120 }
 121
 122 /**************************************************************************************************/
 123
 124
 125