March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

EFFICIENTLY FINDING REGULATORY ELEMENTS
USING CORRELATION WITH GENE EXPRESSIONHIDEO BANNAI1 , SHUNSUKE INENAGA2,3 , AYUMI SHINOHARA2,3 ,MASAYUKI TAKEDA2,3 , SATORU MIYANO1
1 Human Genome Center, Institute of Medical Science, University of Tokyo,

4-6-1 Shirokanedai, Minato-ku, Tokyo 108-8639, Japan{bannai,miyano}@ims.u-tokyo.ac.jp
2 Department of Informatics, Kyushu University 33, Fukuoka 812-8581, Japan
3 PRESTO, Japan Science and Technology Agency (JST){s-ine,ayumi,takeda}@i.kyushu-u.ac.jp

We present an efficient algorithm for detecting putative regulatory elements in the upstream DNA sequences of genes, using gene expression information obtained from microarray experiments. Based on a generalized suffix tree, our algorithm looks for motif
patterns whose appearance in the upstream region is most correlated with the expression
levels of the genes. We are able to find the optimal pattern, in time linear in the total
length of the upstream sequences. We implement and apply our algorithm to publicly
available microarray gene expression data, and show that our method is able to discoverbiologically significant motifs, including various motifs which have been reported previously using the same dataset. We further discuss applications for which the efficiency ofthe method is essential, as well as possible extensions to our algorithm.
Keywords: suffix tree; pattern discovery; gene expression and regulatory elements.

1. Introduction
Gene expression is regulated by transcription factors, which bind to specific sequences usually in the upstream of the coding region of a gene. The binding sites
for a given transcription factor are usually short regions of up to 15 base pairs,
and is fairly conserved across genes which are regulated by the same transcription
factor. In order to fully understand the mechanism underlying gene regulation, it iscritical to discover the motifs that correspond to transcription factor binding sites.
Microarray data provides measurements of gene expression under varying experimental conditions. One popular use of this data is to first cluster the genes
according to their gene expression measurements in several experiments – for example, within a time series – and to subsequently find motifs which are common
to the upstream region of genes within each cluster1,2,3 . However, this clustering
approach has inherent limitations: not all genes in the cluster possess a common
motif; more importantly, any one gene may possess multiple motifs that correspondto different transcription factors, so clustering the genes into disjoint sets may notallow these motifs to be found.
To overcome such problems, a method which does not require a pre-clustering
1March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

2

step has been presented by Bussemaker et al.4 Concerning a single microarray experiment, they assume a linear relation between the log expression ratio of a gene,and the number of times a motif appears in the upstream region of the gene. Theirmethod exhaustively looks for motifs of a certain length which best fits the expression data, using a simple motif model which does not consider any mismatches orambiguity (which we shall call the substring pattern class). In order to enumerate all
possible patterns, a naı̈ve solution for this problem would result in exponential timein the length of the pattern to look for. However, that study presented no rigorous algorithm for finding the optimal pattern. Also, since the length of the pattern
to look for is predefined, optimality of the patterns is not assured. Furthermore,
finding the optimal pattern with such a naı̈ve enumeration is virtually impossible
because it is limited only by the longest sequence in the set of sequences.
The objective of this paper is to present a new algorithm that overcomes these
difficulties: an algorithm that runs in linear time and gives optimal results. The
problem setting described above can be formulated as a case of the string patternregression problem which we have defined in a previous paper5 , provided that we
consider the occurrence of the motif as a binary value - the motif either occurs
or it does not. The rationality behind this simplification is based on the following
observations: 1) we do not know a priori, if the log expression ratio of the gene andthe number of times that a motif appears in the upstream region are linearly related,and 2) a motif of biologically plausible length will only appear a few times, if at all,in each upstream sequence. Thus, a binary indicator is likely to be a good estimateof occurrence counts anyway. We show that our method gives results similar to
those obtained previously4 , but is remarkably fast.
We previously gave an enumerative branch-and-bound algorithm5 which solves
the string pattern regression problem optimally for a general motif model including patterns with don’t care characters, wild-cards, mismatches, etc. However, theproblem is still difficult to solve in general, and the computational cost is still quite
large when considering gene regulatory regions which can be 600 ∼ 800bp long. Ourlinear time algorithm takes a completely different approach, based on a linear timealgorithm to solve the color set size problem, which finds for all possible substring
patterns appearing in a set of strings, the number of strings in the set which possess
the pattern6,7 . We extend this algorithm so that the expression level ratio of thegenes can be taken into account.
In Section 2, we give the problem definition and a brief introduction to background concepts. Section 3 describes the linear time algorithm for optimally solving
the string pattern regression problem for the substring class. In Section 4, we show
results of applying our algorithm to publicly available microarray data8 , and show
its effectiveness, as well as its efficiency. Finally, we discuss other applications andextensions of our algorithm in Section 5.

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

