Practical Algorithms for Pattern Based Linear
Regression
Hideo Bannai1 , Kohei Hatano1 , Shunsuke Inenaga1,2 , and Masayuki Takeda1,3
1

Department of Informatics, Kyushu University,
6-10-1 Hakozaki, Higashi-ku, Fukuoka 812-8581, Japan{bannai, hatano, shunsuke.inenaga, takeda}@i.kyushu-u.ac.jp
2
Japan Society for the Promotion of Science
3
SORST, Japan Science and Technology Agency (JST)Abstract. We consider the problem of discovering the optimal patternfrom a set of strings and associated numeric attribute values. The goodness of a pattern is measured by the correlation between the numberof occurrences of the pattern in each string, and the numeric attributevalue assigned to the string. We present two algorithms based on suﬃx
trees, that can ﬁnd the optimal substring pattern in O(N n) and O(N 2 )time, respectively, where n is the number of strings and N is their total
length. We further present a general branch and bound strategy that canbe used when considering more complex pattern classes. We also showthat combining the O(N 2 ) algorithm and the branch and bound heuristic
increases the eﬃciency of the algorithm considerably.

1

Introduction

Fundamental biological molecules such as DNA, RNA, and proteins can be regarded as strings over a certain alphabet. Although the whole genomic sequencesof many species are now becoming available, there is still much that is unknown
about the information that lie hidden in them. Computational analysis of thesesequences rely on the principle that similarity as strings implies similarity in theirsequence structure, which in turn implies similarity in their functions. Therefore, methods for eﬃciently and eﬀectively discovering meaningful patterns, orsequence elements which are conserved in a given set of strings, is an important
problem in the ﬁeld of Bioinformatics [1].
Earlier work on pattern discovery focus on discovering the most conserved
pattern in a given set of strings, generally preferring longer patterns which occur
in most of the sequences in the set. Another situation is when we are given two
sets of strings, where one set (positive set) consists of sequences known to possess
some biological characteristic, while the other (negative set) consists of sequences
known not to posses these characteristics. The problem is to ﬁnd a discriminating
pattern, that is, a pattern which occurs in most strings of the positive set, butdoes not occur in most of the strings of the negative set [2–6].
Recently, there have been several works which incorporate numeric attributes
which are obtained from other sources, e.g. gene expression data obtained from

microarray experiments, in order to ﬁnd meaningful patterns more eﬀectively [7–11]. The basic idea of these methods is to ﬁnd sequences elements whose occurrences in the sequences are correlated with the numeric attributes. For example,
gene expression is regulated by molecules called transcription factors, which bindto speciﬁc sequences usually in the upstream of the coding region of a gene. The
binding sites for a given transcription factor are fairly conserved across genes
which are regulated by the same transcription factor. Therefore, if we can ﬁnd
sequence elements which occur in upstream regions of genes which are relatively
highly expressed, while not occurring in upstream regions of genes whose expression is relatively low (or vice versa), such patterns are likely to be binding sitesof speciﬁc transcription factors.
In [7], substring patterns of up to length 7 are scored according to the linear
ﬁt between the number of occurrences in the upstream region and the expression
level of the gene. However, they do not consider any algorithm for solving the
problem eﬃciently. Also, the choice of the maximum pattern length is arbitrary
and it is not guaranteed that the optimal pattern will be found. Algorithmic
work for solving a similar problem has been considered in [8, 9]. In this problem
setting, the number of occurrences of a pattern is only considered as an indicatorvalue of 0 or 1, i.e. whether the pattern occurs in the string or not. Based on
the algorithm for solving the color set size problem [12], a very eﬃcient O(N )time algorithm for ﬁnding the optimal substring pattern in this problem setting
is given in [9]. A general branch and bound strategy that can be used for more
complex patterns where the problem can be NP-hard, is given in [8].
Although the algorithms above have been shown to discover similar motifs asin [7], it is generally believed that multiple occurrences of a binding site motif in
the upstream region of a gene will strengthen the function of the transcription
factor for that gene. In this paper, we give an eﬃcient algorithm to discover the
optimal pattern, taking into account the number of occurrences of the patternin each string, as in the problem setting in [7]. We ﬁrst present two simple
algorithms based on the suﬃx tree data structure that ﬁnds the optimal substring
pattern (without a restriction in the length of the pattern), respectively in O(nN )and O(N 2 ) time. We further develop and apply a branch and bound strategy
in order to speed up the algorithm, also allowing the problem to be solved formore complex and descriptive classes of patterns. The algorithms developed areapplied to real biological data to show the eﬃciency and eﬀectiveness of theapproach.

