Finding Optimal Pairs of Cooperative andCompeting Patterns with Bounded DistanceShunsuke Inenaga1 , Hideo Bannai2 , Heikki Hyyrö3 , Ayumi Shinohara3,5,
Masayuki Takeda4,5, Kenta Nakai2, and Satoru Miyano2
1

Department of Computer Science, P.O. Box 26 (Teollisuuskatu 23),
FIN-00014 University of Helsinki, Finland
inenaga@cs.helsinki.fi
2
Human Genome Center, Institute of Medical Science,
The University of Tokyo, Tokyo 108-8639, Japan{bannai,knakai,miyano}@ims.u-tokyo.ac.jp
3
PRESTO, Japan Science and Technology Agency (JST)helmu@cs.uta.fi
4
SORST, Japan Science and Technology Agency (JST)5Department of Informatics, Kyushu University 33, Fukuoka 812-8581, Japan{ayumi,takeda}@i.kyushu-u.ac.jp

Abstract. We consider the problem of discovering the optimal pair of
substring patterns with bounded distance α, from a given set S of strings.
We study two kinds of pattern classes, one is in form p ∧α q that areinterpreted as cooperative patterns within α distance, and the other isin form p ∧α ¬q representing competing patterns, with respect to S. We
show an eﬃcient algorithm to ﬁnd the optimal pair of patterns in O(N 2 )time using O(N ) space. We also present an O(m2 N 2 ) time and O(m2 N )space solution to a more diﬃcult version of the optimal pattern pairdiscovery problem, where m denotes the number of strings in S.

1

Introduction

Pattern discovery is an intensively studied sub-area of Discovery Science. A large
amount of eﬀort was paid to devising eﬃcient algorithms to extract interesting,useful, and even surprising substring patterns from massive string datasets suchas biological sequences [1, 2]. Then this research has been extended to more
complicated but very expressive pattern classes such as subsequence patterns [3,4], episode patterns [5, 6], VLDC patterns [7] and their variations [8].
Another interesting and challenging direction of this research is discoveryof optimal pattern pairs, whereas the above-mentioned algorithms are only forﬁnding optimal single patterns. Very recently, in [9] we developed an eﬃcientO(N 2 ) time algorithm for ﬁnding optimal pairs of substring patterns combinedwith any Boolean functions such as ∧ (AND), ∨ (OR) and ¬ (NOT), where N
denotes the total string length in the input dataset. For instance, the algorithm
allows to ﬁnd pattern pairs in form p ∧ q, which can be interpreted as twosequences with cooperative functions. Some developments have been made forﬁnding cooperative pattern pairs with a certain distance α, denoted, for now, byp ∧α q, in terms of structured motifs [10, 11] and proximity patterns [12, 13].
This paper extends special cases in the work of [9] by producing a genericalgorithm to ﬁnd the optimal pairs of cooperative patterns in form p ∧α q, andcompeting patterns in form p ∧α ¬q. We develop a very eﬃcient algorithm thatruns in O(N 2 ) time using only O(N ) space for any given threshold parameter α ≥0, hence not increasing the asymptotic complexity of [9] which does not deal withany bounded distance between two patterns. The pattern pairs discovered by
our algorithm are optimal in that they are guaranteed to be the highest scoringpair of patterns with respect to a given scoring function. Our algorithm can be
adjusted to handle several common problem formulations of pattern discovery,
for example, pattern discovery from positive and negative sequence sets [1, 14,
8, 15], as well as the discovery of patterns that correlate with a given numeric
attribute (e.g. gene expression level) assigned to the sequences [16–18, 2, 19].
The eﬃciency of our algorithm comes from the uses of helpful data structures, suﬃx trees [20] and sparse suﬃx trees [21, 22]. We also present an eﬃcient
implementation that uses suﬃx arrays and lcp arrays [23] to simulate bottom-up
traversals on suﬃx trees and sparse suﬃx trees.

Comparison to related works. Marsan and Sagot [10] gave an algorithm
which is limited to ﬁnding pattern pairs in form p ∧α q. Although the timecomplexity is claimed to be O(N m), where m denotes the number of strings in
the dataset, a signiﬁcant diﬀerence is that the length of each pattern is ﬁxed
and is regarded as a constant, whereas our algorithm does not impose such a
restriction. Arimura et al. [12, 13] presented an algorithm to ﬁnd the optimal(α, d)-proximity pattern that can be interpreted into a sequence p1 ∧α · · · ∧α pd
of d patterns with bounded distance α. Although the expected running timeof their algorithm is O(αd−1 N (log N )d ), the worst case running time is stillO(αd N d+1 log N ). Since in our case d = 2, it turns out to be O(α2 N 3 log N )which is rather bad compared to our O(N 2 ) time complexity. In addition, it isquite easy to extend our algorithm to ﬁnding (α, d)-proximity pattern in O(N d )time using only O(N ) space. Since the algorithm of [13] consumes O(dN ) space,
our algorithm improves both time and space for extracting proximity patterns.
On the other hand, an algorithm to discover pattern pairs in form p ∧α ¬q was
given recently in [24], in terms of missing patterns. As the algorithm performs in
O(n2 ) time using O(n) space for a single input string of length n, our algorithm
generalizes it for a set of strings with the same asymptotic complexity.
None of the above previous algorithms considers the fact that if one patternof pair p ∧α q is very long, its ending position may be outside of the region
speciﬁed by α, or for pattern pair in form p ∧α ¬q, it may invade the α-region
with respect to some other occurrence of p. Since such pattern pairs seem tobe of less interest for some applications, we design a version of our algorithm
which permits us to ﬁnd pattern pair p ∧α q such that both p and q occur strictly
within the α-region of each occurrence of p, and pattern pair p ∧α ¬q such that