3

2. Preliminaries
Notation
Let R be the set of real numbers. Let Σ be a finite set of characters called thealphabet. A concatenation of the characters of the alphabet is called a string, andlet Σ∗ denote the set of strings of finite length. The length of a string w is denotedby l(w). The empty string is denoted by ε, that is, l(ε) = 0. Strings x, y, and z are
said to be a prefix, substring, and suffix of string w = xyz, respectively. For any set
T , let |T | denote its cardinality.
A pattern class is a pair C = (Π, m), where Π is a set called the pattern set andm : Σ∗ × Π → {−1, 1} is the pattern matching function. An element p ∈ Π is called
a pattern. For a pattern p ∈ Π and string t ∈ Σ∗ , we say p of class C matches t ifm(t, p) = 1, and p of class C does not match t, otherwise.
For example, the substring pattern class that we will use in this paper, is definedwith the pattern set Σ∗ and the pattern matching function strstr given by:½
1 if p is a substring of t
strstr(t, p) =.
−1 otherwise
For a pattern p, we denote by LC (p) = {t ∈ Σ∗ : m(t, p) = 1} the set of strings inΣ∗ which p of class C matches.
Problem Definition
We assume that we are given a set of data having two attributes: a string attributeand an objective numerical attribute. Let D ⊂ Σ∗ × R denote the data set.
For a pattern p ∈ Π, the data set D can be split into two sets: Dp = {(s, r) ∈D | s ∈ LC (p)}, which represents the subset of D for which the pattern p matches
the string attribute, and Dp = D −Dp , where p does not match the string attribute.
The aim of our approach is to find the pattern for which the values of the objective
numeric attribute of Dp is “best distinguished” from that of Dp . There can be
many measures for the goodness of such splits. One possibility is the measure usedin this paper, which is to minimize the mean squared error. We define the problem
as follows:
Definition 1 (string pattern regression5 ). Given a data set D ⊂ Σ∗ × R anda pattern class C = (Π, m), the string pattern regression problem is to find p ∈ Πsuch that the mean squared error defined below is minimized:P
M SE(D, p) =P(s,r)∈Dp (r − µ(Dp ))

2

+

P

(s,r)∈Dp (r − µ(Dp ))2|D|

where µ(D0 ) = (s,r)∈D0 r/|D0 | represents the average of the objective numericalattribute values of the data for any D0 ⊆ D.

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

4