2Preliminaries

Let Σ 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 len(w). The empty string is denoted byε, that is, len(ε) = 0. For any set S, let |S| denote the cardinality of the set.
The empty set is denoted by ∅, that is, |∅| = 0. Let R represent the set of real
numbers.
2

Let Π be a set of patterns. We call a function deﬁned over a text string andpattern ψ : Σ ∗ × Π → R a matching function. Let ψp : Σ ∗ → R represent
the matching function for a ﬁxed p ∈ Π, that is, ψ(s, p) = ψp (s) for any text
string s ∈ Σ ∗ . For the matching function value, we shall consider the numberof occurrences of a given pattern in the text string. A substring pattern p is apattern p ∈ Π = Σ ∗ , where the matching function value ψp (s) is deﬁned as thenumber of substrings in s which is equal to p. A don’t care pattern p is a patternp ∈ Π = ({.}∪Σ)∗ , where “.” is a don’t care symbol, and the matching function
value ψp (s) is deﬁned as the number of substrings in s which can be obtained
from p by appropriate substitution of the don’t care symbols with characters of
the alphabet Σ. For the above two pattern classes, we shall refer to Σ or {.} ∪ Σas the pattern alphabet.
For the rest of the paper, we assume that we are given as input, a sequenceof ordered pairs consisting of a string and an associatednumeric attribute value:∑n{(s1 , y1 ), . . . , (sn , yn )} ⊂ Σ ∗ × R. Let N =len(si ) represent the total
i=1
T
length of the input strings. We denote by y = (y1 , . . . , yn ) ∈ Rn , the vector
consisting of the numeric attribute values. Further, for a given pattern p, weTdenote by ψp (s) = (ψp (s1 ), . . . , ψp (sn )) ∈ Rn , the vector consisting of thematching function values for the input text strings. We deﬁne for later use, apreorder over patterns as follows:
Deﬁnition 1. For any p0 , p ∈ Π, denote p0 ¹ p if for all {si | i = 1, . . . , n},ψ(si , p0 ) ≤ ψ(si , p).
We now consider how to score the goodness of a given pattern. For a givenT
x = (x1 , . . . , xn ) ∈ Rn , let X = (1, x), and deﬁne
RSS (y|x) = ||y − X β̂||2 =n∑(yi − (β̂0 + β̂1 xi ))2
i=1

T

where β̂ = (β̂0 , β̂1 ) are the least square estimates of β in the linear modely = Xβ + ².
We consider the problem of ﬁnding the pattern which can best ﬁt the numeric
attribute values yi with respect to ψp (si ).
Deﬁnition 2 (Pattern Based Linear Regression). We deﬁne the patternbased linear regression problem as follows. Given {(s1 , y1 ), . . . , (sn , yn )} ⊂ Σ ∗ ×
R, and a matching function ψ, ﬁnd the pattern p ∈ Π that minimizes
RSS (y|ψp (s)) = ||y − (1, ψp (s))β̂||2 =n∑(yi − (β̂0 + β̂1 ψp (si )))2
i=1

T

where β̂ = (β̂0 , β̂1 ) are the least square estimates of β in the linear modely = (1, ψp (s))β + ².
3

We note that the consistency problem [13,14] is a special case of our problem.
Since the consistency problem is shown to be NP-complete for several patternclasses (e.g. subsequence patterns), the above problem is NP-hard for such cases.
An exception is the case for the substring pattern class, for which we shall present
eﬃcient solutions in Section 3.