α

α

α

α

s
q1

q3

p

q3

q4

q2

p

q1

Fig. 1. From pattern p in string s, q1 begins within α-distance, q2 begins outside αdistance, q3 occurs within α-distance, and q4 occurs outside α-distance.

q is strictly outside the α-region of any occurrences of p. This version of ouralgorithm runs in O(m2 N 2 ) time using O(m2 N ) space.

2

Preliminaries

2.1

NotationLet Σ be a ﬁnite alphabet. An element of Σ ∗ is called a string. Strings x, y, andz are said to be a preﬁx, substring, and suﬃx of string w = xyz, respectively.
The length of a string w is denoted by |w|. The empty string is denoted by ε,that is, |ε| = 0. The i-th character of a string w is denoted by w[i] for 1 ≤ i ≤ n,
and the substring of a string w that begins at position i and ends at position j isdenoted by w[i : j] for 1 ≤ i ≤ j ≤ n. For any set S, let |S| denote the cardinalityof S. Let [i, j] denote the set of consecutive integers i, i + 1, . . . , j.
For strings p and s, let Beg(p, s) denote the set of all beginning positions of
p in s, e.g.,
Beg(p, s) = {i | s[i : i + |p| − 1] = p}.
Also, let C ov(p, s) denote the set of all positions in s covered by p, e.g.,
C ov(p, s) = ∪j∈Beg(p,s) [j, j + |p| − 1].
For strings p, s and threshold value α ≥ 0, let us deﬁne set D(p, s, α) and of
integers by
D(p, s, α) = (∪j∈Beg(p,s) [j − α, j + α]) ∩ [1, |s|].
That is, D(p, s, α) consists of the regions inside s that are closer than or equalto α positions from the beginning position of some occurrence of p. The shaded
regions illustrated in Fig. 1 are the regions in D(p, s, α).
If |B eg(q, s) ∩ D(p, s, α)| ≥ 1, that is, if one or more occurrences of q in s
are within α positions from at least one occurrence of p in s, then q is said tobegin within α-distance from p, and otherwise begin outside α-distance. Similarly,
if |C ov(q, s) ∩ D(p, s, α)| ≥ 1, that is, if there is at least one overlap between
an occurrence of p and a region in D(p, s, α), then q is said to occur within
α-distance from p, and otherwise occur outside α-distance. See also Fig. 1.
Let ψ(p, s) be a Boolean matching function that has the value true if p isa substring of s, and false otherwise. We deﬁne p, F, q as a Boolean patternpair (or simply pattern pair) which consists of two patterns p, q and a Booleanfunction F : {true, false} × {true, false} → {true, false}. We say that asingle pattern or pattern pair π matches string s if and only if ψ(π, s) = true.
For α-distance, we deﬁne four Boolean matching functions F1 , F2 , F3 , and F4
as follows: F1 is such that ψ(p, F1 , q, s) = true if ψ(p, s) = true and q begins
within α-distance from p, and ψ(p, F1 , q, s) = false otherwise. F2 is such thatψ(p, F2 , q, s) = true if ψ(p, s) = true and q begins outside α-distance fromp, and ψ(p, F2 , q, s) = false otherwise. F3 is such that ψ(p, F3 , q, s) = trueif ψ(p, s) = true and q occurs within α-distance from p, and ψ(p, F3 , q, s) =false otherwise. F4 is such that ψ(p, F4 , q, s) = true if ψ(p, s) = true and q
occurs outside α-distance from p, and ψ(p, F4 , q, s) = false otherwise. In thesequel, let us write as:
p, F1 , q = p ∧b(α) q, p, F2, q = p ∧b(α) ¬q,
p, F3, q = p ∧c(α) q, p, F4 , q = p ∧c(α) ¬q.
Note that, by setting α ≥ |s|, we always have D(p, s, α) = [1, |s|]. Therefore,
the above pattern pairs are generalizations of pattern pairs such as (p ∧ q) and(p ∧ ¬q) that are considered in our previous work [9].
Given a set S = {s1 , . . . , sm } of strings, let M (π, S) denote the subset of
strings in S that π matches, that is,
M (π, S) = {si ∈ S | ψ(π, si ) = true}.
We suppose that each si ∈ S is associated with a numeric attribute value ri . Fora single pattern or pattern pair π, let M (π,S) ri denote the sum of ri over all
si in S such that ψ(π, si ) = true, that is,ri =(ri | ψ(π, si ) = true).
M (π,S)si ∈SAs Sis usually ﬁxed in our case, we shall omit
S where possible and let M (π)
and M (π) ri be a shorthand for M (π, S) and M (π,S) ri , respectively.
2.2

Problem Definition

In general, the problem of ﬁnding a good pattern from a given set S of strings
refers to ﬁnding a pattern π that maximizes some suitable scoring function scorewith respect to the 
strings in S. We concentrate on score that takes parametersof type |M (π)| and M (π) ri , and assume that the score value computation itself
can be done in constant time if the required parameter values are known. We
formulate the pattern pair discovery problem as follows:
Problem 1. Given a set S = {s1 , . . . , sm } of m strings, where each string si
is assigned a numeric attribute value ri , a threshold parameter α ≥ 0, a scoring function score : R × R → R, and a Boolean Function F ∈ {F1 , . . . , F4 },∗ﬁnd the Boolean pattern pair π ∈ {p, F, q | p, q ∈ Σ } that maximizes
score(|M (π)|, M (π) ri ).