Minimizing the mean squared error can be achieved by maximizing the interclassvariance (ICV )9 , defined by:
ICV (D, p) = |Dp |(µ(D) − µ(Dp ))2 + |Dp |(µ(D) − µ(Dp ))2 .
Without loss of generality, we can subtract the average of all the numeric attributes
from each numeric attribute so that their average will be 0, that is, µ(D) = 0.
P
Under this conversion, let xp = |Dp | and yp = (s,r)∈Dp r for a pattern p. Sincey−ypµ(Dp ) = xpp and µ(Dp ) = |D|−x, the interclass variance can be written as:p(0³´ if xp = 0 or xp = |D|ICV (D, p) = f|D| (xp , yp ) =1
1
2
yp xp + |D|−xp otherwise(1)

This measure is convenient for our calculations, and will be used for the remainingof the paper.
The string pattern regression problem is NP-hard in general, since it can containthe consistency problem10,11 as a special case. A branch-and-bound algorithm whichsolves the string pattern regression problem exactly has been given in a previous
work5 . It is applicable to various classes of patterns such as subsequence patterns,
approximate patterns, patterns with wild-cards (VLDC patterns12,13 ), etc. However,
we will show in Section 3 that for the substring pattern class, the problem can besolved in linear time in the total length of the string attributes, by a clever use ofgeneralized suffix trees7 .
Suffix Tree
A suffix tree7 of string s is a rooted tree consisting of nodes (internal node or leaf)and edges, where each edge is labeled by a non-empty substring of s. Each internal
node has at least 2 children, and the first character of the label on the edges to
its children must be distinct. For a node v, the string obtained by concatenating
the labels on the path from the root to v is denoted by p(v). For any leaf v, p(v)must represent a distinct suffix of s. Note that for any substring of s, the suffix tree
will contain such a path from the root, which may end at a node, or perhaps onan edge. If v is an internal node, p(v) represents a common prefix of suffixes of s(corresponding to leaves in the subtree below the node).
Figure 1 (A) shows an example of a suffix tree for the string “cocoa”. The
following algorithmic result is well known, assuming an alphabet of fixed size:
Lemma 1 (suffix tree construction14,15 ). A suffix tree for a string s can beconstructed in time linear in the length of the string.
It is also possible to construct a generalized suffix tree (GST), that represents
all suffixes of a set of strings, in time linear in the total length of the strings.
A conceptually simple way to do this would be to create a suffix tree for a stringobtained by concatenating each of the strings in the set with delimiters (e.g. for S ={s1 , . . . , sk }, create the suffix tree for string s1 $1 s2 $2 · · · sk $k , where $i (i = 1, . . . , k)

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

5

c
o

a
c
o
a

a
o

a

$
b

1,2

a
$

b
b $

c
o
a

a b

2

1

a
$
1

$
2

b
b
$

$

$
a
$

2

1
1

2

(A)

(B)

Fig. 1. (A) A suffix tree for the string “cocoa”. (B) A generalized suffix tree (GST) for the strings“abba$”(1) and “babb$”(2). The character ‘$’ represents a unique delimiter character that does
not appear in any other part of the strings. Each leaf on the GST contains a list of numbers
corresponding to the strings which contain that suffix. Note that the actual representation of thelabels on an edge of a suffix tree is not as depicted, but is a pair of integers that identifies asubstring of the given string.

are delimiters which do not appear in any of the strings). Since all strings end witha distinct character, and we do not need to consider substrings spanning across thedelimiters, the tree is not constructed beyond the delimiters (i.e. paths end once a$i appears). Therefore, each leaf in the GST corresponds to a $i ending the path(which in turn corresponds to a string in the set) and is thus labeled. Although weomit details, it is possible to avoid the apparent increase in the alphabet size by
considering all $i ’s as the same character $ and marking i on the leaves. An example
of a GST is shown in Figure 1 (B). Note that the number of nodes and edges of
the generalized suffix tree is linear in the length of the concatenated string, andtherefore traversals on the tree can also be done in time linear in the total length.

Lowest Common Ancestor
The lowest common ancestor of any two nodes in a tree is the first node at which
the paths from the nodes (to the root) intersect. The following is also a well known
result concerning the computational complexity of lca queries:

Lemma 2 (efficient lca query16,17 ). A query for the lca of any two nodes ofa given tree can be answered in constant time, with linear (in the size of the tree)
time preprocessing.

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

6

3. Linear Time Algorithm for Solving the String Pattern
Regression Problem for the Substring Pattern ClassFor the given data set D ⊂ Σ∗ × R, let S = {s | (s, r) ∈ D} denote the k stringsP
in D (k = |D|), N =s∈S l(s) the total length of the strings in S, and T thegeneralized suffix tree of S. The string pattern regression problem is equivalent tofinding the node v̂ in T , which corresponds to the pattern p(v̂) giving the maximum
ICV , that isv̂ = arg max{ICV (D, p(v)) | v is a node of T }.
v

It is sufficient to consider only the O(N ) patterns represented by the internal
nodes of T , because:• patterns which do not have a corresponding path on the tree do not matchany of the strings in S, and therefore do not need to be considered.
• patterns corresponding to paths which end on an edge of the tree can beextended to the next node at the end of the edge, and will still match thesame set of strings, and therefore have the same ICV score.
From Equation (1), the ICV for a node v in T can be calculated if we know thenumber of data xp(v) , as well as the sum of the numeric attributes yp(v) of the data,
for which the pattern p(v) matches. It has been shown6 that xp(v) for all v canbe computed in linear time. This fact alone does not give us yp(v) in linear time,
but careful examination of the algorithm reveals that it can indeed be extended forincorporating such sums of real valued data associated with each string.
We first describe this algorithm by Hui6,7 that calculates xp(v) in O(N ) time.
xp(v) is the total number of distinctly labeled leaves in the subtree of v. To calculate
this value, one can imagine a single bottom-up traversal from all of the leaves, where
each leaf node will have a value of 1, and each internal node will have the summed
value of its child nodes. However, if we simply sum the values as we ascend, wewould be adding “too much”, that is, the value calculated for each node would not
be representing the number of distinctly labeled leaves in the subtree of the node.
This is because if there is more than one leaf with the same label in the subtree,
it would be accounted for more than once. A simple way to prevent this would
be to prepare a k bit vector at each node, where each bit corresponds to a leaf
label, representing whether or not the subtree contains a leaf with the label (i.e.
whether the pattern represented by the node occurs in the string corresponding to
the label). At each ascent, we would ‘or’ the bit vectors. However, this results in a
O(kN ) algorithm. We can do away with the k bit vector if we do the bottom-up
traversal k times, where, for the ith (1 ≤ i ≤ k) traversal, we would start from only
the leaves corresponding to si , and only add the values once to a given node in eachtraversal. However, this requires k traversals of O(N ) each, which is still O(kN )time. O(kN ) can be written as O(k 2 ¯l), where ¯l is the average length of the strings,
and can be computationally intensive when k becomes large.

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

7

The idea to avoid the extra factor of k is to calculate exactly how much was“too much”, and subtract it from the total sum including the redundancy. We shallcall this the correction factor7 . Let L be the list of leaves in the order that appears
in a depth-first traversal of the generalized suffix tree, and let Li (1 ≤ i ≤ k) be thesubsequence of leaves in L corresponding to the string si , in the same order. Notethat with the depth-first traversal, the leaves in the subtree of a node v appear as asingle contiguous sub-list in L. This means that the total number of times a subtreeof v contains the lca of a consecutive pair of leaves in Li , is equal to one less thenumber of leaves of Li which are in the subtree of v. Therefore, by calculating thelca for each consecutive pair of leaves in Li for each i (1 ≤ i ≤ k) and counting how
many times each node was selected as an lca, the correction factor at each node canbe computed with a single bottom-up traversal from the leaves. This is illustrated
in Figure 2.
lca(3,6).=?
.=?

(10,8)

lca(6,7)lca(5,8)(7,5)
lca(1,2)lca(7,9)

(3,2)

(4,2)
lca(4,5)lca(2,3)12

3

lca(8,10)

(3,1)

(2,1)

(2,0)
4

5

6

7

8

(3,1)

9

10

Fig. 2. Calculating the correction factors with lca queries on a GST. The GST represents two
strings, identified with the box and oval labels at the leaves, and the lca for consecutive nodesof each string is depicted. The pair of numbers (x1 , x2 ) at each node is calculated by a single
bottom-up traversal of the tree, where x1 represents the total number of leaves, and x2 represents
the total number of nodes that were counted as lca’s, in the subtree starting with the node. The
number of distinct strings that the pattern corresponding to the node matches, is obtained byx1 − x2 .

Now that we are ready, the algorithm can be summarized as follows: Assume
each node holds a pair of integers (x1 , x2 ), where x1 is the total number of leaves
in its subtree, and x2 is the correction factor for that node. A leaf should have
(x1 , x2 ) = (1, 0), while internal nodes are initialized to (x1 , x2 ) = (0, 0).
(1) Construct GST T for the set of strings S = {s1 , . . . , sk } (O(N )).
(2) Do a depth first traversal on T to make k lists of leaves (L1 , . . . , Lk ) , one list

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

8

for each string si (1 ≤ i ≤ k) (O(N )).
(3) Preprocess T for lca queries (O(N )).
(4) For each list Li (1 ≤ i ≤ k) obtained in step 2, calculate the lca’s for consecutive
nodes in the list. Each time a node is chosen as an lca, increment the correction
factor x2 of the node by 1 (O(N ) since the total length of the lists is O(N ),and each query is constant time).
(5) Do a bottom-up traversal from the leaves, adding x1 and x2 of the children to
those of its parents (O(N )).
Since all steps can be conducted in O(N ), the total complexity is also O(N ).
We next show that the above algorithm can be extended to incorporate numeric values associated with each string. Each node will hold another pair of values(y1 , y2 ), where y1 is the sum of all numeric values for all leaves in the subtree ofthe node, and y2 is the correction factor for the sum. Steps 4 and 5 are modified asfollows:(4’) Calculate the lca’s for consecutive nodes in the list obtained in step 2. Eachtime a node is chosen as an lca, 1) increment the correction factor x2 of thenode by 1, and 2) add the numeric attribute value of the corresponding string
to y2 of the node.
(5’) Do a bottom-up traversal from the leaves, adding x1 , x2 , y1 , and y2 of thechildren to those of its parents.
We can calculate the ICV for each node with fk ((x1 − x2 ), (y1 − y2 )), and we
can find the node which maximizes this value, by a final traversal of all nodes inthe tree which can be done again in O(N ). Altogether, we have solved the string
pattern regression problem for the substring pattern class in O(N ) time.
Theorem 1. The string pattern regression problem can be solved in O(N ) time,
for the substring pattern class, where N is the total length of the string attributesof the input.
Note that although we defined the score to optimize as ICV in Equation (1),any other scoring function of xp and yp can be used.
Finding Multiple Motifs
Assuming a model where each motif contributes additively to the numeric value, we
can perform a greedy iteration for finding multiple motifs4 , by first finding the best
pattern which contributes most to the numeric value, then subtracting the effectthat the pattern has on the numeric attributes. The procedure is then repeated withthe modified numeric attributes to find the next pattern with the most influence.
More precisely, for the best pattern p̂, we subtract µ(Dp̂ ) from the numeric
attributes of data in Dp̂ , and µ(Dp̂ ) from numeric attributes of data in Dp̂ , andfinally adjust the values so that the µ(D) = 0 again. Notice that the suffix tree

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