3

Methods

3.1

Finding the Optimal Substring Pattern

When considering substring patterns, it can be shown that the number of possible patterns which give distinct RSS scores is linear in the total length ofstrings, i.e. O(N ). We make use of a very convenient and well studied data
structure called suﬃx trees.
Generalized Suﬃx Trees A suﬃx tree [15] for a given string s is a rooted tree
whose edges are labeled with substrings of s, satisfying the following characteristics. For any node v in the 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
node v, l(v) is a distinct suﬃx of s, and for each suﬃx of s, there exists such aleaf v. Furthermore, each node has at least two children, and the ﬁrst characterof the labels on the edges to its children are distinct. A generalized suﬃx tree
(GST) for a set of n strings S = {s1 , . . . , sn } is basically a suﬃx tree for thestring s1 $1 · · · sn $n , where each $i (1 ≤ i ≤ n) 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 is labeled with idi . An example of a
GST is shown in Fig. 1. It is well known that suﬃx trees (and generalized suﬃx
trees) can be represented in linear space and constructed in linear time [15] withrespect to the length of the string (total length of the strings for GST).
Notice that candidate substring patterns may be restricted to those represented by nodes of the generalized suﬃx tree. This is because, for any substring
pattern that does not correspond to a path in the suﬃx tree, the pattern does
not occur in any of the strings in the set. Also, note that for a given pattern that
does correspond to a path in the suﬃx tree, all occurrences of the pattern in thestrings are represented by the leaves of the suﬃx tree in the subtree below thispath. This means that for any substring pattern that corresponds to a path that
ends in the middle of an edge of the suﬃx tree, its occurrences in the strings are
identical to the occurrences of the substring pattern corresponding to the pathextended to the next node.
As stated in the Introduction, the pattern based linear regression problem
has been shown to be solvable in O(N ) time if the matching function is considered to be an indicator function returning the value 0 or 1 [9]. However, it
is assumed that the score of a pattern is a function of the sum of the matching
function values and the sum of the numeric attribute values of the strings thatthe pattern occurs in. The algorithm cannot be applied to our case since we
4

a

b
b

id

$1

b a
b
$1 b
$2
$1

a
$2

b

a
b
$1 b
$2

$2

$2
$1

6

1

7

3

2

8

4

9

5

10

2

1

2

1

1

2

1

2

1

2

Fig. 1. Example of a generalized suﬃx tree for the string set { abab, aabb }.

require the matching function values of each string in order to compute the RSS
score. Below, we give two algorithms which calculate the score for each candidate
pattern p corresponding to a node in the generalized suﬃx tree.
An O(N 2 ) Algorithm Calculating the score for each of the O(N ) candidate
patterns requires the calculation of ψp (s) for each p, as well as the calculationof the least square estimates. The former can be calculated in O(N ) time using
well known linear time string pattern matching algorithms such as the KnuthMorris-Pratt algorithm [16]. The latter can be calculated by ﬁrst obtaining theleast square estimate of the parameters:β̂ = (X T X)−1 X T y,
where X = (1, ψp (s)). It is not diﬃcult to see that this can be calculated in
O(n) time. Further, RSS (y|ψp (s)) can be calculated in O(n) time, and therefore
the resulting time complexity is O(N ) · O(N + n) = O(N 2 ).
An O(N n) Algorithm Consider assigning a vector of length n at each nodeand leaf of the suﬃx tree. The vectors are initialized as follows: for an internalnode, all values are set to 0. For a leaf labeled with idi , the value at the ith
position is set to 1, and the rest is set to 0. Then, with a bottom-up (postorder)
traversal on the suﬃx tree, we add up the values in the vector at each node,5element-wise, into the vector of its parent node. The result is that we obtain, at
each node, a vector of length n where each position i of the vector represents thenumber of leaves in the subtree of the suﬃx tree, with idi . This corresponds to thenumber of times the substring pattern occurs in string si , and consequently, thevector represents ψp (s). This means that ψp (s) can be calculated in O(n) timefor each pattern, resulting in a total of O(N n) time for the score calculations.
3.2