The speciﬁc choice of the scoring function depends highly on the particularapplication, for example, as follows: The positive/negative sequence set discrimination problem [1, 14, 8, 15] is such that, given two disjoint sets of sequences S1and S2 , where sequences in S1 (the positive set) are known to have some biological function, while the sequences in S2 (the negative set) are known not to, ﬁnd
pattern pairs which match more sequences in S1 , and less in S2 . Common scoring
functions that are used in this situation include the entropy information gain,
the Gini index, and the chi-square statistic, which all are essentially functions of|M (π, S1 )|, |M (π, S2 )|, |S1 | and |S2 |. The correlated patterns problem [16–18, 2,
19] is such that we are given a set S of sequences, with a numeric attribute value
ri associated with each sequence si ∈ S, and the task is to ﬁnd pattern pairs
whose occurrences in the sequences correlate with their numeric attributes. In
this framework,
scoring functions such as mean squared error, which is a functionof |M (π)| and M (π) ri , can be used.
2.3

Data Structures

We intensively use the following data structures in our algorithms. The eﬃciencyof our algorithms, in both time and space consumptions, comes from the use of
these data structures.
The suﬃx tree [20] for a string s is a rooted tree whose edges are labeledwith substrings of s, satisfying the following characteristics. For any node v inthe suﬃx tree, let l(v) denote the string spelled out by concatenating the edge
labels on the path from the root to v. For each leaf v, l(v) is a distinct suﬃx
of s, and for each suﬃx in s, there exists such a leaf v. Each leaf v is labeledwith pos(v) = i such that l(v) = s[i : |s|]. Also, each internal node has at least
two children, and the ﬁrst character of the labels on the edges to its children are
distinct. It is well known that suﬃx trees can be represented in linear space andconstructed in linear time with respect to the length of the string [20].
The generalized suﬃx tree (GST) for a set S = {s1 , . . . , sm } of strings isbasically the suﬃx tree for the string s1 $1 · · · sm $m , where each $i (1 ≤ i ≤ m)
is a distinct character which does not appear in any of the strings in the set.
However, all paths are ended at the ﬁrst appearance of any $i , and each leaf islabeled with idi that speciﬁes the string in S the leaf is associated with. For set{abab, aabb}, the corresponding GST is shown in Fig. 2. GSTs are also linear
time constructible and can be stored in linear space with respect to the total
length of the strings in the input set [20].
The sparse suﬃx tree (SST) [21, 22] for a string s is a rooted tree which
represents a subset of the suﬃxes of s, i.e., a suﬃx tree of s from which the deadleaves (and corresponding internal nodes) are removed, where the dead leaves
stand for the leaves that correspond to the suﬃxes not in the subset. On theother hand, the leaves existing in the SST are called the living leaves. An SST
can be constructed in linear time by once building the full suﬃx tree and pruning
dead leaves from it. (Direct construction of SSTs was also considered in [21, 22]
but the pruning approach is suﬃcient for our purposes.)

As’ lcps’

a
b
b a
b
$1 b
$2
$1

a
b
b
$2

6

1

a
b
$1 b
$2
7

$1

3

2

$2

$2
$1

8

4

9

5

10

6
1
7
3
2
8
4
9
5
10

0
1
2
2
0
1
1
1
0
0

Fig. 2. The tree is the GST for set {abab, aabb}. Each leaf y is labeled by pos(y) withrespect to the concatenated string s = abab$1 aabb$2 . The left and right arrays are
the suﬃx and lcp arrays of string s , respectively.

The suﬃx array [23] As for a string s of length n, is a permutation of theintegers 1, . . . , n representing the lexicographic ordering of the suﬃxes of s.
The value As [i] = j in the array indicates that s[j : n] is the ith suﬃx in thelexicographic ordering. The lcp array for a string s, denoted lcps , is an array ofintegers representing the longest common preﬁx lengths of adjacent suﬃxes inthe suﬃx array, that is,
lcps [i] = max{k | s[As [i − 1] : As[i − 1] + k − 1] = s[As [i] : As[i] + k − 1]}.
Note that each position i in As corresponds to a leaf in the suﬃx tree of s, andlcps [i] denotes the length of the path from the root to leaves of positions i − 1
and i for which the labels are equivalent. Fig. 2 shows the suﬃx array As andthe lcp array lcps of string s = abab$1 aabb$2 .
Recently, three methods for constructing the suﬃx array directly from astring in linear time have been developed [25–27]. The lcp array can be constructed from the suﬃx array also in linear time [28]. It has been shown thatseveral (and potentially many more) algorithms which utilize the suﬃx tree
can be implemented very eﬃciently using the suﬃx array together with its lcp
array [28, 29]. This paper presents yet another good example for eﬃcient implementation of an algorithm based conceptually on full and sparse suﬃx trees, butuses the suﬃx and lcp arrays.
The lowest common ancestor lca(x, y) of any two nodes x and y in a tree isthe deepest node which is common to the paths from the root to both of thenodes. The tree can be pre-processed in linear time to answer the lowest common
ancestor (lca-query) for any given pair of nodes in constant time [30]. In termsof the suﬃx array, the lca-query is almost equivalent to a range minimum query(rm-query) on the lcp array, which, given a pair of positions i and j, rmq(i, j)
returns the position of the minimum element in the sub-array lcp[i : j]. The lcp
array can also be pre-processed in linear time to answer the rm-query in constanttime [30, 31].