9

does not change between iterations, and can be reused for the next iteration. Also,
since the tree does not change, the preprocessing for the lca queries need only be
conducted once. Therefore, subsequent iterations can be computed much faster than
the first.
4. Application to Gene Expression Data
We show the effectiveness and efficiency of our algorithm by applying it to publicly
available gene expression data by Spellman et al.8 , which are microarray data obtained from time course experiments concerning the cell cycle of S. cerevisiae. We
implemented the algorithm in the Objective Caml language18 . We use the on-line
suffix tree construction algorithm of Ukkonen15 , and an efficient lca query method
by Bender and Farach-Colton17 , based on range minimum queries. Runtime performance was measured on a dual Xeon 2.2GHz PC with 2GB of memory.
Motifs Detected
The same dataset has been analyzed previously by several other methods, and weevaluate the motifs discovered by our algorithm by comparing it to the output of apreviously proposed, very similar method4 mentioned in Section 1. We give here theresults for the 14-minute time point in the α-synchronized cell-cycle experiment8 .
There were 6177 genes in the combined data, but we discard genes with no expression ratio measurements, as well as those whose sequence we could not retrieve
from SCPD19 . For the 5976 genes that could be used, we take the 600 (−600 ∼ −1)
basepair sequence upstream of the gene from SCPD, and pair them with their log
expression ratios. We run the algorithm for 20 iterations (i.e. output 20 motifs),
only considering motifs whose lengths are less than or equal to 7 (the same length
as used in the previous analysis4 ). Although we need not limit the motif length since
we are able to obtain the optimal pattern in O(N ) time, we limited our output sinceallowing longer motifs seemed to make the algorithm capture very long motifs common in very few genes (2 or 3) whose expression level ratios were similarly relatively
higher or lower than the rest.
Table 1 shows the first 20 motifs discovered through the iterative procedure.
Notice that a positive value for yp suggests that the motif enhances transcription,
whereas a negative value suggests that it inhibits transcription. Comparing withthe 11 motif result4 obtained from the same expression data, we mark the motifs
with edit distance less than or equal to 1, from the closest motif of the 11. We also
mark sup (resp. sub) if the string is a superstring (resp. substring) of a motif in the
11. We can see that 4 exactly equivalent motifs appear, as well as 6 others whichwere very close (although AAGGGG and CAGGGG correspond to the same motif
AGGGG), showing that the motifs obtained from our method are comparable tothose obtained from the previous method4 . Motifs corresponding to known elements,
MCB (ACGCGT), STRE (AGGGG and CCCCT), and SFF (GTMAACAA), arealso shown in the table.

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