Branch and Bound Strategy

Since the pattern based linear regression problem can be NP-hard when considering more complex pattern classes, we propose an enumerative branch andbound framework for ﬁnding the optimal pattern. The basic idea of the enumeration is similar to previous works [2–6, 8]. The main contrivance of this paper
is in the method for calculating the lower bound of the RSS score for speciﬁc
subspaces of the pattern space.
The algorithm proceeds by traversing and enumerating nodes on a search tree,
where each node in the tree represents some pattern in Π. For any pattern p in
the search tree, let p0 be a pattern represented by a node in the subtree rooted
at the node for p. While traversing the search tree at the node corresponding
to p, suppose that we are able to calculate a lower bound for the RSS for any
pattern p0 . If this lower bound is greater than the current best RSS found in
the traversal, we know that the score for p0 cannot be below the current best
RSS . This allows us to prune the search space by disregarding the subtree ofthe search tree rooted at the node corresponding to p.
Below, we show how such a lower bound can be calculated. The assumptions
for our calculations below is that p0 ¹ p. For the case of string patterns anddon’t care patterns, this assumption can be fulﬁlled by considering the searchtree described as follows. The root corresponds to the empty string ε, and eachnode v will have child nodes whose pattern corresponds to the pattern obtainedby extending the pattern at node v by one character of the pattern alphabet. Although we do not elaborate in this paper, the same branch and bound approach
can be applied to a variety of other patterns such as approximate patterns anddegenerate patterns.
Problem 1 (lower bound of the residual sum of squares for p0 ¹ p). Given somepattern p ∈ Π, ﬁnd a lower bound of the score function RSS (y|ψp0 (s)) for any
pattern p0 ¹ p.
T

T

opt
Below, let x = (x1 , . . . , xn ) = ψp (s). Also, let xopt = (xopt
1 , . . . , xn ) =Targ minz∈Dx RSS (y|z) where Dx = {(z1 , . . . zn ) | 0 ≤ zi ≤ xi , i = 1, . . . , n},and let β̂ opt = (β̂0opt , β̂1opt ) be the least square estimates for the linear modely = (1, xopt )β + ². Our objective now is to ﬁnd RSS (y|xopt ). This can beconsidered as a relaxed version of the problem stated above, since we do not
require that there exists a pattern p0 ¹ p such that ψp0 (s) = xopt . The following
theorem gives a simple lower bound on the RSS score.

6

Algorithm 1: Simple algorithm for calculating a lower bound of RSS .
Input: x = (x1 , . . . , xn )T ,y = (y1 , . . . , yn )TOutput: lb ≤ min{RSS (y|z) | z = (z1 , . . . , zn )T , 0 ≤ x0i ≤ xi }1 X ← {i | xi = 0};
2 if X = ∅ then return 0;/* Can move all xi onto a single line. */
∑
y
/|X|;
/* mean of values {yi | xi = 0} */
3 m ←
i
∑ i∈X
2
(m
−
x
)
;
/*
residualsum
of squares for points {xi = 0} */4 lb =ii∈X
5 return lb

Theorem 1 (simple lower bound). Let X = {i | xi = 0}. If X 6= ∅, then∑(m − yi )2 ≤ RSS (y|xopt )i∈Xwhere m =∑i∈X yi /|X|.