In the sequel, we will show a clever implementation of our algorithms which
simulates bottom-up traversals on SGSTs using suﬃx and lcp arrays. Our algorithms need no explicit prunings of leaves from the SGSTs which means thatthe suﬃx and lcp arrays do not have to be explicitly modiﬁed or reconstructed.
Therefore our algorithms are both time and space eﬃcient in practice.

3

Algorithms

In this section we present our eﬃcient algorithms to solve Problem 1 for Booleanfunctions F1 , . . . , F4. Let N = |s1 $1 · · · sm $m |, where S = {s1 , . . . , sm }. The
algorithms calculate scores for all possible pattern pairs, and output a patternpair of the highest score.
3.1

Algorithm for Single Patterns

For a single pattern, we need to consider only patterns of form l(v), where v
is a node in the GST for S: If a pattern has a corresponding path that ends
in the middle of an edge of the GST, it will match the same set of strings as
the pattern corresponding to the end node of that same edge, and hence thescore would be the same. Below, we recall the O(N ) time algorithm of [9] thatcomputes M (l(v)) ri for all single pattern candidates of form l(v). The algorithm
is derived from the technique for solving the color set size problem [32]. Notethat we do not need to consider separately
how to compute |M (l(v))|: If we giveeach attribute ri the value 1, then M (l(v)) ri = |M (l(v))|.
For any node v in the GST of S, let LF (v) be the set of all leaves in thesubtree rooted by v, ci (v) be the number of leaves in LF (v) with the label idi ,
and


ri =
(ci (v)ri | ψ(l(v), si ) = true).
LF (v)For any such v, ψ(l(v), si ) = true for at least one string si . Thus, we haveri =(ri | ψ(l(v), si ) = true)
M (l(v))

=



ri −



((ci (v) − 1)ri | ψ(l(v), si ) = true).

LF (v)We deﬁne the 
above subtracted sum to be a correction factor, denoted bycorr(l(v), S) = ((ci (v) − 1)ri | ψ(l(v), si ) = true). The following recurrence

 
ri =
(
ri | v is a child of v)
LF (v)

LF (v  )allows us to compute the values LF (v) ri for all v during a linear time bottomup traversal of the GST.
 Now we need to remove the redundancies, represented
by corr(l(v), S), from LF (v) ri .

Let I(idi ) be the list of all leaves with the label idi in the order they appear
in a depth-ﬁrst traversal of the tree. Clearly the lists I(idi ) can be constructed inlinear time for all labels idi . We note four properties in the following proposition:
Proposition 1. For the GST of S, the following properties hold:(1) The leaves in LF (v) with the label idi form a continuous interval of lengthci (v) in the list I(idi ).
(2) If ci (v) > 0, a length-ci (v) interval in I(idi ) contains ci (v) − 1 adjacent(overlapping) leaf pairs.
(3) If x, y ∈ LF (v), the node lca(x, y) belongs to the subtree rooted by v.
(4) For any si ∈ S, ψ(l(v), si ) = true if and only if there is a leaf x ∈ LF (v)with the label idi .
We initialize each node v to have a correction value 0. Then, for each adjacentleaf pair x, y in the list I(idi ), we add the value ri into the correction value ofthe node lca(x, y). It follows from properties (1) - (3) of Proposition 1, that now
the sum of the correction values in the nodes of the subtree rooted by v equals(ci (v) − 1)ri . After repeating the process for each of the lists I(idi ), property (4)tells that the preceding total sum of the correction values in the subtree rooted
by v becomes((ci (v) − 1)ri | ψ(l(v), si ) = true) = corr(l(v), S).

Now a single linear time bottom-up traversal of the tree enables us to cumulate
the correction values corr(l(v), S) from the subtrees
into each node v, and atthe same time we may record the ﬁnal values M (l(v)) ri .
The preceding process involves a constant number of linear time traversalsof the tree, as well as a linear number ofconstant time lca-queries (after a linear
time preprocessing). Hence the values M (l(v)) ri are computed in linear time.
3.2

Algorithms for Pattern Pairs with α-Distance
The algorithm described in Section 3.1 permits us to computer inM (l(v))
 i
O(N ) time. In this section, we concentrate on how to compute the values M (π) ri
for various types of pattern pair π with bounded distance α. In so doing, we goover the O(N ) choices for the ﬁrst pattern, and for each ﬁxed pattern l(v1 ), we
compute M (π) ri where π consists of l(v1 ) and second pattern candidate l(v).
The result of this section is summarized in the following theorem:
Theorem 1. Problem 1 is solvable in O(N 2 ) time using O(N ) space for F ∈{F1 , F2}, and O(m2 N 2 ) time and O(m2 N ) space for F ∈ {F3 , F4}.
Below, we will give the details of our algorithm for each type of pattern pair.