10

Table 1. The first 20 motifs detected from the 14-minute time point in the α-synchronized cellcycle microarray experiment of Spellman et al. The superscript on the patterns represent the edit
distance from (sup (resp. sub) represents that it is a superstring (resp. substring) of) the closest
motif in the 11 motif result of Bussemaker et al., based on the same data.

ICV

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

AAATTTT1
ACGCGT0(MCB)
AAGGGGsup (STRE)CGATGAG0
CTCATCG0TGACGCGsup
TGAAAAAsup
ATAAGGA
CCCCsub(STRE)
CCTGGAA
TTCAAAA
AAAAGGTAAACAA0(SFF)
CAGGGGsup (STRE)
TTTTTCGCTGGGTCGAACCACTGGGCTGGCACACCATCTCA

18.12
12.54
9.553
8.27
5.54
4.28
3.45
3.30
3.00
2.62
2.41
2.40
2.29
2.28
2.20
2.10
1.96
1.98
1.90
1.87

total
occurrence

pattern p

iteration

1845
414
662
348
326
106
2018
396
4877
189
1134
2845
761
263
6582
95
121
94
98
363

xp

yp

1549
371
616
331
315
104
1678
387
2669
186
1020
2227
698
247
3751
94
121
94
97
321

−144.2+66.06+72.65−50.86−40.67+20.91−64.55+34.54+66.56+21.74−45.18+57.97−37.56+23.24−55.40−13.94+15.25−13.52+13.47−23.82

From gene expressions of other microarray experiments at different time points(details not shown), we were able to find other known motifs, such as the SCB
element (CGCGAAA). We also obtained ATGCGAA, a motif similar to the ‘histone’ motif (ATGCGAAR)8 in the 8th iteration of the 35-minute time point of
α-synchronized cell-cycle experiment, while it is reported that it was not detectedin the previous analysis4 .
We note that the sequence data we used seems to be slightly different fromthose used in the previous analysis4 , since the same patterns showed different total
occurrence counts.
Runtime Performance
For the experiments described in the previous subsection, the time for the first
iteration took around 102 seconds: the construction of the suffix tree around 65
seconds, and preprocessing for the lca queries around 15 seconds. The rest of theMarch 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

