]> git.donarmstrong.com Git - mothur.git/blob - suffixtree.cpp
test change
[mothur.git] / suffixtree.cpp
1 /*
2  *  suffixtree.cpp
3  *  
4  *
5  *  Created by Pat Schloss on 12/15/08.
6  *  Copyright 2008 Patrick D. Schloss. All rights reserved.
7  *
8  *      This is my half-assed attempt to implement a suffix tree.  This is a cobbled together algorithm using materials that
9  *      I found at http://marknelson.us/1996/08/01/suffix-trees/ and:
10  *
11  *              Ukkonen E. (1995). On-line construction of suffix trees. Algorithmica 14 (3): 249--260
12  *              Gusfield, Dan (1999). Algorithms on Strings, Trees and Sequences: Computer Science and Computational Biology. 
13  *                      USA: Cambridge University Press
14  *
15  *      The Ukkonen paper is the seminal paper describing the on-line method of constructing a suffix tree.
16  *
17  *      I have chosen to store the nodes of the tree as a vector of pointers to SuffixNode objects.  The root is stored at
18  *      nodeVector[0].  Each tree also stores the sequence name and the string that corresponds to the actual sequence. 
19  *      Finally, this class provides a way of counting the number of suffixes that are needed in one tree to generate a new
20  *      sequence (countSuffixes).  This method is used to determine similarity between sequences and was inspired by the
21  *      article and Perl source code provided at http://www.ddj.com/web-development/184416093.
22  *
23  */
24
25 #include "sequence.hpp"
26 #include "suffixnodes.hpp"
27 #include "suffixtree.hpp"
28
29
30 //********************************************************************************************************************
31
32 inline bool compareParents(SuffixNode* left, SuffixNode* right){//      this is necessary to print the tree and to sort the
33         return (left->getParentNode() < right->getParentNode());        //      nodes in order of their parent
34 }
35
36 //********************************************************************************************************************
37
38 SuffixTree::SuffixTree(){ m = MothurOut::getInstance(); }
39
40 //********************************************************************************************************************
41
42 SuffixTree::~SuffixTree(){
43         for(int i=0;i<nodeVector.size();i++){   delete nodeVector[i];   }       
44         nodeVector.clear();
45 }
46
47 //********************************************************************************************************************
48
49 void SuffixTree::loadSequence(Sequence seq){
50         nodeCounter = 0;                                                        //      initially there are 0 nodes in the tree
51         activeStartPosition = 0;
52         activeEndPosition = -1;                                         
53         seqName = seq.getName();
54         sequence = seq.convert2ints();
55         sequence += '5';                                                        //      this essentially concatenates a '$' to the end of the sequence to
56         int seqLength = sequence.length();                      //      make it a cononical suffix tree
57         
58         nodeVector.push_back(new SuffixBranch(-1, 0, -1));      //      enter the root of the suffix tree
59         
60         activeNode = root = 0;
61         string hold;
62         for(int i=0;i<seqLength;i++){
63                 addPrefix(i);                                                   //      step through the sequence adding each prefix
64         }
65 }
66
67 //********************************************************************************************************************
68
69 string SuffixTree::getSeqName() {
70         return seqName;         
71 }
72
73 //********************************************************************************************************************
74
75 void SuffixTree::print(){
76         vector<SuffixNode*> hold = nodeVector;
77         sort(hold.begin(), hold.end(), compareParents);
78         m->mothurOut("Address\t\tParent\tNode\tSuffix\tStartC\tEndC\tSuffix"); m->mothurOutEndLine();
79         for(int i=1;i<=nodeCounter;i++){
80                 hold[i]->print(sequence, i);
81         }
82 }
83
84 //********************************************************************************************************************
85
86 int SuffixTree::countSuffixes(string compareSequence, int& minValue){   //      here we count the number of suffix parts 
87                                                                                                                         //      we need to rewrite a user supplied sequence.  if the 
88         int numSuffixes = 0;                                                                    //      count exceeds the supplied minValue, bail out.  The
89         int seqLength = compareSequence.length();                               //      time complexity should be O(L)
90         int position = 0;
91         
92         int presentNode = 0;
93         
94         while(position < seqLength){            //      while the position in the query sequence isn't at the end...
95                 
96                 if(numSuffixes > minValue)      {       return 1000000;         }       //      bail if the count gets too high
97                 
98                 int newNode = nodeVector[presentNode]->getChild(compareSequence[position]);     //      see if the current node has a
99                                                                                                                                 //      child that matches the next character in the query
100                 if(newNode == -1){                                                                              
101                         if(presentNode == 0){   position++;             }                       //      if not, go back to the root and increase the count
102                         numSuffixes++;                                                                          //      by one.
103                         presentNode = 0;
104                 }
105                 else{                                                                                                   //      if there is, move to that node and see how far down
106                         presentNode = newNode;                                                          //      it we can get
107                         
108                         for(int i=nodeVector[newNode]->getStartCharPos(); i<=nodeVector[newNode]->getEndCharPos(); i++){
109                                 if(compareSequence[position] == sequence[i]){
110                                         position++;                                                                     //      as long as the query and branch agree, keep going
111                                 }
112                                 else{
113                                         numSuffixes++;                                                          //      if there is a mismatch, increase the number of 
114                                         presentNode = 0;                                                        //      suffixes and go back to the root
115                                         break;
116                                 }
117                         }
118                 }
119                 //      if we get all the way through the node we'll go to the top of the while loop and find the child node
120                 //      that corresponds to what we are interested in           
121         }
122         numSuffixes--;                                                                                          //      the method puts an extra count on numSuffixes
123         
124         if(numSuffixes < minValue)      {       minValue = numSuffixes; }       //      if the count is less than the previous minValue,
125         return numSuffixes;                                                                                     //      change the value and return the number of suffixes
126         
127 }
128 //********************************************************************************************************************
129
130 int SuffixTree::countSuffixes(string compareSequence){  //      here we count the number of suffix parts 
131                                                                                                                         //      we need to rewrite a user supplied sequence.  if the 
132         int numSuffixes = 0;                                                                    //      count exceeds the supplied minValue, bail out.  The
133         int seqLength = compareSequence.length();                               //      time complexity should be O(L)
134         int position = 0;
135         
136         int presentNode = 0;
137         
138         while(position < seqLength){            //      while the position in the query sequence isn't at the end...
139                 
140                 int newNode = nodeVector[presentNode]->getChild(compareSequence[position]);     //      see if the current node has a
141                                                                                                                                 //      child that matches the next character in the query
142                 if(newNode == -1){                                                                              
143                         if(presentNode == 0){   position++;             }                       //      if not, go back to the root and increase the count
144                         numSuffixes++;                                                                          //      by one.
145                         presentNode = 0;
146                 }
147                 else{                                                                                                   //      if there is, move to that node and see how far down
148                         presentNode = newNode;                                                          //      it we can get
149                         
150                         for(int i=nodeVector[newNode]->getStartCharPos(); i<=nodeVector[newNode]->getEndCharPos(); i++){
151                                 if(compareSequence[position] == sequence[i]){
152                                         position++;                                                                     //      as long as the query and branch agree, keep going
153                                 }
154                                 else{
155                                         numSuffixes++;                                                          //      if there is a mismatch, increase the number of 
156                                         presentNode = 0;                                                        //      suffixes and go back to the root
157                                         break;
158                                 }
159                         }
160                 }
161                 //      if we get all the way through the node we'll go to the top of the while loop and find the child node
162                 //      that corresponds to what we are interested in           
163         }
164         numSuffixes--;                                                                                          //      the method puts an extra count on numSuffixes
165         
166         return numSuffixes;                                                                                     //      change the value and return the number of suffixes
167 }
168 //********************************************************************************************************************
169
170 void SuffixTree::canonize(){    //      if you have to ask how this works, you don't really want to know and this really
171                                                                 //      isn't the place to ask.
172         if ( isExplicit() == 0 ) {      //      if the node has no children...
173                 
174                 int tempNodeIndex = nodeVector[activeNode]->getChild(sequence[activeStartPosition]);
175                 SuffixNode* tempNode = nodeVector[tempNodeIndex];
176                 
177                 int span = tempNode->getEndCharPos() - tempNode->getStartCharPos();
178                 
179                 while ( span <= ( activeEndPosition - activeStartPosition ) ) {
180                         
181             activeStartPosition = activeStartPosition + span + 1;
182                         
183                         activeNode = tempNodeIndex;
184                         
185             if ( activeStartPosition <= activeEndPosition ) {
186                                 tempNodeIndex = nodeVector[tempNodeIndex]->getChild(sequence[activeStartPosition]);
187                                 tempNode = nodeVector[tempNodeIndex];
188                                 span = tempNode->getEndCharPos() - tempNode->getStartCharPos();
189             }
190                         
191         }
192     }
193 }
194
195 //********************************************************************************************************************
196
197 int SuffixTree::split(int nodeIndex, int position){     //      leaves stay leaves, etc, to split a leaf we make a new interior 
198                                                                                                         //      node and reconnect everything
199         SuffixNode* node = nodeVector[nodeIndex];                                       //      get the node that needs to be split
200         SuffixNode* parentNode = nodeVector[node->getParentNode()];     //      get it's parent node
201         
202         parentNode->eraseChild(sequence[node->getStartCharPos()]);      //      erase the present node from the registry of its parent
203         
204         nodeCounter++;
205         SuffixNode* newNode = new SuffixBranch(node->getParentNode(), node->getStartCharPos(), node->getStartCharPos() + activeEndPosition - activeStartPosition);      //      create a new node that will link the parent with the old child
206         parentNode->setChildren(sequence[newNode->getStartCharPos()], nodeCounter);//   give the parent the new child
207         nodeVector.push_back(newNode);
208         
209         node->setParentNode(nodeCounter);       //      give the original node the new node as its parent
210         newNode->setChildren(sequence[node->getStartCharPos() + activeEndPosition - activeStartPosition + 1], nodeIndex);
211         //      put the original node in the registry of the new node's children
212         newNode->setSuffixNode(activeNode);//link the new node with the old active node
213         
214         //      recalculate the startCharPosition of the outermost node
215         node->setStartCharPos(node->getStartCharPos() + activeEndPosition - activeStartPosition + 1 );
216         
217         return node->getParentNode();
218 }
219
220 //********************************************************************************************************************
221
222 void SuffixTree::makeSuffixLink(int& previous, int present){
223         
224 //      here we link the nodes that are suffixes of one another to rapidly speed through the tree
225         if ( previous > 0 ) {   nodeVector[previous]->setSuffixNode(present);   }
226         else                            {       /*      do nothing                                                              */      }
227         
228     previous = present;
229 }
230
231 //********************************************************************************************************************
232
233 void SuffixTree::addPrefix(int prefixPosition){
234         
235         int lastParentNode = -1;        //      we need to place a new prefix in the suffix tree
236         int parentNode = 0;
237         
238         while(1){
239                 
240                 parentNode = activeNode;
241                 
242                 if(isExplicit() == 1){  //      if the node is explicit (has kids), try to follow it down the branch if its there...
243                         if(nodeVector[activeNode]->getChild(sequence[prefixPosition]) != -1){   //      break out and get next prefix...
244                                 break;                                                                                          
245                         }
246                         else{                           //      ...otherwise continue, we'll need to make a new node later on...
247                         }
248                 }
249                 else{                                   //      if it's not explicit (no kids), read through and see if all of the chars agree...
250                         int tempNode = nodeVector[activeNode]->getChild(sequence[activeStartPosition]);
251                         int span = activeEndPosition - activeStartPosition;
252                         
253                         if(sequence[nodeVector[tempNode]->getStartCharPos() + span + 1] == sequence[prefixPosition] ){
254                                 break;                  //      if the existing suffix agrees with the new one, grab a new prefix...
255                         }
256                         else{
257                                 parentNode = split(tempNode, prefixPosition);   //      ... otherwise we need to split the node
258                         }
259                         
260                 }
261                 
262                 nodeCounter++;  //      we need to generate a new node here if the kid didn't exist, or we split a node
263                 SuffixNode* newSuffixLeaf = new SuffixLeaf(parentNode, prefixPosition, sequence.length()-1);
264                 nodeVector[parentNode]->setChildren(sequence[prefixPosition], nodeCounter);
265                 nodeVector.push_back(newSuffixLeaf);
266                 
267                 makeSuffixLink( lastParentNode, parentNode );           //      make a suffix link for the parent node
268                 
269                 if(nodeVector[activeNode]->getParentNode() == -1){      //      move along the start position for the tree
270             activeStartPosition++;
271         } 
272                 else {
273             activeNode = nodeVector[activeNode]->getSuffixNode();
274                 }
275                 canonize();                                                                                     //      frankly, i'm not entirely clear on what canonize does.
276         }
277         
278         makeSuffixLink( lastParentNode, parentNode );
279         activeEndPosition++;                                                                    //      move along the end position for the tree
280         
281         canonize();                                                                                             //      frankly, i'm not entirely clear on what canonize does.
282         
283 }
284
285 //********************************************************************************************************************
286