Proof. For any i, if xi = 0, then since 0 ≤ xopt
≤ xi we have that xopt
= 0.
ii
Ignoring the residuals for all data points x 6= 0, the minimum possible residualsum of squares for data where xi = xopt
= 0 must be smaller than the residuali
sum of squares for the entire data set.
ut
The simple lower bound can be calculated with the pseudo-code shown in Algorithm 1.
Next, we try to improve on this lower bound. The following lemma gives theconditions between xopt and the regression line.
Lemma 1. For all i = 1, . . . , n, if β̂1opt > 0 then{0
if xi = 0 or yi ≤ β̂0opt
xopt
=
opt
opt
i
(yi − β̂0 )/β̂1
otherwise (xi > 0 and yi > β̂0opt )and if β̂1opt < 0,
{
0
xopt
=
i
(yi − β̂0opt )/β̂1opt

if xi = 0 or yi ≥ β̂0opt
otherwise (xi > 0 and yi < β̂0opt ).

Proof. We will prove the case for β̂1opt > 0. The case for β̂1opt < 0 can be done
opt
similarly. The lemma states that the points (xopt
1 , y1 ), . . . , (xn , yn ) lie either
opt
opt
on the regression line y = β̂0 + β̂1 x, or on the y-axis. Suppose there exist
opt
points (xopt> 0 and the point is not on thei , yi ) to the contrary, that is, xi
opt
opt
regression line. If yi < β̂0 , then since β̂1 > 0, considering point (0, yi ) instead
optof (xopt.
i , yi ) would give a smaller residual, contradicting the deﬁnition of x
opt
opt
opt
If yi > β̂0 , then again since β̂1 > 0, (xi , yi ) cannot lie to the right of theopt
opt
regression line, or we can replace (xopti , yi ) with a point ((yi − β̂0 )/β̂1 , yi )on the regression line with a smaller residual. We can also say that the points
7

cannot lie left of the regression line, since we can construct a new regression
line passing the point (0, β0 ) that lies left of the points, and replace all points(xopt
i , yi ) with points on the new regression line. The residual of the new points
are clearly smaller, and again contradicts the deﬁnition of xopt .
utCorollary 1. Let X opt = {i | xopt
= 0}. If X opt = ∅, RSS(y|xopt ) = 0. Ifiopt
X
6= ∅, then∑β̂0opt =yi /|X opt |i∈X

andRSS(y|xopt ) =∑2(β̂0opt − xopti ) .

i∈X opt

Proof. As a consequence of Lemma 1. If X = ∅, then all points (xopti , yi ) lie on
the regression line. Otherwise, since the residual for all points xopt
> 0 is 0, β0iis chosen to minimize the residual sum of squares of points where xopt
= 0. uti
In order to calculate the lower bound, the problem now is how to obtain thevalue β̂0opt . Although the exact value of β̂0opt is not known beforehand, Algorithm 2 shows how to calculate the minimum residual sum of squares for anyTz = (z1 , . . . , zn ) satisfying the constraint: 0 ≤ zi ≤ xi for all i = 1, . . . , n.
Corollary 2. Let
opt= {i ∈ X opt | yi > β̂0opt }X+
opt= {i ∈ X opt | yi ≤ β̂0opt }.
X−

If β̂1opt > 0, thenopt
= {i | xi = 0, yi > β̂0opt }X+
opt= {i | yi ≤ β̂0opt }.
X−

If β̂1opt < 0, thenopt
X+
= {i | xi = 0, yi < β̂0opt }optX−= {i | yi ≥ β̂0opt }.

Proof. As a consequence of Lemma 1. When β̂1opt > 0 and yi > β̂0opt , xopt
= 0 if
i
opt
opt
opt
and only if xi = 0. Similarly, when β̂1 < 0 and yi < β̂0 , xi = 0 if and only
if xi = 0.
utTheorem 2. Algorithm 2 correctly outputs RSS (y|xopt ).
8

Algorithm 2: Algorithm for calculating the lower bound of RSS .
Input: x = (x1 , . . . , xn )T , y = (y1 , . . . , yn )TOutput: min{RSS (y|z) | z = (z1 , . . . , zn )T , 0 ≤ x0i ≤ xi }1 X ← {i | xi = 0}; k ← 1;
2 if X = ∅ then return 0 ;/* Can move all xi onto a single line. */opt// assumingβ̂1 > 0 ////////////////////////////////////////∑
y /|X|;/* mean of values {yi | xi = 0} */
3 m ←
i∈X i
4 (i1 , . . . , in−|X| ) ← indices sorted in increasing order of {yi | i 6∈ X} ;5 while yik ≤ m doX←∑
X ∪ {ik };/* xik → 0 */6m←
y /|X|;/* update mean */
7
i∈X i
8
k ← k + 1;
9 endw∑(m − xi )2 ;/* m = β̂0opt and X = X opt if β̂1opt > 0 */10 lb =i∈X
opt// assuming β̂1 < 0 ////////////////////////////////////////
11 X ← {i | xi = 0}; k ← 1;∑
y /|X|;/* mean of values {yi | xi = 0} */12 m ←i∈X i
13 (i1 , . . . , in−|X| ) ← indices sorted in decreasing order of {yi | i 6∈ X} ;14 while yik ≥ m do
15X←∑
X ∪ {ik };/* xik → 0 */
m←
y
/|X|;
/*
updatemean */
16ii∈X
17k ← k + 1;
18 endw∑(m − xi )2 };/* m = β̂0opt and X = X opt if β̂1opt < 0 */19 return min{lb,i∈X

Proof. Consider the case where {i | xopt
= 0} =
6 ∅ and β̂1opt > 0. The claim
iopt
is that m = β̂0 , and X = X opt , after the while loop of lines 5–9. If thiscan be proved, the result follows from Corollary 1. Let us split the set X thus
calculated into two disjoint sets, X+ = {i ∈ X | yi > m} and X− = {i ∈X | yi ≤ m}. Since the algorithm does not add indices i where yi > m toX, we have that X+ = {i | xi = 0, yi > m} from the initial construction of
X at line 1. Also, since all indices i where yi ≤ m are added to X, we havethat X− = {i | yi ≤ m}. Suppose m < β̂0opt . From Corollary 2, we have thatopt
opt
X+ ⊇ X +, and X− ⊆ {i | yi ≤ β̂0opt } = X−
. However, this contradicts the∑∑opt
assumption, since i∈X yi /|X| = m ≥ β̂0 = i∈X opt yi /|X opt |. Next, suppose
opt
opt
and X− ⊇ X−.
m > β̂0opt . We have the opposite situation and X+ ⊆ X+
However, due to the way X is constructed, there must have been a point in thewhile loop (lines 5–9) where all indices i where yi ≤ β̂0opt are added to X, andno index i where yi > β̂0opt is added to X. From Corollary 2, at such point,
opt
opt, which implies that m = β̂0opt . Due to the condition
and X− = X−
X+ = X+of the while loop, no other index i where yi > β̂0opt could have been added to X
afterwards, contradicting the assumption m > β̂0opt . Therefore, m = β̂0opt , andconsequently X = X opt .
u
t
9

y

avg at 0yshift to 0xx

Fig. 2. Case of Algorithm 2 where X = ∅ (left) and X 6= ∅, β1 > 0 (right).

Fig. 2 shows the basic idea of Algorithm 2. The time complexities of Algorithm 1 and 2 are both O(n), since they just conduct a constant number of scanson the data set, provided that the data is initially sorted and ranked according
to yi , so that the sorting at lines 4 and 13 of Algorithm 2 can be computed in
O(n) time.
Combining Suﬃx Tree Traversal and Branch and Bound It is easy to
see that for substring patterns, the generalized suﬃx tree itself can be usedas the search tree, and we can combine the O(N 2 ) algorithm and the pruning
strategy described in this section. Combining the O(N n) algorithm and thepruning strategy is not readily realizable, since the direction of the traversal
over the search tree is in the reverse direction. We discuss this issue further in
Section 5.

4

Computational Experiments

We implement our algorithms using the C++ language, and measure the runningtimes of our algorithm using a Sun Fire 15K (UltraSPARC III Cu 1.2GHz x
96 CPUs) using a single processor for each run. We note that the suﬃx tree
algorithm is simulated by a suﬃx array structure, using the method presented
in [17, 18].
For the numeric attribute values, we use the S. cerevisiae gene expression
data obtained from microarray experiments given in [19]. The data consists ofnormalized log expression level ratios of genes at speciﬁc time points of the yeast
cell cycle. For the string data, we use the 600 nucleotides from the upstream ofthe start codon of each gene.
Table 1 shows the running times of the algorithms for ﬁnding the optimal
substring pattern applied to the expression data of the 14-minute time point in
10

Table 1. Comparison of running times of algorithm for ﬁnding optimal substring
patterns from the 14-minute time point in the α-synchronized cell-cycle microarrayexperiment [19].
nO(N n)

100
500
1000
1500
2000
2500
3000
3500
4000
4500
5000
5500
5907

05.96s00m37s02m03s04m09s
07m53s12m30s
18m20s
25m17s
33m08s
42m48s55m03s67m01s75m55sO(N 2 )+O(N 2 )+bb
simple bb
12.43s 05.90s
05.71s01m47s 00m30s
00m26s05m28s 01m02s
00m53s
11m16s 01m35s
01m19s
19m55s 02m00s
01m38s30m05s 02m35s
02m07s42m17s 03m29s
02m51s
56m34s 04m11s
03m23s
73m06s 04m52s
03m55s
92m16s 05m34s
04m28s113m20s 06m13s
04m59s136m25s 07m08s
05m41s157m27s 07m57s
06m21sO(N 2 )

the α-synchronized cell-cycle microarray experiment. The times are measured forvarious sizes of n by a random sampling from the entire set of 5907 genes whichwere available for this time point. Note that since all strings are of ﬁxed length,
N = 600n. From the table, we see that the O(N n) and O(N 2 ) time algorithms
are able to ﬁnd the optimal pattern in a reasonable amount of time. However, we
can see that O(N 2 )+bb (the O(N 2 ) algorithm with the branch and bound strategy Algorithm 2) is much faster for all input sizes. Although O(N 2 )+simple bb(the O(N 2 ) algorithm with the simple branch and bound strategy Algorithm 1)
is fairly eﬃcient as well, the extra work invested in the O(N 2 ) + bb algorithm ispayed oﬀ by a ∼ 20% reduction in the overall computation time.
To show that the algorithm also allows for the discovery of more complex
patterns in a practical amount of time, we searched for the optimal don’t carepattern on the same data set. The search took 765 minutes and 529 minutes,
respectively, using simple bb and bb with a simple enumeration of don’t carepatterns, limiting the maximum length of the pattern to 15 and the number of
don’t cares characters in the pattern to 20% of its length.

5

Discussion

We presented eﬃcient algorithmic solutions to the problem of discovering the
optimal pattern in terms of a linear least squares ﬁtting of the numeric attributevalues associated with strings, and the matching function values. The eﬃciencyof the algorithms are conﬁrmed through computational experiments conductedon actual biological data.
11

In [20], the branch-and-bound enumerative search was applied in the reverse
direction for ﬁnding the optimal degenerate pattern that discriminates between
a positive string set and negative string set, where the matching function isan indicator function. In their search, the search tree is essentially traversed
bottom-up. A bound on the score is computed in a similar way as for the original
direction, and the traversal on the search can be pruned. This reverse directionis, however, diﬃcult to achieve for the problem considered in this paper. Thisis because in the original direction, the bound is calculated from the numeric
attribute values of strings which do not match the pattern (ψp (si ) = 0), whichis possible due to the condition ψp (si ) ≥ 0. In order to calculate a bound for thereverse direction, we would need to assume a maximum value c where ψp (si ) ≤ c,and we would calculate a bound for the residual sum of squares from the numeric
attribute values of strings which give ψp (si ) = c. However, the matching functionused in this paper is not suitable, since the maximum value would change foreach string.
The lower bound calculated in this paper underestimates the actual minimum
value that can be achieved with the matching function and numeric attributevalues. This is because we did not require that there exist a pattern whosematching function values would be equal to xopt . Notice that when calculating
the lower bound, xopt is obtained by considering points on the regression line.
However, we know that the matching function will only take discrete integer
values, and it may not be possible for some xopt
to lie on the regression line.
i
Residuals for such points will not be zero, and would therefore increase thelower bound. Finding an eﬃcient way to calculate a better lower bound for thediscretized version of the problem is an interesting open problem.

Acknowledgements
Computation time was provided in part by the Super Computer System, Human
Genome Center, Institute of Medical Science, University of Tokyo.

References
1. Brazma, A., Jonassen, I., Eidhammer, I., Gilbert, D.: Approaches to the automatic
discovery of patterns in biosequences. J. Comput. Biol. 5 (1998) 279–305
2. Hirao, M., Hoshino, H., Shinohara, A., Takeda, M., Arikawa, S.: A practical algorithm to ﬁnd the best subsequence patterns. Theoretical Computer Science 292(2002) 465–479
3. 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
4. Takeda, M., Inenaga, S., Bannai, H., Shinohara, A., Arikawa, S.: Discovering mostclassiﬁcatory patterns for very expressive pattern classes. In: 6th International
Conference on Discovery Science (DS 2003). Volume 2843 of LNCS., SpringerVerlag (2003) 486–493

125. 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 (DS2001). Volume 2226 of LNAI., Springer-Verlag (2001)
435–440
6. Inenaga, S., Bannai, H., Shinohara, A., Takeda, M., Arikawa, S.: Discovering
best variable-length-don’t-care patterns. In: Proceedings of the 5th International
Conference on Discovery Science. Volume 2534 of LNAI., Springer-Verlag (2002)
86–97
7. Bussemaker, H.J., Li, H., Siggia, E.D.: Regulatory element detection using correlation with expression. Nature Genetics 27 (2001) 167–171
8. 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
9. 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 2 (2004) 273–288
10. 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)
11. Bannai, H., Hyyrö, H., Shinohara, A., Takeda, M., Nakai, K., Miyano, S.: AnO(N 2 ) algorithm for discovering optimal Boolean pattern pairs. IEEE/ACM
Transactions on Computational Biology and Bioinformatics 1 (2004) 159–170 (special issue for selected papers of WABI 2004).
12. Hui, L.: Color set size problem with applications to string matching. In: Proceedings of the Third Annual Symposium on Combinatorial Pattern Matching (CPM92). Volume 644 of LNCS., Springer-Verlag (1992) 230–243
13. Miyano, S., Shinohara, A., Shinohara, T.: Which classes of elementary formal
systems are polynomial-time learnable? In: Proceedings of the 2nd Workshop on
Algorithmic Learning Theory. (1991) 139–150
14. Miyano, S., Shinohara, A., Shinohara, T.: Polynomial-time learning of elementary
formal systems. New Generation Computing 18 (2000) 217–242
15. Gusﬁeld, D.: Algorithms on Strings, Trees, and Sequences. Cambridge University
Press (1997)
16. Knuth, D.E., Morris, J.H., Pratt, V.R.: Fast pattern matching in strings. SIAM
Journal on Computing 6 (1977) 323–350
17. 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 2001). Volume 2089
of LNCS., Springer-Verlag (2001) 181–192
18. Kasai, T., Arimura, H., Arikawa, S.: Eﬃcient substring traversal with suﬃx arrays.
Technical Report 185, Department of Informatics, Kyushu University (2001)
19. Spellman, P.T., Sherlock, G., Zhang, M.Q., Iyer, V.R., Anders, K., Eisen, M.B.,
Brown, P.O., Botstein, D., Futcher, B.: Comprehensive identiﬁcation of cell cycleregulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization.
Mol. Biol. Cell 9 (1998) 3273–3297
20. Shinozaki, D., Akutsu, T., Maruyama, O.: Finding optimal degenerate patterns in
DNA sequences. Bioinformatics 19 (2003) ii206–ii214

13