11

time was spent calculating the correction factors and finding the node giving the
optimal pattern. Each subsequent iteration takes around 19 seconds each, for a totalof 467 seconds for the 20 iterations.

Runtime of linear algorithm with different # of data
500
1-iteration
20-iterations
b&b: 1-iteration

450
400

runtime (seconds)

350
300
250
200
150
100
50
0
0

1000

2000

3000
# of data

4000

5000

6000

Fig. 3. The running times on a random subset of the total data, varying the sizes of the subset. The
times for one iteration (filled box), 20 iterations (cross) are shown for the linear time algorithm(finding the optimal pattern), and the time for one iteration (clear box) using the branch-andbound algorithm by Bannai et al. with the string pattern class (limiting the pattern length to7) is shown. Our linear time algorithm is much faster for one iteration. The difference grows for
subsequent iterations, since we can reuse the suffix tree structure between iterations.

To demonstrate that our algorithm actually runs in linear time, we also measured
the runtimes of our algorithm when applied to random subsets of the total set ofgenes. Figure 3 shows the runtime of our algorithm when we vary the sample size,
and we can see that the runtime grows linearly, as expected.
The runtimes using a naı̈ve string matching algorithm with the branch-andbound enumeration of patterns5 is also shown. The maximum length of the patternin the enumeration was limited to 7. This provides an optimistic lower bound forthe computation time of the previous method4 when using an enumerative strategy.
Note that this does not include the time for calculation that would be required forthe linear fitting between the occurrence counts and expression ratios. Also, it is
not apparent whether the branch-and-bound strategy can be applied to the originalMarch 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

12

scoring function.
The branch-and-bound algorithm is fairly efficient, but our linear time algorithm
is at least 4 times faster for one iteration. Furthermore, since we can reuse the suffix
tree between iterations, we can do 20 iterations with the linear time algorithm while
we only do one iteration with the branch-and-bound algorithm, meaning that our
new algorithm is almost 20 times faster.
5. Discussion
We present several possible extensions to our algorithm which can be made without
largely sacrificing its efficiency.
5.1. Scoring Function
In Section 4 we noted that when not limiting the pattern length, the algorithm
seemed to come up with very long patterns common to the upstream of 2 or 3 genes.
In fact, some genes seemed to have exactly the same upstream regions. Although
the biological meanings of such long common sequences remain to be investigated,
perhaps choosing a more appropriate scoring function could prevent such effects if
it is deemed undesirable. Note that this is an interesting observation which could
be made since our method is able to find the optimal pattern with respect to acertain scoring function, having no limit in the length of the pattern. As noted in
Section 3, although we defined the score to optimize as ICV in Equation (1), anyother function of xp and yp can be used. Also notice that for each node v of thegeneralized suffix tree, the length of the pattern lp = |p(v)| can be recorded on eachnode in linear time by a simple top down traversal of the tree. The scoring function
could also incorporate this information and take the form f (xp , yp , lp ). Furthermore,
if some assumptions on the monotonicity of f with respect to lp can be made (that is,
if for any l : l1 ≤ l ≤ l2 , we have max f (xp , yp , l) = max{f (xp , yp , l1 ), f (xp , yp , l2 )}),
it is not difficult to see that the string pattern regression problem for scoring function
f can also be solved in linear time, provided that f can be computed in constant
time.
Another possible extension could be in the number of microarray experiments.
We have only considered the case where we use one gene expression data per run,but by a suitable extension to the scoring function, our method can use multiple
microarrays obtained at different conditions. The data for each gene would consistof the upstream sequence together with a vector of gene expression ratios, and thealgorithm would find the pattern p with optimal score f (xp , yp ), where yp is thesum of the gene expression ratio vectors of genes which contain p. The runtime
would increase by a factor of the number of microarray data used.
5.2. Substring Patterns which Match in Reverse Complement
It is known that for several transcription factors, the binding site can occur oneither strand of the upstream DNA region. Therefore, it can be biologically mean-

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

13