For Pattern Pairs in form p∧b(α) q. Let D(l(v1 ), S, α) = ∪si ∈S D(l(v1 ), si , α).
It is not diﬃcult to see that D(l(v1 ), S, α) is O(N ) time constructible (see Theorem 9 of [24]). We then mark as ‘living’ every leaf of the GST whose position belongs to D(l(v1 ), S, α). This aims at conceptually constructing the sparse generalized suﬃx tree (SGST) for the subset of suﬃxes corresponding to D(l(v1 ), S, α),and this can be done in O(N ) time. Then, we additionally label each stringsi ∈ S, and the corresponding leaves in the SGST, with the Boolean value
ψ(l(v1 ), si ). This can be done in O(N ) time. Now the trick is to cumulate thesums and correction factors only for the nodes existing in the SGST, and separately for diﬀerent values of theadditional label. The end result is that we willobtain the values M (π) ri = (ri | ψ(π, si ) = true) in linear time, where πdenotes the pattern pair l(v1 ) ∧b(α) l(v) for all nodes v in the conceptual SGST.
Since there are only O(N ) candidates for l(v1 ), the total time complexity isO(N 2 ). The space requirement is O(N ), since we repeatedly use the same GST
for S and the additional information storage is also linear in N at each phase of
the algorithm.
For Pattern Pairs in form p ∧b(α) ¬q. Let us denote any pattern pair inform p ∧b(α) ¬q by π. For the pattern pairs of this form, it stands thatri =(ri | ψ(π, si ) = true)
M (π)

=



(ri | ψ(p, s) = true) −



ri ,

M (π)



where π denotes p ∧b(α) q. For each ﬁxed ﬁrst pattern l(v1 ), M (π) ri and (ri |ψ(l(v1 ), s) = true) can be computed in O(N ) time, and hence M (π) ri canbe computed in O(N ) time as well. Since we have O(N ) choices for the ﬁrst
pattern, the total time complexity is O(N 2 ), and the space complexity is O(N ).
For Pattern Pairs in form p∧c(α) q. There are two following types of patterns
that have to be considered as candidates for the second pattern q:(1) pattern q which has an occurrence beginning at some position in D(l(v1 ), S, α).
(2) pattern q which has an occurrence beginning at some position in [1, N ] −D(l(v1 ), S, α), but overlapping with some positions in D(l(v1 ), S, α).
Now, we separately treat these two kinds of patterns. In so doing, we consider
two SGSTs, one consists only of the leaves x such that pos(x) ∈ D(l(v1 ), S, α),and the other consists only of the leaves y such that pos(y) ∈ [1, N ]−D(l(v1 ), S, α).
Denote these by SGST1 and SGST2 , respectively.
For integer i with 1 ≤ i ≤ m, let Ni = |s1 $1 · · · si $i |. For any leaf y of SGST2
with idi , we deﬁne l en(y) as follows:
Ni−1 + h − α − pos(y) if pos(y) < Ni−1 + k − α,l en(y) =Ni − pos(y)otherwise,

α

Ni-1
S1

$1

..........

αlen(y)ααlen(y’)

Si

$i
l(v1)

..........

$m

l(v1)
pos(y)pos(y’)

Fig. 3. For two leaves y, y of SGST2 , l en(y) and l en(y  ) with respect to string si ∈ S
and the ﬁrst pattern l(v1 ).

where h denotes the minimum element of Beg(l(v1 ), si ) satisfying Ni−1 + h >pos(y), and k = max(B eg(l(v1 ), si )). Fig. 3 illustrates l en(y) with respect tosi and l(v1 ). Computing l en(y) can be done in constant time for each leaf y,
after a linear time preprocessing of scanning D(l(v1 ), S, α) from right to left.
For any node v of SGST2 , let S LFi (v) denote the list of all leaves that are in
the subtree of SGST2 rooted at v and are labeled with idi . Then, we computel eni (v) = miny∈SLFi (v)(l en(y)) for every internal node v in SGST2 . For all nodesv and i = 1, . . . , m, l eni (v) can be computed in O(mN) time using O(mN)
space. Notice that only nodes v such that for some i, l eni (v) < |l(v)| in SGST2 ,can give candidates for the second pattern, since the occurrences of the substrings
represented by the other nodes are completely outside D(l(v1 ), S, α).
Consider an internal node v and its parent node v in SGST2 . Assume thelength of the edge between v and v is more than one. Suppose the subtree of
SGST2 rooted at v contains a leaf with idi . Now consider an ‘implicit’ nodet that is on the edge from v to v . Then, even if l eni (v ) < |l(v )|, it might
stand that l eni (t) ≥ |l(t)| since l eni (t) = l eni (v ) and |l(t)| < |l(v )| always
hold. Note that in this case we have to take into account such implicit node t
for the second pattern candidate, independently from v . On the other hand, if
we locate t so that l eni (t) = |l(t)|, we will have l eni (t ) > |l(t )| for any implicit
node t upper than t in the edge, and will have l eni (t ) < |l(t )| for any implicit
node t lower than t in the edge. Thus this single choice of t suﬃces. Since thesubtree can contain at most m distinct id’s, each edge can have at most m such
implicit nodes on it, and because SGST2 has O(N ) edges, the total numberof the implicit nodes we have to treat becomes O(mN). What still remains is
how to locate these implicit nodes in O(mN) time. This is actually feasible byusing the so-called suﬃx links and the canonization technique introduced in [33].
For each of these implicit nodes we have to maintain m number of information
as mentioned in the above paragraph, and hence in total we are here requiredO(m2 N ) time and space.
Cumulating the sums and correction factors is done in two rounds: For simplicity, let us denote pattern pair l(v1 ) ∧c(α) l(v) by π. First we traverse SGST1
and for each encountered node v, calculate the value
M (π),1 ri and store it
in v. Then,wetraverse
SGST
and
for
each
encounterednode w calculate2the value M (π),2 ri if necessary, while checking whether l eni (w) < |l(w)| ornot, and store it in w but separately from
M (π),1 ri . Then, the ﬁnal value



