Common substrings of more than 2 strings (VI) By traversing T, for each internal node v, compute C(v). [C(v) 3. is defined as the number of distinct termination symbols in the subtree rooted at v] This step takes O(Kn) time. 3 c a % # $ g 3 3 4 5 5 g c g a 3 3 2 c g # % a c a c c $ $ # $ $ # # % % 4 3 3 1 2 1 4 2 3 2 1
Common substrings of more than 2 strings (VII) Traverse T and visit every internal node v. For each v, if V(C(v)) < 4. string-depth of v, set V(C(v)) = string-depth of v. [After step 4, V(k) = the length of the longest substring common to exactly k of these strings.] l(k)=V(k). For i=k-1 downto 2, l(i)=max{l(i+1), V(i)} . 5. This two steps take O(n) time. For our example, V(2) = 3, V(3) = 2. Thus, l(3) = 2, l(2) = 3. In total, this algorithm takes O(Kn) time. Actually, we can improve this algorithm to O(n) time by mean of lcp!
Linear time algorithm for constructing suffix tree
Straightforward construction of suffix tree Consider S = s 1 s 2 … s n where s n =$ Algorithm: Initialize the tree with only a root For i = n to 1 Includes S[i..n] into the tree Time: O(n 2 )
Example of construction S=acca$ Init For-loop c c a c a a a $ a $ $ $ c $ $ $ c c $ $ a c a a a $ $ a $ $ $ $ 4 4 5 5 5 3 4 3 2 4 1 3 2 5 5 I 4 I 3 I 5 I 2 I 1
Construction of generalized suffix tree S ’ = c# Init For-loop c c c a a a c c c $ # $ # $ c c c c a c a c a a a a # $ $ a $ $ a $ $ a $ $ $ $ $ $ 4 1 3 2 2 4 1 3 2 2 4 1 3 2 1 5 5 5 I 1 J 2 J 1
Can we construct a suffix tree in o(n 2 ) time? Yes. We can construct it in O(n) time. Weiner ’ s algorithm [1973] Linear time for constant size alphabet, but much space McGreight ’ s algorithm [JACM 1976] Linear time for constant size alphabet, quadratic space Ukkonen ’ s algorithm [Algorithmica, 1995] Online algorithm, linear time for constant size alphabet, less space Farach ’ s algorithm [FOCS 1997] Linear time for general alphabet Hon,Sadakane, and Sung ’ s algorithm [FOCS 2003] O(n) bit space O(n log e n) time for 0<e<1 O(n) bit space O(n) time for suffix array construction We will discuss Farach ’ s algorithm later.
Idea Build Odd Suffix Tree and Even Suffix Tree Then, merge odd and even suffix tree. $ a g $ a g ca ca $ $ 7 6 7 6 g g c c a a $ $ c c a a g g g g 5 5 $ $ c c $ $ a a g g g g $ $ 4 4 $ $ 2 2 1 3 1 3 Even Suffix Tree Odd Suffix Tree
Idea Input: a string S of length n Recursively compute the suffix tree T o of all 1. suffixes beginning at the odd positions. T o is of size n/2. From T o , compute T e which is the suffix tree 2. for all suffixes beginning at the even positions. Merge T o and T e to form the suffix tree for 3. S.
Stage 1: Constructing odd suffix tree Given a string S[1..n], we generate a new string S ’ [1..n/2] as follows. we map pairs of characters into single characters as follows: S[1..2], S[3..4], S[5..6], … , S[n-1..n]. Remove the duplicates from the pairs of characters and sort them by radix sort. S ’ [i] = rank of S[2i-1..2i] in the sorted list, for i=1, 2, … , n/2. By recursion, we get the suffix tree T ’ for S ’ Convert T ’ to the odd suffix tree T o .
Example (I) S = aaabbbabbaba$ S[1..2]=aa, S[3..4]=ab, S[5..6]=bb, S[7..8]=ab, S[9..10]=ba, S[11..12]=ba. By stable sort, aa < ab < ba < bb. Rank(aa)=1, Rank(ab)=2, Rank(ba)=3, Rank(bb)=4. So, S ’ =124233$.
Example (II) By recursion, construct the suffix tree T ’ for S ’ : $ 2 3 7 3 $ 6 5 1 4 2
Example (III) Convert T ’ to the odd tree: $ a 13 5 b i 2i-1 $ 11 9 7 1 This is not a suffix tree 3
Example (IV) Refine the odd tree T o : $ b a 13 b 5 a b $ 11 9 7 1 3
Time complexity for building the odd tree Let Time(n) be the time to build a suffix tree for a string of length n. Stable sorting and refinement of the odd trees take O(n) time. Build suffix tree for S ’ takes Time(n/2). So, Stage 1 takes Time(n/2)+O(n) time.
Stage 2: Build the even tree Generate the lex-ordering of the 1. leaves in T e . For any two adjacent leaves 2i and 2j, 2. we find lcp(2i, 2j). Construct the even tree T e from left to 3. right (according to the lex-ordering).
Build the even tree (Step 1) We get the lex-ordering of the leaves in T o . Generate the lex-ordering of the leaves in T e . For each leaf i in T o , get the preceding character c=S[i-1] and form a pair (c,i). Each pair represents a even suffix i-1. Perform stable sorting on those pairs. We get the lex-ordering of the leaves in T e .
Example S = aaabbbabbaba$ Lex-ordering of the leaves in T o : 13 < 1 < 7 < 3 < 11 < 9 < 5 The pairs are: (a,13), ($,1), (b,7), (a, 3), (a, 11), (b, 9), (b, 5). After stable sorting, we have ($, 1), (a, 13), (a, 3), (a, 11), (b, 7), (b, 9), (b, 5). Hence, the lex-ordering of the leaves of T e : 12 < 2 < 10 < 6 < 8 < 4
Build the even tree (Step 2) For any two adjacent leaves 2i and 2j, we first find lcp(2i, 2j). Observation: lcp(2i, 2j) = lcp(2i+1, 2j+1)+1 if S[2i]=S[2j] 0 otherwise Proof: If S[2i] ≠ S[2j], lcp(2i,2j)=0. Otherwise, lcp(2i,2j)=1+lcp(2i+1,2j+1).
Example Recall that the lex- $ b a ordering of leaves: 13 b 12 < 2 < 10 < 6 < 8 < 4. 5 a b By the previous observation, we have $ lcp(8,4)=lcp(9,5)+1=2 Similarly, we have 11 9 lcp(12,2)=1, lcp(2,10)=1, 7 lcp(10,6)=0, lcp(6,8)=1, 1 lcp(8,4)=2 3
Build the even tree (Step 3) Construct the even tree T e from left to right. a a $ $ a a b b 12 12 b b b b a 12 a 10 b b b b a a b b a a $ $ 2 2
Build the even tree (Step 3) b b a a a b $ $ $ a a a b b b 12 12 12 b b b b b b a a a 10 10 8 10 b b b 8 b b b 6 6 a a a b b b 6 a a a 4 $ $ $ 2 2 2
Time complexity for building the even tree Step 1: O(n) time Step 2: O(n) time Step 3: O(n) time
Stage 3: Merge odd and even trees We can merge T o and T e by DFS. However, it takes O(n 2 ) time. $ $ b b a a babbaba$ babbaba$ 13 13 b b a,1 a,1 b,1 b,1 a a 5 5 $,1 $,1 a a b b a a b b 13 13 b,1 b,1 b b a,1 a,1 b,1 b,1 b b b b b b b b $,1 $,1 a a b b $ $ a a a a a a b b a a b b a,1 a,1 a,2 a,2 12 12 $ $ b b a a b b b b b,1 b,1 a a b b $ $ a,2 a,2 $ $ b,1 b,1 b b 11 11 9 9 a a a a $,1 $,1 a,2 a,2 a a b b 10 10 b,8 b,8 b b $ $ 8 8 a a 12 12 7 7 a,2 a,2 b b $ $ a,11 1 1 b b 11 11 9 9 b,5 b,5 a a 1 1 3 3 10 10 b,10 b,5 b,5 b b 4 4 8 8 a,4 b,9 a,4 b,9 Odd tree 5 5 b b a a 7 7 3 3 6 6 2 2 b b 6 6 a a 4 4 $ $ 2 2 Even tree
Stage 3: Merge odd and even trees We merge T o and T e by DFS. We merge two edges as long as they start with the same character. The merge is ended when one edge is longer than the other. $ b a b a 13 5 a b $ a b 12 b $ b a 10 b 8 b 11 9 a b 6 7 a 4 1 $ 2 3
Merge odd and even trees We merge T o and T e by DFS. We merge two edges as long as they start with the same character. The merge is ended when one edge is longer than the other. a,1 b,1 $,1 b,1 13 a,1 b,2 $,1 12 a,4 a,11 $,1 $,1 b,3 b,8 10 8 a,4 11 9 b,3 2 b,9 4 7 b,3 5 $,1 6 1 3
Merge odd and even trees The merging may over-merged some nodes. To correct the tree, we need to unmerge some nodes. a,1 b,1 $,1 b,1 13 a,1 b,2 $,1 12 a,4 a,11 $,1 $,1 b,3 b,8 10 8 a,4 11 9 b,3 2 b,9 4 7 b,3 5 $,1 6 1 3
Definition of L() and d() For every node u which may be over-merged, there exist two leaves 2i and 2j-1 such that u=lca(2i, 2j-1). Denote L(u) be the correct depth of u, that is, lcp(2i,2j-1). Note that lcp(2i,2j-1) = 1+lcp(2i+1,2j) if S[2i]=S[2j-1]; 0 otherwise. Let v be lca(2i+1,2j). Denote d(u) = v. Note that d() is equivalent to suffix link!
Example of d() a,1 b,1 $,1 b,1 13 a,1 b,2 $,1 12 a,4 a,11 d() $,1 $,1 b,3 b,8 10 8 a,4 11 9 b,3 2 b,9 4 7 b,3 5 $,1 6 1 3
Relationship between L() and d() Suppose u = lca(2i, 2j-1). Note1: if u is not the root, then S[2i]=S[2j-1]. Note2: lcp(2i,2j-1) = 1+lcp(2i+1,2j) if S[2i]=S[2j-1]; 0 otherwise Note3: d(u) = lcp(2i+1,2j) Hence, L(u) = 1 + L(d(u)) if u is not the root. Otherwise, L(u)=0. Lemma: L(u) = the length of the purple path from u to the root.
Example of L() L( )=1 a,1 b,1 $,1 L( )=1 b,1 13 a,1 b,2 $,1 L( )=2 L( )=2 12 L( )=2 a,4 a,11 $,1 $,1 b,3 b,8 10 8 L( )=4 a,4 11 9 b,3 2 b,9 L( )=2 L( )=3 4 7 b,3 5 $,1 6 1 3
Unmerge the border nodes based on L() (I) L( )=1 a,1 b,1 $,1 L( )=1 b,1 13 a,1 b,2 $,1 L( )=2 L( )=2 12 L( )=2 a,4 a,11 $,1 $,1 b,3 b,8 10 8 L( )=4 a,4 11 9 b,3 2 b,9 L( )=2 L( )=3 4 7 b,3 5 $,1 6 1 3
Unmerge the border nodes based on L() (II) a,1 b,1 $,1 b,1 13 a,1 b,1 $,1 a,1 a,2 12 b,1 a,2 b,1 $,1 a,2 b,8 10 8 a,2 a,11 11 9 b,5 1 b,10 b,5 b,9 a,4 4 5 7 3 6 2
Time complexity for merging Merge the tree using DFS takes O(n) time. Compute the links d() takes O(n) time. Compute L() takes O(n) time. Unmerge takes O(n) time.
Total time complexity of Farach ’ s algorithm Stage 1: Time(n/2)+O(n) Stage 2: O(n) Stage 3: O(n) Thus, Time(n) = Time(n/2)+O(n). By solving the equation, Time(n)= O(n).
Disadvantage of suffix tree Suffix tree is space inefficient. It requires O(n| Σ |log n) bits. Manber and Myers (SIAM J. Comp 1993) proposes a new data structure, called suffix array, which has a similar functionality as suffix tree. Moreover, it only requires O(n log n) bits.
Suffix Array (I) It is just sorted suffixes. E.g. consider S = acacag$ Suffix Position SA[i] Suffix acacag$ 1 7 $ cacag$ 2 1 acacag$ acag$ 3 => 3 acag$ 5 cag$ 4 Sort ag$ ag$ 5 2 cacag$ g$ 6 4 cag$ $ 7 6 g$ Suffix array is an array of n indices. Thus, it takes O(n log n) bits.
Observation The leaves of a suffix tree is in suffix array order. SA[i] Suffix $ a g ca 7 $ $ 7 6 1 acacag$ g c 3 acag$ a $ c 5 ag$ a g g 2 cacag$ 5 $ c $ a g 4 cag$ g $ 2 4 $ 6 g$ 1 3
Linear time construction of suffix array from suffix tree Recall that the suffix tree T of S[1..n] can be constructed in O(n) time. Then, by “ lexical ” depth-first traversal of T , the suffix array of S is obtained. This takes O(n) time. However, the space used during construction is the same as that for suffix tree! This defeats the purpose of suffix array. Today, we can build suffix array using O(n) bit space and O(n) time.
range(T,Q) For a pattern Q, its SA[i] Suffix occurrences in T form a 1 7 $ consecutive SA range. 2 1 acacag$ 3 3 acag$ Example: For T= acacag$, ca 4 5 ag$ occurs in SA[5] and SA[6]. 5 2 cacag$ Definition: 6 4 cag$ We called range(T,Q)=[st..ed] 7 6 g$ if Q is a prefix of every T j for j=SA[st], SA[st+1], … , SA[ed] where T j = j suffix of T = T[j..n]. Example: range(T,ca)=[5..6]
Find occurrence of query Q in a string S using suffix array Input: (1) the suffix array of a string T of length n and (2) a query Q of length m Aim: check if Q occurs in T Idea: binary search!
Algorithm
Example Consider T = acacag$ i SA[i] Suffix Pattern Q = acag 1 7 $ 2 1 acacag$ 3 3 acag$ L=1 4 5 ag$ 5 2 cacag$ R=7 6 4 cag$ 7 6 g$ M=(L+R)/2=4
Example Consider T = acacag$ i SA[i] Suffix Pattern Q = acag 1 7 $ 2 1 acacag$ L=1 3 3 acag$ 4 5 R=7 ag$ 5 2 cacag$ M=(L+R)/2=4 6 4 cag$ 7 6 g$ suffix-SA[M] > Q. Set R=M=4.
Example Consider T = acacag$ i SA[i] Suffix Pattern Q = acag 1 7 $ 2 1 acacag$ L=1 3 3 acag$ 4 5 R=4 ag$ 5 2 cacag$ M=(L+R)/2=2 6 4 cag$ 7 6 g$ suffix-SA[M] < Q. Set L=M=2.
Example Consider T = acacag$ i SA[i] Suffix Pattern Q = acag 1 7 $ 2 1 acacag$ 3 3 acag$ L=2 4 5 ag$ 5 2 cacag$ R=4 6 4 cag$ M=(L+R)/2=3 7 6 g$ The pattern Q is found at SA[M]=3.
Can we do better? During each step of binary search, we need to compare Q with a suffix using O(m) time, which is time consuming. Can we do better? We have the following observation. Suppose LCP(Q, suffix-SA[L]) is l and LCP(Q, suffix-SA[R]) is r. Then, LCP(Q, suffix-SA[M]) > min{l,r}. Below, we describe how to utilize this observation to speedup the computation.
Algorithm
Example Consider T = acacag$ i SA[i] Suffix Pattern Q = acag 1 7 $ 2 1 acacag$ 3 3 acag$ L=1, l=0 4 5 ag$ 5 2 cacag$ R=7, r=0 6 4 cag$ 7 6 g$ mlr = min(l,r)=0 M=(L+R)/2=4
Example Consider T = acacag$ i SA[i] Suffix Pattern Q = acag 1 7 $ L=1, l=0 2 1 acacag$ R=7, r=0 3 3 acag$ mlr = min(l,r)=0 4 5 ag$ M=(L+R)/2=4, m=1 5 2 cacag$ 6 4 cag$ The (m+1) char of suffix-SA[M] is g. The (m+1) char of Q is c. 7 6 g$ So, suffix-SA[M] > Q. Set R=M=4 and r=m=1.
Example Consider T = acacag$ i SA[i] Suffix Pattern Q = acag 1 7 $ L=1, l=0 2 1 acacag$ R=4, r=1 3 3 acag$ mlr = min(l,r)=0 4 5 ag$ M=(L+R)/2=2, m=3 5 2 cacag$ 6 4 cag$ The (m+1) char of suffix-SA[M] is c. The (m+1) char of Q is g. 7 6 g$ So, suffix-SA[M] < Q. Set L=M=2 and l=m=3.
Example Consider T = acacag$ i SA[i] Suffix Pattern Q = acag 1 7 $ 2 1 acacag$ L=2, l=3 3 3 acag$ 4 5 R=4, r=1 ag$ 5 2 cacag$ mlr = min(l,r)=1 6 4 cag$ M=(L+R)/2=3, m=4 7 6 g$ The pattern Q is found at SA[M]=3.
Time analysis Binary search will perform log n comparisons Each comparison takes at most O(m) time In the worst case, O(m log n) time. Myers and Manber report that, in practice, the time is O(m + log n).
Suffix array and suffix tree We show one example of replacing suffix tree by suffix array Note that most applications related to suffix tree can be solved using suffix array with some time blow up! When space is limited, replacing suffix tree by suffix array is a good choice.
The size is still too big! Why? DNA sequences can be very long! E.g. Fly: ~100M bases, Human: ~3G bases, Tree: ~9G bases Storage to store indexing data structure for human genome Suffix Tree: ~40G bytes Suffix Array: ~13G bytes Can we further reduce the space?
Solution Grossi, Vitter (STOC2000) Compressed suffix array (CSA) Ferragine, Manzini (FOCS2000) FM-index Both of them can be stored in O(n) bit space For Human Genome Both CSA and FM-index can be stored within 2G bytes.
FM-index Consider a text T=acacag$ FM-index stores: SA[i] Suffix T[SA[i]-1] The string BW=g$ccaaa A. 7 $ g C[x] = total no. of 1 acacag$ $ B. occurrences of each symbol 3 acag$ c less than x. E.g. C[a]=1, 5 ag$ c C[c]=4, C[g]=6, C[t]=7 2 cacag$ a A data-structure occ(x, i) C. 4 cag$ a which tells us the number of 6 g$ a occurrences of x in BW[1..i] using O(1) time.
Data-structure for answering the occ(x, i) query? n/log 2 n 1 2 3 4 BW ……… log 2 n log n ……… 1 2 3 log n 00x0xxx0x0 Given the text BW[1..n], we divide BW[1..n] into buckets of size log 2 n. For each bucket i = 1, … , n/log 2 n, we store P[i] = number of x ’ s in BW[1.. ilog 2 n]. Each bucket is further subdivided into log n sub-buckets of size log n. For each sub-bucket j of the bucket i, we store Q[i][j] = number of x ’ s in the first j sub-buckets.
Data-structure for answering the occ(x, i) query? We also need a lookup table rank(b,k) b is any string of length (log n)/2 1 ≤ k ≤ (log n)/2 rank(b,k) = number of x in the first k characters of b. 1 2 3 ……………… (log n)/2 Number of entries 00 … 0 log n log n log n = = 2 n 2 2 2 … 1 1 2 ……………………… x0x..0 Each entry take s O (log log n ) bits … Total space xx … x = = O ( n log n log log n ) o ( n ) bits
Space complexity of the occ() data-structure P[1..n/log 2 n] uses O(n/log n) bits Q[1..n/log 2 n][1..log n] uses O(n log log n / log n) bits rank(b,k) uses 2 logn/2 (log n/2) = o(n) bits In total, we use O(n log log n / log n) bits.
How to compute occ(x,i)? n/log 2 n 1 2 3 4 BW ……… log 2 n log n ……… 1 2 3 log n 00x0xxx0x0 Suppose log n = 10. To compute occ(x, 327), The result is P[3]+Q[4][2]+rank(00x0x,5)+rank(xx0x0,2) Hence, O(1) time to compute occ(x,i).
Size of FM-index Structure A can be store in 2n bits Structure B can be store in O(log n) bits Structure C can be store in O(n log log n/log n) bits In total, the size of FM-index is O(n) bits.
Observation C[x]+occ(x,i) is the number of SA[i] Suffix T[SA[i]-1] suffixes smaller than of xT SA[i] . 7 $ g 1 acacag$ $ Example: 3 acag$ c C[c]=4 5 ag$ c occ(c, 6)=2 2 cacag$ a Number of suffixes smaller than 4 cag$ a cT[SA[6]..n]=ccag$ is 6. 6 g$ a
Lemma Suppose range(T,Q) is [st..ed]. Then, range(T,xQ) = [p..q] where p = C[x] + occ(x,st-1) + 1 q = C[x] + occ(x,ed) Proof: p = 1+number of suffixes strictly smaller than xQ. The latter term = number of suffixes smaller than or equal to xT SA[st-1] . q = number of suffixes smaller than or equal to xT SA[ed] .
Backward search Given the text T and the FM-index, we want to determine if Q exists in T. Algorithm BW_exist(Q[1..m]) x=Q[m], i=m-1; 1. /* find range(T,Q[m]) */ 2. st = C[x]+1, ed = C[x+1]; while (st ≤ ed and i>1) { 3. /* find range(T, Q[i-1..m]) */ x = Q[i-1]; st = C[x] + occ(x, st-1) + 1; ed = C[x] + occ(x, ed); i = i – 1; } if st > ed, then pattern not found else pattern found. 4.
Example T = acacag$ SA[i] Suffix T[SA[i]-1] 7 $ g Q = aca 1 acacag$ $ 3 acag$ c 5 ag$ c Q[3..3]=a 2 cacag$ a sp=C[a] +1 4 cag$ a =1+1=2 6 g$ a ep=C[c] =4 g $ c c a a a
Example T = acacag$ SA[i] Suffix T[SA[i]-1] 7 $ g Q = aca 1 acacag$ $ 3 acag$ c 5 ag$ c Q[2..3]=ca 2 cacag$ a st=C[c]+occ(c,st old -1)+1 4 cag$ a =4+0+1=5 6 g$ a ed=C[c]+occ(c,ed old ) =4+2=6 g $ c c a a a
Example T = acacag$ SA[i] Suffix T[SA[i]-1] 7 $ g Q = aca 1 acacag$ $ 3 acag$ c 5 ag$ c Q[1..3]=aca 2 cacag$ a st=C[a]+occ(a,st old -1)+1 4 cag$ a =1+0+1=2 6 g$ a ed=C[a]+occ(a,ed old ) =1+2=3 g $ c c a a a
Example T = acacag$ SA[i] Suffix T[SA[i]-1] 7 $ g Q = aca 1 acacag$ $ 3 acag$ c 5 ag$ c Q[1..3]=aca 2 cacag$ a st=C[a]+occ(a,st old -1)+1 4 cag$ a =1+0+1=2 Q occurrences in T! 6 g$ a ed=C[a]+occ(a,ed old ) =1+2=3 g $ c c a a a
Time complexity of Backward Search To find a pattern Q[1..m] Step 1, 2, and 4 can be computed in O(1) time. For step 3, We need to iterate the loop for m-1 times. Each iteration of the loop can be computed in O(1) time. The loop takes O(m) time. In total, O(m) time for backward search.
Conclusion Suffix tree is a powerful data-structure which has a lot of applications in Computational Biology. Problems: Suffix tree is too big! Can be solved using CSA and FM-index Suffix tree can only solve exact match problem! (Most of the biology problems are approximate match!) Many works have been done on this area! But still not pratical. One of the important area to explore!
Recommend
More recommend