ingful to consider pattern matching functions where a string matches if the string
itself is contained as a substring, or its reverse complement (i.e. exchange A ↔ T,
C ↔ G, and reverse the direction of the pattern) is contained as a substring of thesequence. This is partially indicated from the results of Table 1, where we can see,for example, that (2:“ACGCGT”) is self-complement, and also (4:“CGATGAG”)
and (5:“CTCATCG”), are complements of each other.
Such reverse complement patterns can be handled in our algorithm as follows:when given a set of strings S = {s1 , . . . , sk }, construct the generalized suffix tree
for the string s1 $1 s01 $1 s2 $2 s02 $2 · · · sk $k s0k $k , where s0i represents the reverse complement string of si , and $i (i = 1, . . . , k) are delimiters which do not appear inany of the strings. The tree is not built deeper than any $i . Since xp(v) and yp(v)of a node v depends on the number of distinct $i ’s which appear in the subtree ofv, running the algorithm on the new GST will produce the correct results as in theoriginal case. The time complexity will increase by a factor of 2, but is still linear
in the length of the total length of the input sequences.
Table 2 shows the first 20 motifs obtained with this extension, run on the samedata as in Table 1. Similar motifs are found again, but with some notable additions:the motif (13: “CGGGTAA/TTACCCG”) which is known as the REB1 binding site,and (16: “AACCCA/TGGGTT”) which is known as the RAP1 binding site.
5.3. Other Applications
The algorithm described in this paper was successfully employed20 as a subroutine
to extend a previous Bayesian network inference method21 that estimates a networkof genes from microarray data.
We roughly describe the method below: First, estimate an initial network using
the previous method21 . Consider a certain gene g as a possible transcription factor,
and use our algorithm to find a motif in the upstream regions that is most correlatedwith the likelihood (which is a numeric which depends on the structure of the current
network) that a given gene is a child of g. The rationale for this is that the initial
gene network estimation is fairly good, and genes which are truly the child of gshould have relatively higher likelihoods, while other genes should have relatively
lower likelihoods. Our motif detection algorithm finds a sequence common to genes
with higher likelihood of being regulated by g, and not appearing in genes withlower likelihood. After the motif is obtained, increase the prior probability thata gene which possesses this motif is actually a direct child of g. After updating
the prior probabilities, the network is re-estimated, and the process repeated untilconvergence.
The integration of network inference and motif discovery is shown to improve
the accuracy of the inferred network, compared to the inference method which does
not use the sequence information. The key to this integration is that our method
can consider the inaccuracy of the network structure with the likelihoods. Also,
since the inference of the network involves numerous iterations of motif discovery

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

14

Table 2. The first 20 motifs detected from the 14-minute time point in the α-synchronized cell-cycle
microarray experiment of Spellman et al, allowing reverse complement matches. The superscripton the patterns represent the edit distance from (sup represents that it is a superstring of) the
closest motif in the 11 motif result of Bussemaker et al., based on the same data.

iteration
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

pattern p

CGATGAG0 /CTCATCG0ACGCGT0 /ACGCGT0 (MCB)
AAAATTT0 /AAATTTT1
AGGGG0 /CCCCT0(STRE)CGCGTCA/TGACGCGsup
ATTTTTCsup /GAAAAATCATCTCA/TGAGATG
ATGCAGG/CCTGCATCATCGCA/TGCGATGGGCACAC/GTGTGCC
AAAACAA/TTGTTTTAAAAAAT/ATTTTTT1
CGGGTAA/TTACCCG (REB1)CAGGGAG/CTCCCTG
TAAACAA0 /TTGTTTA (SFF)
AACCCA/TGGGTT(RAP1)AAGCTCG/CGAGCTTGAGGAGA/TCTCCTCCAGTGAG/CTCACTG
CCTGGAA/TTCCAGG

ICV

total
occurrence

xp

yp

20.14
12.42
10.82
8.98
5.75
5.56
3.52
3.24
3.11
2.80
2.77
2.59
2.53
2.50
2.30
2.31
2.39
2.20
2.00
1.92

674
828
3755
3155
185
2990
645
286
483
187
2820
5560
676
157
1449
1371
260
509
238
370

630
371
2252
2346
170
2359
564
278
450
183
2140
3497
639
156
1259
1211
254
486
230
351

−106.5+65.73−123.3+113.2+30.81−89.12−42.37+29.31−35.97+22.27−61.73−61.36−38.00+19.47−47.80−47.26+24.11+31.32−21.03+25.19and network structure inference, the efficiency of the motif discovery algorithm isessential.
6. Conclusion
We presented an efficient, linear time algorithm which can be used for finding putative regulatory elements in the upstream regions of genes, by combining sequence
information together with the information of gene expression obtained from microarray experiments. This simple method seems to balance algorithmic efficiency
with biological relevance: we were able to rediscovery several biologically meaningful
motifs very efficiently. We have shown that the same algorithm can be extended inseveral ways without sacrificing its efficiency, and have successfully demonstrated avariation of the algorithm considering reverse complement matches. We also mention that the algorithm has also been shown to be effective in other situations where
it is required to find, efficiently, patterns whose appearance is most correlated with

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

15