r
=
r
+
r
.
Since
we
go
over O(N ) choices for the ﬁrst
M (π) iM (π),1 i
M (π),2 i
2 2pattern, the resulting algorithm requires O(m N ) time and O(m2 N ) space.
For Pattern Pairs in form p ∧c(α) ¬q. As similar arguments to pattern pairs
in form p ∧b(α) ¬q hold, it can also be done O(m2 N 2 ) time and O(m2 N ) space.
3.3

Implementation

We show an eﬃcient implementation of the algorithms for pattern pairs in formp ∧b(α) q and p ∧b(α) ¬q, both run in O(N 2 ) time using O(N ) space as previously
explained. We further modify the implementation of the algorithm in [9], which
is an extension of the Substring Statistics algorithm in [28], that uses the suﬃx
and lcp arrays to simulate the bottom up traversals of the generalized suﬃx tree.
We note that the simulation does not increase the asymptotic complexity of thealgorithm, but rather it is expected to increase the eﬃciency of the traversal, as
conﬁrmed in previous works [34, 28, 29].
Fig. 4 shows a pseudo-code for solving the Color Set Size problem with pruned
edges using a suﬃx array, which is used in the algorithm of the previous subsection. There are two diﬀerences from [9]: (1) the correction factors are set tolca(i, j) of consecutive ‘living’ leaves i, j (i < j) corresponding to the same string.
Since the simulation via suﬃx arrays does not explicitly construct all internal
nodes of the suﬃx tree, we use another array CF to hold the correction factors.
The correction factor for these leaves i, j is summed into CF [rmq(i + 1, j)]. (2)‘Dead’ leaves are ignored in the bottom-up traversal. This can be done in several
ways, the simplest which assigns 0 to the weights of each ‘dead’ leaf. However,
although the worst-case time complexity remains the same, such an algorithm
would report redundant internal nodes that only have at most a single ‘living’ leafin its subtree. A more complete pruning can be achieved, as shown in the pseudocode, in the following way: for any consecutive ‘living’ leaves i, j (i < j) in thesuﬃx array, we add up all CF [k] for i < k ≤ j which is then added to the nodecorresponding to lca(i, j). The correctness of this algorithm can be conﬁrmedby the following argument: For a given correction factor CF [k] (i < k ≤ j),since the query result came from consecutive ‘living’ leaves of the same string,
they must have been from leaves i , j  , such that i ≤ i and j ≤ j  . Thismeans that lcp[k] = lcp[rmq(i + 1, j  )] ≤ lcp[rmq(i + 1, j)] ≤ lcp[k]. Therefore
lcp[rmq(i + 1, j  )] = lcp[rmq(i + 1, j)] = lcp[k], and lca(i , j  ) = lca(i, j), whichshows that CF [k] should be summed into the node with depth lcp[rmq(i + 1, j)].

4Computational Experiments

We tested the eﬀectiveness of our algorithm using the two sets of predicted
3’UTR processing site sequences provided in [35], which are constructed basedon the microarray experiments in [36] that measure the degradation rate of yeast

1 Let Stack = {(0, −1, 0)} be the stack;2 cc := 0; iplus1 = 0;
3 foreach j = 1, . . . , n + 1 do:
4
cc := cc + CF [j];5
(L, H, C) = top(Stack);
6
if j is a living leaf
7
(Lj , Hj , Cj ) := (j, lcp[rmq(iplus1, j)], 0);8while (H > Hj ) do:
9(L, H, C) := pop(Stack);
10
report (L, H, C); /* s[A[L] : A[L] + H − 1] has count C */11
Cj := C + Cj ;12
(L, H, C) := top(Stack);
13
if (H < Hj ) then14push((Lj , Hj , Cj + cc), Stack);
15
else /* H = Hj */16(L, H) := pop(Stack);17
push((L, H, C + Cj + cc), Stack);18
push((j, N − A[j] + 1, 1), Stack);
19
cc := 0; iplus1 = j + 1;
Fig. 4. Core of the algorithm for solving the Color Set Size problem with pruned edges
using a suﬃx array. We assume each leaf (position in array) has been labeled ‘dead’or ‘living’, and the correction factors for ‘living’ leaves are stored in the array CF . A
node in the suﬃx tree together with |M (l(v))| is represented by a three-tuple (L,H,C).

mRNA. One set Sf consists of 393 sequences which have a fast degradation rate
(t1/2 < 10 minutes), while the other set Ss consists of 379 predicted 3’UTR
processing site sequences which have a slow degradation rate (t1/2 > 50 minutes).
Each sequence is 100 nt long, and the total length of the sequences is 77, 200
nt. For the scoring function, we used the chi-squared statistic, calculated by(|Sf |+|Ss |)(tp∗tn−fp∗fn)2 /(tp + fn)(tp + fp)(tn + fp)(tn + fn) where tp =|M (π, Sf )|, fp = |Sf | − tp, tn = |Ss | − fn, and fn = |M (π, Ss )|.
For b(α) and 1 ≤ α ≤ 30, the best score was 48.3 (p < 10−11 ) for the pair
AUA∧b(10) ¬UGUA, which matched 159/393 and 248/379 fast and slowly degrading
sequences respectively, meaning that this pattern pair is more frequent in theslowly degrading sequences. This result is similar to the best pattern of the formp ∧ ¬q obtained in [9], which was UGUA ∧ ¬AUCC with score 33.9 (240/393 and152/379 (p < 10−8 )), appearing more frequently in the fast degrading sequences.
In fact, UGUA is known as a binding site of the PUF protein family which plays
important roles in mRNA regulation [37]. Of the 268/393 and 190/379 sequencesthat contain both AUA and UGUA, their distances were farther than α = 10 in only
37/268 of fast degrading sequences, while the number was 67/190 for slowly
degrading sequences. This could mean that AUA is another sequence element
whose distance from UGUA inﬂuences how eﬃciently UGUA functions.
We did not observe any notable diﬀerence in computational time comparedto the algorithm in [9], that the marking of each leaf as ‘dead’ or ‘living’ could

α

α

α

α

s
q4

q1

q3

p

q4

q2

q2

p

q3

q1

Fig. 5. From pattern p in string s, q1 begins within α-gap, q2 begins outside α-gap, q3
occurs within α-gap, and q4 occurs outside α-gap.

cause. For a given value of α, computation took around 2140 seconds for theabove dataset on a PC with Xeon 3 GHz Processor running Linux.

5

Concluding Remarks

In this paper we studied the problem of ﬁnding the optimal pair of substring
patterns p, q with bounded distance α, from a given set S of strings, which isan extension of the problem studied in [9]. We developed an eﬃcient algorithm
that ﬁnds the best pair of cooperative patterns in form p ∧b(α) q, and competing
patterns in form p ∧b(α) ¬q, which performs in O(N 2 ) time using O(N ) space,
where N is the total length of the strings in S. For the more diﬃcult versions of
the problem, referred to ﬁnding p ∧c(α) q and p ∧c(α) ¬q, we gave an algorithm
running in O(m2 N 2 ) time and O(m2 N ) space, where m = |S|. An interesting
open problem is if the m-factors can be removed from the complexities.
Any pairs output from our algorithms are guaranteed to be optimal in thepattern classes, in the sense that they give the highest score due to the scoring function. Our algorithms are adapted to various applications such as the
positive/negative sequence set discrimination problem for which the entropy information gain, the Gini index, and the chi-square statistic are commonly usedas the scoring function, and the correlated patterns problem for which the meansquared error can be used.
It is possible to use our algorithms for bounded gaps α (see Fig. 5). Since theending positions of the ﬁrst patterns have to be considered in this case, we will
naturally have O(N 2 ) choices for the ﬁrst pattern, and this fact suggests the useof suﬃx tries rather than suﬃx trees for the ﬁrst patterns. Still, there are only
O(N ) choices for the second pattern, and thus for the pattern pairs with respect
to b(α), we can modify our algorithm to solve the problems in O(N 3 ) time using
O(N 2 ) space. We remark that for c(α), a modiﬁed algorithm ﬁnds the optimal
pattern pair in O(m2 N 3 ) time using O(N 2 + m2 N ) space.

References
1. Shimozono, S., Shinohara, A., Shinohara, T., Miyano, S., Kuhara, S., Arikawa,
S.: Knowledge acquisition from amino acid sequences by machine learning system
BONSAI. Transactions of Information Processing Society of Japan 35 (1994)
2009–2018
2. Bannai, H., Inenaga, S., Shinohara, A., Takeda, M., Miyano, S.: Eﬃciently ﬁnding
regulatory elements using correlation with gene expression. Journal of Bioinformatics and Computational Biology (2004) in press.

3. Baeza-Yates, R.A.: Searching subsequences (note). Theoretical Computer Science78 (1991) 363–376
4. Hirao, M., Hoshino, H., Shinohara, A., Takeda, M., Arikawa, S.: A practical algorithm to ﬁnd the best subsequence patterns. In: Proc. 3rd International Conference
on Discovery Science (DS’00). Volume 1967 of LNAI., Springer-Verlag (2000) 141–
154
5. Mannila, H., Toivonen, H., Verkamo, A.I.: Discovering frequent episode in sequences. In: Proc. 1st International Conference on Knowledge Discovery and Data
Mining, AAAI Press (1995) 210–215
6. Hirao, M., Inenaga, S., Shinohara, A., Takeda, M., Arikawa, S.: A practical algorithm to ﬁnd the best episode patterns. In: Proc. 4th International Conference on
Discovery Science (DS’01). Volume 2226 of LNAI., Springer-Verlag (2001) 435–440
7. Inenaga, S., Bannai, H., Shinohara, A., Takeda, M., Arikawa, S.: Discovering
best variable-length-don’t-care patterns. In: Proc. 5th International Conference on
Discovery Science (DS’02). Volume 2534 of LNCS., Springer-Verlag (2002) 86–97
8. Takeda, M., Inenaga, S., Bannai, H., Shinohara, A., Arikawa, S.: Discovering mostclassiﬁcatory patterns for very expressive pattern classes. In: Proc. 6th International Conference on Discovery Science (DS’03). Volume 2843 of LNCS., SpringerVerlag (2003) 486–493
9. Bannai, H., Hyyrö, H., Shinohara, A., Takeda, M., Nakai, K., Miyano, S.: Finding
optimal pairs of patterns. In: Proc. 4th Workshop on Algorithms in Bioinformatics(WABI’04). (2004) to appear.
10. Marsan, L., Sagot, M.F.: Algorithms for extracting structured motifs using a suﬃx
tree with an application to promoter and regulatory site consensus identiﬁcation.
J. Comput. Biol. 7 (2000) 345–360
11. Palopoli, L., Terracina, G.: Discovering frequent structured patterns from string
databases: an application to biological sequences. In: 5th International Conference
on Discovery Science (DS’02). Volume 2534 of LNCS., Springer-Verlag (2002) 34–
46
12. Arimura, H., Arikawa, S., Shimozono, S.: Eﬃcient discovery of optimal wordassociation patterns in large text databases. New Generation Computing 18 (2000)
49–60
13. Arimura, H., Asaka, H., Sakamoto, H., Arikawa, S.: Eﬃcient discovery of proximity
patterns with suﬃx arrays (extended abstract). In: Proc. the 12th Annual Symposium on Combinatorial Pattern Matching (CPM’01). Volume 2089 of LNCS.,
Springer-Verlag (2001) 152–156
14. Shinohara, A., Takeda, M., Arikawa, S., Hirao, M., Hoshino, H., Inenaga, S.: Finding best patterns practically. In: Progress in Discovery Science. Volume 2281 of
LNAI., Springer-Verlag (2002) 307–317
15. Shinozaki, D., Akutsu, T., Maruyama, O.: Finding optimal degenerate patterns in
DNA sequences. Bioinformatics 19 (2003) 206ii–214ii
16. Bussemaker, H.J., Li, H., Siggia, E.D.: Regulatory element detection using correlation with expression. Nature Genetics 27 (2001) 167–171
17. Bannai, H., Inenaga, S., Shinohara, A., Takeda, M., Miyano, S.: A string pattern regression algorithm and its application to pattern discovery in long introns.
Genome Informatics 13 (2002) 3–11
18. Conlon, E.M., Liu, X.S., Lieb, J.D., Liu, J.S.: Integrating regulatory motif discovery and genome-wide expression analysis. Proc. Natl. Acad. Sci. 100 (2003)
3339–3344