a numeric value.
Acknowledgments
The authors would like to acknowledge the anonymous referees for their helpfulcomments. This research was supported in part by Grant-in-Aid for Encouragementof Young Scientists (B), and Grant-in-Aid for Scientific Research on Priority Areas(C) “Genome Biology” from the Ministry of Education, Culture, Sports, Scienceand Technology of Japan.
References
1. S. Tavazoie, J. D. Hughes, M. J. Campbell, R. J. Cho and G. M. Church, “Systematic
determination of genetic network architecture,” Nat. Genet. 22, 281–285 (1999).
2. J. Vilo, A. Brazma, I. Jonassen, A. Robinson and E. Ukkonen, “Mining for putative
regulatory elements in the yeast genome using gene expression data,” In Proc. of the
8th International Conference on Intelligent Systems fo Molecular Biology (ISMB),
384–394 (2000).
3. U. Ohler and H. Niemann, “Identification and analysis of eukaryotic promoters: recent
computational approaches,” Trends Genet. 17, 50–60 (2001).
4. H. J. Bussemaker, H. Li and E. D. Siggia, “Regulatory element detection using correlation with expression,” Nat. Genet. 27, 167–171 (2001).
5. H. Bannai, S. Inenaga, A. Shinohara, M. Takeda and S. Miyano, “A string patternregression algorithm and its application to pattern discovery in long introns,” Genome
Informatics 13, 3–11 (2002).
6. L. Hui, “Color set size problem with applications to string matching,” In Proc. ofthe 3rd Annual Symposium on Combinatorial Pattern Matching (CPM) LNCS 644,
230–243 (1992).
7. D. Gusfield, Algorithms on Strings, Trees, and Sequences. Cambridge University Press,
1997.
8. P. T. Spellman, G. Sherlock, M. Q. Zhang, V. R. Iyer, K. Anders, M. B. Eisen,
P.O. Brown, D. Botstein and B. Futcher, “Comprehensive identification of cell cycleregulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization,”
Mol. Biol. Cell 9, 3273–3297 (1998).
9. Y. Morimoto, H. Ishii and S. Morishita, “Efficient construction of regression trees withrange and region splitting,” Machine Learning 45, 235–259 (2001).
10. T. Jiang and M. Li, On the complexity of learning strings and sequences, In Proc. 4th
ACM Conf. on Computational Learning Theory, 367–371, (1991).
11. S. Miyano, A. Shinohara, T. Shinohara, Which classes of elementary formal systems
are polynomial-time learnable?, In Proc. 2nd Workshop on Algorithmic Learning Theory, 139–150 (1991).
12. S. Inenaga, H. Bannai, A. Shinohara, M. Takeda and S. Arikawa, “Discovering best
variable-length-don’t-care patterns,” In Proc. 5th International Conference on Discovery Science LNCS 2534, 86–97 (2002).
13. M. Takeda, S. Inenaga, H. Bannai, A. Shinohara, and S. Arikawa, “Discovering Most
Classificatory Patterns for Very Expressive Pattern Classes”, In Proc. 6th International Conference on Discovery Science, LNCS 2843, 486–493 (2003).
14. P. Weiner, “Linear pattern matching algorithms,” In Proc. 14th IEEE Annual Symp.
on Switching and Automata Theory, 1–11 (1973).

March 3, 2004

0:3

WSPC/INSTRUCTION FILE

jbcb

16

15. E. Ukkonen, “On-line construction of suffix trees,” Algorithmica 14, 249–260 (1995).
16. D. Harel, and R. Tarjan, “Fast algorithms for finding nearest common ancestors,”
SIAM J. Comput. 13, 338–355 (1984).
17. M. A. Bender and M. Farach-Colton, “The LCA problem revisited,” In Proc. Latin
American Theoretical Informatics (LATIN) LNCS 1776, 88–94 (2000).
18. X. Leroy, D. Doligez, J. Garrigue and J. Vouillon, The Objective Caml system:
http://caml.inria.fr/ (1996–).
19. J. Zhu and M. Zhang, “SCPD: A promoter database of the yeast Saccharomyces
cerevisiae,” Bioinformatics 15, 607–11 (1999).
20. Y. Tamada, S. Kim, H. Bannai, S. Imoto, K. Tashiro, S. Kuhara and S. Miyano, “Estimating gene networks from gene expression data by combining Bayesian networkmodel with promoter element detection,” Bioinformatics, 19 (Supplement: Proceedings of the European Conference on Computational Biology 2003), ii227–ii236 (2003).
21. S. Imoto, T. Goto and S. Miyano, “Estimation of genetic networks and functional
structures between genes by using Bayesian networks and nonparametric regression,”
In Proc. Pacific Symposium on Biocomputing (PSB) 7, 175–186 (2002).