19. Zilberstein, C.B.Z., Eskin, E., Yakhini, Z.: Using expression data to discover RNA
and DNA regulatory sequence motifs. In: The First Annual RECOMB Satellite
Workshop on Regulatory Genomics. (2004)
20. Gusﬁeld, D.: Algorithms on Strings, Trees, and Sequences. Cambridge University
Press (1997)
21. Andersson, A., Larsson, N.J., Swanson, K.: Suﬃx trees on words. Algorithmica
23 (1999) 246–260
22. Kärkkänen, J., Ukkonen, E.: Sparse suﬃx trees. In: Proceedings of the Second
International Computing and Combinatorics Conference (COCOON’96). Volume1090 of LNCS., Springer-Verlag (1996) 219–230
23. Manber, U., Myers, G.: Suﬃx arrays: a new method for on-line string searches.
SIAM Journal on Computing 22 (1993) 935–948
24. Inenaga, S., Kivioja, T., Mäkinen, V.: Finding missing patterns. In: Proc. 4thWorkshop on Algorithms in Bioinformatics (WABI’04). (2004) to appear.
25. Kim, D.K., Sim, J.S., Park, H., Park, K.: Linear-time construction of suﬃx arrays. In: 14th Annual Symposium on Combinatorial Pattern Matching (CPM’03).
Volume 2676 of LNCS., Springer-Verlag (2003) 186–199
26. Ko, P., Aluru, S.: Space eﬃcient linear time construction of suﬃx arrays. In: 14th
Annual Symposium on Combinatorial Pattern Matching (CPM’03). Volume 2676
of LNCS., Springer-Verlag (2003) 200–210
27. Kärkkäinen, J., Sanders, P.: Simple linear work suﬃx array construction. In: 30th
International Colloquium on Automata, Languages and Programming (ICALP’03).
Volume 2719 of LNCS., Springer-Verlag (2003) 943–955
28. Kasai, T., Arimura, H., Arikawa, S.: Eﬃcient substring traversal with suﬃx arrays.
Technical Report 185, Department of Informatics, Kyushu University (2001)
29. Abouelhoda, M.I., Kurtz, S., Ohlebusch, E.: The enhanced suﬃx array and its
applications to genome analysis. In: Second International Workshop on Algorithms
in Bioinformatics (WABI’02). Volume 2452 of LNCS., Springer-Verlag (2002) 449–
463
30. Bender, M.A., Farach-Colton, M.: The LCA problem revisited. In: Proc. 4th Latin
American Symp. Theoretical Informatics (LATIN’00). Volume 1776 of LNCS.,
Springer-Verlag (2000) 88–94
31. Alstrup, S., Gavoille, C., Kaplan, H., Rauhe, T.: Nearest common ancestors: asurvey and a new distributed algorithm. In: 14th annual ACM symposium on
Parallel algorithms and architectures (SPAA’02). (2002) 258–264
32. Hui, L.: Color set size problem with applications to string matching. In: Proc. 3rd
Annual Symposium on Combinatorial Pattern Matching (CPM’92). Volume 644
of LNCS., Springer-Verlag (1992) 230–243
33. Ukkonen, E.: On-line construction of suﬃx trees. Algorithmica 14 (1995) 249–260
34. Kasai, T., Lee, G., Arimura, H., Arikawa, S., Park, K.: Linear-time longestcommon-preﬁx computation in suﬃx arrays and its applications. In: 12th Annual Symposium on Combinatorial Pattern Matching (CPM’01). Volume 2089 of
LNCS., Springer-Verlag (2001) 181–192
35. Graber, J.: Variations in yeast 3’-processing cis-elements correlate with transcript
stability. Trends Genet. 19 (2003) 473–476 http://harlequin.jax.org/yeast/
turnover/.
36. Wang, Y., Liu, C., Storey, J., Tibshirani, R., Herschlag, D., Brown, P.: Precisionand functional speciﬁcity in mRNA decay. Proc. Natl. Acad. Sci. 99 (2002) 5860–
5865
37. Wickens, M., Bernstein, D.S., Kimble, J., Parker, R.: A PUF family portrait:
3’UTR regulation as a way of life. Trends Genet. 18 (2002) 150–157