bioinformatics vol 19 no 15 2003 pages 1952 1963doi 10 1093 bioinformatics btg266
background rareness-based iterative multiple
sequence alignment algorithm for regulatory
element detection
chandrasegaran narasimhan philip locascio and edward
uberbacher∗
life sciences division oak ridge national laboratory po box 3480 oak ridge tn 37830 usa
received on november 8 2002 revised on march 13 2003 accepted on april 25 2003
abstract
motivation experimental methods capable of generating
sets of co-regulated genes have become commonplace how-
ever recognizing the regulatory motifs responsible for this
regulation remains difficult as a result computational detec-
tion of transcription factor binding sites in such data sets
has been an active area of research most approaches
have utilized either gibbs sampling or greedy strategies to
identify such elements in sets of sequences these exist-
ing methods have varying degrees of success depending on
the strength and length of the signals and the number of
available sequences we present a new deterministic iter-
ative algorithm for regulatory element detection based on a
markov chain background as in other methods sequences
in the entire genome and the training set are taken into
account in order to discriminate against commonly occurring
signals and produce patterns which are significant in the
training set results the results of the algorithm compare favorably with
existing tools on previously known and newly compiled data
sets the iteration based search appears rather rigorous not
only finding the binding sites but also showing how the bind-
ing site stands out from genomic background the approach
used to score the results is critical and a discussion of various
scoring schemes and options is also presented benchmark-
ing of several methods shows that while most tools are
good at detecting strong signals gibbs sampling algorithms
give inconsistent results when the regulatory element signal
becomes weak a markov chain based background model
alleviates the drawbacks of map maximum a posteriori log
likelihood scores availability available on request from the authors contact uberbacherec@ornl gov
supplementary information data and the results presented
in this paper are available on the web at http compbio ornl gov mira index html
∗to whom correspondence should be addressed introduction
with the annotation of many genomes the stage is set for
understanding the complex inter-relationships between genes
and the ways in which they act in concert in recent years with new experimental tools for studying gene expression more attention is being paid to the information present in
the non-coding regions of the genome particularly the reg-
ulatory regions gene regulation is one of the fundamental
mechanisms by which an organism responds to external and
internal biochemical signals and facilitates differentiation and
important developmental events the importance of understanding the mechanisms of
gene regulation has been well known in life sciences one of the theories for the cambrian explosion is that
the developments of new regulatory networks lead to
prototypes of virtually all-modern species innovative
experimental techniques like microarrays bucher 1999 http cmgm stanford edu pbrown sporulation chu et al 1998 derisi et al 1997 lockhart et al 1996 spellman
et al 1998 wolfsberg et al 1999 zhang 1999 and chip-
array experiments ren et al 2000 iyer et al 2001 lieb
et al 2001 simon et al 2001 have provided a flood of data
that can be mined for conserved regulatory motifs within a
species and between species methods to detect such elements
embedded in non-coding sequences are bound to accelerate
our understanding of transcription factors and the vital role
they play in complex organisms there have been two approaches for finding binding sites one approach is based on dictionary building or word count-
ing jensen and knudsen 2000 vanet et al 2000 van
helden et al 1998 2000 and another is based on finding
over-represented words within a data set by using mul-
tiple sequence alignments stormo and hartzell proposed
a greedy nearest neighbor algorithm stormo and hartzell 1989 hertz et al 1990 which tries to find putative motifs
on the basis of information content schneider et al 1986 which is based on information theory kullback and leibler 1951 consensus hertz and stormo 1999 is the most
1952 bioinformatics 19 15 oxford university press 2003 all rights reserved at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from algorithm for regulatory element detection
recent version of their greedy algorithm which can also estim-
ate the statistical significance of information content for a
given quantity of input data lawrence and reilly proposed an expectation maximization em algorithm lawrence and reilly 1990 which represents
a maximum likelihood method to solve the problem simil-
arly meme bailey and elkan 1994 1995 is a popular em
algorithm for aligning 5′ upstream regions for cis-regulatory
elements in 1993 lawrence and colleagues proposed a
gibbs sampling algorithm lawrence et al 1993 which is a
stochastic variant of the expectation maximization approach em and gibbs sampling algorithms have been successfully
used to search for binding motifs in gene expression data chu
et al 1998 http cmgm stanford edu pbrown sporulation spellman et al 1998 alignace hughes et al 2000 http atlas med harvard edu a gibbs sampling algorithm uses maximum a posteriori
log likelihood map scores to detect binding sites alig-
nace also incorporates two different features into the scoring
scheme to predict a motif they are group specificity a meas-
ure of how well a given motif is localized to certain upstream
regions and positional bias which measures the concentra-
tion of motifs within certain distance from transcriptional start
site in 2001 liu and colleagues proposed a new scoring scheme which is based on a multi-order markov chain background liu et al 2001a the tool called bioprospector is a gibbs
sampling-based multiple sequence alignment program thijs
et al proposed a higher order background model derived from
intergenic sequences thijs et al 2001 mdscan liu et al 2002 has been developed for binding site detection where a
confidence estimate can be applied to the upstream regions of
potentially co-regulated genes in expression array experiment all of the above methods have power for locating reg-
ulatory binding sites but further significant improvement
in methods is needed in this paper a new iterative
algorithm called mira markov-background iterative rare-
ness algorithm is presented and evaluated using existing
microarray data as well as data from the saccharomyces
cerevisiae promoter database zhu and zhang 1999 http cgsigma cshl org jian its performance is compared
with existing methods a discussion of advantages and
disadvantages of various scoring scheme is also presented the problem statement
from a given set of genes the goal is to find the subset of
genes that are likely to be co-regulated as well as the likely
binding site responsible for this observed co-regulation the
premise of the current work is that such subsequences are
over-represented in the upstream regions of the co-regulated
genes but generally rare in other parts of the genome this
makes sense from an evolutionary selection and a practical
standpoint since only a portion of the sequences in the set
may actually be co-regulated and some sequences may have
multiple sites for the same transcription factor motif detection
software must be able to accommodate these uncertainties given n sequences likely to contain the binding site the width
w for the site and r a rough estimate of number of sequences
likely to have the binding site mira returns top ten putative
binding sites the algorithm
the mira algorithm uses an iterative procedure to find a
group of ‘words’ or motifs that are over-represented in the
upstream sequence data sets yet are rare in the overall gen-
omic background the basic premise of the algorithm is that
transcription factor binding sites tend to occur exclusively in
the upstream regions of certain genes and are not likely to be
common in the remaining part of the genome such spurious
occurrences could lead to miss-regulation to achieve this
objective we use a scoring scheme termed ‘least likely
consensus’ llc which is composite of the information
content of a set of words or motifs found combined with the
background rareness br or logarithm of inverse background
frequency of the word or motif in the overall genome back-
ground the premise is that sets of words which are consistent
with each other in an upstream data set but have low frequency
in the overall genome are potentially biologically meaningful the algorithm chooses the first sequence as the root
sequence the first word with width w user input parameter is chosen and a candidate set of words is collected from rest
of the sequences based on a certain pre-defined mismatch dis-
tance the number of mismatches allowed is dependent on
w and is shown in table 1 if the mismatch distance is too
large then there would be more noise in the initial word set leading to possible omission of binding sites but it was also
observed that tolerance for mismatch increases with increas-
ing w this mismatch distance can also be defined by the user
as an optional input parameter a position specific transition
matrix pstm is constructed explained later in the article from all the candidate words the score for each word in this
initial set transition ratio or tr can be calculated as the ratio
of the log of pstm score for the word to that of background
model score this is given by
tr = log2 q b1 w
∑
i=2
log2 q bi bi−1 − log2 p b1 bm w
∑
l=m 1
log2 p bl bl−1 bl−m where q b1 is the probability of occurrence of nucleotide
b1 at position 1 in the foreground q bi bi−1 the trans-
ition probability to position i based on position i − 1 in the
foreground p b1 bm the probability of occurrence of the
1953 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from c narasimhan et al table 1 the corresponding mismatches allowed for the specified width
of the word width 7 8 9 10 11 12 13 14 15 16 17 18
mismatch 1 2 3 3 4 4 5 5 6 7 7 8
dna segment in the background and p bl bl−1 bl−m is
the transition probability based on previous m positions in the
background
low scoring words are removed from the initial candidate
set and pstm is recomputed on the remaining words until
convergence the candidate set is reduced 20 lowest scor-
ing words based on the trs to the expected number of binding
sites r based on user input the llc score for this final set
is given by
llc = log2 n w
∑
i=1
t
∑
j=a
pij log2 pij
−
1
n
n
∑
k=1 log2 pk b1 bm w
∑
l=m 1
log2 pk bl bl−1 bl−m where n is the number of words in the set pij the frequency
of nucleotide j at position i pk b1 bm the probability
of occurrence of kth word segment in the background and
pk bl bl−1 bl−m the transition probability of kth word
at position l in the background the first major term in the llc score above represents the
consistency of the word set how well they are represented as
a consensus while the latter terms represent the background
rareness of the word set information content ic and br are also calculated separ-
ately to visualize and study the individual contributions of ic
and br on the overall score for this application the ic for
the candidate set is given by
ic =
w
∑
i=1
t
∑
j=a
pij log2 pij
qj where pij is the frequency of nucleotide j at position i and
qj is the background frequency of nucleotide j the background rareness br or logarithm of inverse
background frequency above is given by
br = −
1
n
n
∑
k=1 log2 pk b1 bm w
∑
l=m 1
log2 pk bl bl−1 bl−m where n is the number of words in the set pk b1 bm the
probability of occurrence of kth word segment in the back-
ground pk bl bl−1 bl−m the transition probability of kth
word at position l in the background this procedure is repeated for all the words of width w in
the root sequence after this procedure one has an llc score
for each of the root word in the root sequence throughout
this analysis 1000 base pair upstream sequences were used thus one will have 1000−w llc ic and br scores at the
end of this step the average information content of the candidate word
sets which are in the top one-fifth in terms of llc scores
are computed the candidate word sets which have higher
information content than this average are re-seeded with the
consensus obtained from the first iteration as the root word after this second iteration the top ten word sets with the
highest llc scores is returned as putative binding sites the
program also returns all the root words with the correspond-
ing llc ic and br scores for generating ic versus br plots this data can also be used to study the rank of oligonucleotide
of interest to the biologist position specific transition matrix
for iteration
the algorithm’s iteration process reduces the initial set of t
words to a set of r words generally for typical data sets it
was observed that t is 10 20 times the value of r previous
methods represent the binding site as a w × 4 position spe-
cific frequency matrix psfm recently higher order profile
matrices have been shown to separate signal from noise and
give better results keich and pevzner 2002 mira uses first
order position specific transition matrix pstm of dimen-
sion w × 4 × 4 for the purpose of iterative convergence it
was observed that a first order pstm had much better con-
vergence behavior compared to a position-independent psfm
for micro array data sets especially where only a fraction of
sequences have a binding site for a given word set pstm stores all possible nucleotide
transitions between positions i − 1 and i the transition score t of a word using pstm is given by the expression
t = log2 q b1 w
∑
i=2
log2 q bi bi−1 where q b1 is the frequency of base b1 at position 1 and
q bi bi−1 is the frequency of transition from bi−1 to bi using a first order model builds dependencies in adja-
cent positions that can be poorly represented in a position-
independent model for example when the size of candidate
set is large compared to the number of binding sites in the data
set the most frequently observed nucleotides in neighboring
positions may not actually occur simultaneously in a binding
site with significant probability this is illustrated in figure 1 1954 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from algorithm for regulatory element detection
i 2 i
i-2 i-1 i 1
0 00
0 470 16 a 0 00
0 530 42 c 0 21
g 0 000 26 0 000 16 0 79 t 0 5
0 4 0 5
1 0 t g c a
t
g
c
a
t
g
c
a
t
g
c
a
t g c a
0 6 1 0
1 0
1 0
0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 1 00 0 00 0 00 0 60 0 40 0 00 0 00 aa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt 0 00 0 00 0 00 1 00 0 00 0 00 0 50 0 50 0 00 0 00 0 00 1 00 0 00 0 00 0 00 1 00 fig 1 psfm and pstm for a set of words from positions i − 1
to i 1 of a word set of width w the word set for this example
consists of 5-gta- 4-cta- 4-cgc- 3-atc- and 3-ttc- while the highest
scoring word according to psfm is ctc- the highest scoring ‘path’
according to pstm is gta- pstm performed better than psfm
in the data set under investigation for simplicity the tables do not
contain pseudo-counts for the purpose of pattern recognition pseudo-counts which reflect the general background behavior are added to
the known observations this can be understood as adding b
‘random’ words to the existing word set the first position can
be treated as a psfm and the pseudo-counts b words are
distributed proportional to the mononucleotide background
frequency from the whole genome see lawrence et al 1993 footnote 20 for all other transition positions pseudo-counts
are calculated as a fraction of the first order transition probab-
ilities calculated from whole genome for practical purposes choosing b =
√
r gave favorable results lawrence et al 1993 tr defined earlier is used to remove low scoring words
during the iteration tr is a combination of t as a foreground
model with br background model
unlike gene prediction tools which can rely on strong codon
biases to distinguish coding from non-coding regions fickett
and tung 1992 the foreground signals for regulatory motifs
are very weak as a result it becomes critical to accurately
model the background so that foreground signals will have
contrast traditionally a mononucleotide background model
has been used for contrasting binding sites from non-binding
sites markov chains and hmms have been effectively used
in a number of different applications in biological sequence
analysis xu et al 1994 burge and karlin 1997 salzberg
et al 1998 bioprospector liu et al 2001a uses a fixed
order background model with variable order in the initial
chain positions within the bayesian framework in biopro-
spector the background score is computed from the following
expression
pbr = − log2 p b1 m
∑
l=2
log2 p bl b1 bl−1 w
∑
l=m 1
log2 p bl bl−1 bl−m thus a word of width w is represented by a markov chain
of order m with orders from 0 to m used for the initial
m 1 positions of the chain in 2002 mdscan liu et al 2001a utilized a background model suggested by thijs et al 2001 the background model is constructed from intergenic
sequences and computes the background probability of a word
of width w by calculating the transition probabilities at all pos-
itions with a single markov chain of order m for calculating
the transition probability at first m positions m bases preced-
ing the word is used thus the background score for a word
computed with a markov chain order m is given by
pbr = − w
∑
l=1
log2 p bl bl−1 bl−m mira uses a background-scoring scheme where the first m
positions are computed based on the frequency of that m-mer
in the background and the transitions that follow are pre-
dicted from an m-order chain this makes the background
score independent of the surrounding dna the expression
for background score is
pbr = − log2 p b1 bm w
∑
l=m 1
log2 p bl bl−1 bl−m it should be noted that as m tends toward w the score pbr tends to change from an estimate to the actual number of
occurrences of that segment in the background using the
1955 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from c narasimhan et al expression for pbr the value of br can be calculated from the
expression
br = −
1
n
n
∑
k=1 log2 pk b1 bm w
∑
l=m 1
log2 pk bl bl−1 bl−m br represents the normalized background frequency of occur-
rence estimate of a collection of n segments the normaliza-
tion is required for comparing the br scores of two different
word sets that contain different number of segments to verify previous studies and to determine an effective
background model order for mira markov chains from order
one to order six were developed and tested using two data
sets the yeast whole genome sequence and yeast intergenic
sequence the intergenic background model performed the
best with a 4th order model optimal for longer binding sites and a 6th order model optimal for short binding sites w 9 the br score measures a quantity similar in character to
group specificity score s hughes et al 2000 this work
incorporates group specificity measurements into the scoring
scheme for alignace this measurement with a threshold of
s 10−10 was used to identify 25 known motifs in the yeast
genome and was shown to be effective for identifying yeast
transcription factor binding sites the group specificity score
is given by
s =
min s1 s2 ∑
i=x s1
i m − s1
s2 − i m
s2 where m is the total number of orfs in the yeast genome s1
is the set of orfs in which the motif is present s2 is the set
of orfs given as input to the software s is the probability of
finding x or greater number of motifs out of s2 sequences as
the number of binding sites in the input sequences becomes
comparable to the number of motifs present in the whole gen-
ome the overlap of sets s1 and s2 increases as well and the
numerator in the expression approaches 1 and s approaches 0 likewise if a motif is present only in the input data set then
it will have the highest possible br compared to any other
collection of words from the data set data sets
the yeast genome saccharomyces cerevisiae was used as
the basis for the current study because of the availability of
expression and regulatory motif data sets related to different
cellular processes such as cell cycle sporulation etc some
data sets have been derived from the collection of publica-
tions in pubmed http www ncbi nlm nih gov and from
individual web sites maintained by researchers some data
sets abf1 ebf1 reb1 & mcm1 were extracted from
the s cerevisiae promoter database scpd zhu and zhang 1999 the data set for pdr was obtained from the alignace
web server hughes et al 2000 the data set for gal4 is a
combination of scpd data and alignace data data sets for
csre kratzer and schuller 1997 roth and schuller 2001 walther and schuller 2001 and hap1 deckert et al 1995 kennedy et al 1999 lodi et al 1999 is a combination of
scpd data and manually compiled data from recent publica-
tions in pubmed the mig1 data set was extracted from the
table compiled by christopher klein klein et al 1998 http mips gsf de proj yeast cygd db index html the
rap1p data was compiled from pavlidis et al 2001 http www cs columbia edu compbio prom-svm and lascaris
et al 1999 the data set for metabolic phase in sporula-
tion is compiled from the work of chu chu et al 1998 http cmgm stanford edu pbrown sporulation the fkh1
microarray zhu et al 2000 data set was compiled from
the clb2 cluster listed in spellman et al 1998 the fkh1 pic et al 2000 data was compiled from supplementary
website of the work of zhu et al 2000 http genome-
www stanford edu fkh fkhtab1 htm the mac1 data set was
extracted from a microarray study by gross et al 2000 for
all the genes in each cluster group 1000 base pair upstream
sequence was extracted and used as input for in silico analysis results
the stability of mira was tested on biological sequences by
gradually increasing the noise in the system to a factor of 2 a comparison of the performance of mira with other search
algorithms was also performed these results are presented in
this section the present study concludes that when 1000 base
pair long upstream sequences are used at least five sequences
actually containing binding sites are required for correct iden-
tification of regulatory motifs when the number of 1000 base
pair sequences with motifs is increased to 10 the prediction
method became less susceptible to noise the possible reasons
for limitations are also presented and discussed specificity and sensitivity analysis
to check the performance of mira on noisy data sets and
to see the dependence of results on the input parameter r two different sets of runs were conducted two curated data
sets for four different binding sites with r = 5 and r = 10
were made the binding sites were chosen in such a way as to
represent real biological situation both short and long motifs as well as very conserved fkh1 motifs and motifs with very
variable positions abf1 were included in this analysis an
increasing number of randomly picked upstream sequences 1000 base pairs long were added to these data sets the results
of mira on these data sets are shown in tables 2 and 3 and
figure 2 specificity is defined as the number of sequences
identified with the correct consensus divided by the number
1956 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from algorithm for regulatory element detection
table 2 results of specificity and sensitivity analysis for five sequences with motif motif name n = 5 r = 5 n = 6 r = 5 n = 7 r = 5 n = 8 r = 5 n = 9 r = 5 n = 10 r = 5
fkh1 rank 4 10 w = 7 sensitivity 5 5 5 5
specificity 5 5 5 5
abf1 rank 5 3 6 w = 12 sensitivity 5 5 5 5 5 5
specificity 5 5 5 5 5 5
rap1p rank 1 1 1 2 4 3 w = 14 sensitivity 5 5 5 5 5 5 4 5 4 5 4 5
specificity 5 5 5 5 5 5 4 5 4 5 4 5
mcm1 rank 4 5 5 6 5 w = 16 sensitivity 5 5 5 5 4 5 4 5 4 5
specificity 5 5 5 5 4 5 4 5 4 5
the table shows the rank along with specificity and sensitivity for four motifs with increasing noise in the data set for all the results shown in this table five sequences in each data
set had the known motifs r = 5 dashes indicate that the motif was not ranked in top 10 table 3 results of specificity and sensitivity analysis for ten sequences with motif
motif name n = 10 r = 10 n = 12 r = 10 n = 14 r = 10 n = 16 r = 10 n = 18 r = 10 n = 20 r = 10
fkh1 rank 1 1 2 1 2 5 w = 7 sensitivity 10 10 10 10 10 10 10 10 10 10 10 10
specificity 10 10 10 11 10 11 10 13 10 13 10 13
abf1 rank 2 1 2 3 5 w = 12 sensitivity 10 10 10 10 9 10 9 10 8 10
specificity 10 10 10 10 9 10 9 10 8 10
rap1p rank 1 1 2 1 5 1 w = 14 sensitivity 10 10 10 10 9 10 9 10 8 10 9 10
specificity 10 10 10 10 9 10 9 10 8 10 9 10
mcm1 rank 1 2 2 1 1 1 w = 16 sensitivity 10 10 9 10 9 10 9 10 9 10 10 10
specificity 10 10 9 10 9 10 9 10 9 10 10 12
the table shows the rank along with specificity and sensitivity for four motifs with increasing noise in the dataset for all the results shown in this table five sequences in each data
set had the known motifs r = 10 dashes indicate that the motif was not ranked in top 10 of sequences that were identified sensitivity is defined as the
number of sequences identified with the correct consensus
divided by the number of sequences which actually had the
binding site the results in table 2 indicate that for short binding sites the quality of prediction decreases rapidly with increasing
level of noise in the data sets on the other hand the results
obtained for long binding sites were quite satisfactory even
when the noise in the data was doubled n = 10 r = 5 on further analysis it was observed that the possibility for
obtaining ic higher than that of the binding site increased with
increasing noise in the data sets this was especially true for
short motifs for example all the top 10 ranked consensus for
fkh1 n = 10 r = 5 data set had a higher ic than the known
motif this explains the decreasing quality of prediction with
added noise if the number of binding sites in the data sets
were increased to 10 the overall prediction quality improved see table 3 yet there was a general decrease in quality of
solution in terms of rank specificity and sensitivity of the
motif compared to the situation with less noise the overall
results of tables 2 and 3 have been summarized as a graph in
figure 2 tables 2 and 3 show that the quality of prediction depends
on the number of binding sites in the data set investigation of
microarray data sets showed that mira grew more tolerant
of improper specification of r with increasing number of
binding sites in the data set for a completely new data set a user can run mira with four different values of r r =
60 70 80 and 90 of n and see the variation in results if
a putative motif occurs consistently in some of these runs 1957 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from c narasimhan et al fig 2 the figure shows the average rank of the motifs from table 3 when r was set to 5 and table 4 when r was set to 10 for all
the cases when the motif was not ranked in the top 10 the rank was
taken as 15 the figure shows that the prediction quality improved
with increase in number of sequences with binding sites then they are the potential binding sites which occur in that
fraction of input sequences width of the motif
the current algorithm requires width of the motif to be spe-
cified this limits the complexity and dimensionality of the
search investigation of the effects generated by varying the
width on curated data sets gave favorable results however the quality of prediction decreased when there was significant
noise in the data sets results show that in general w can be
varied ±3 without affecting the solution figure 3 shows a
plot of ics versus the corresponding width from the data sets
of all the motifs used in this paper as the width of a binding
site increases there tend to be more positions that are variable even though a 17 mer has an ic of 19 0 the corresponding
ic scaled to 7 will be 7 82 which is far lower than the value
of 13 0 which is the ic of a binding motif of corresponding
width for an ic of 7 82 for a 7 mer the specificity of scoring
system decreases rapidly because many random word sets which are not binding sites can generate such a score see
table 2 fkh1 this explains the difficulty in making the
search independent of w most other methods also require w
to be specified and the current study indicates that this seems
to be general limitation comparison with other methods
the current study performed analysis on a series of published
microarray data sets as well as on well-documented binding
sites from scpd analysis was performed on 14 different
fig 3 the figure shows the plot of ic versus the corresponding
width of the 13 motifs used in the current study the dotted line is a
linear approximation of the plot data sets and the results were compared with existing motif
detection tools the results of these analyses are presented in
table 4 for the current purpose a tool was considered to have
found a binding site if it recognized at least half of the site and
ranked this within the top 10 scoring sets all methods were
configured to detect motifs on both strands and report at least
10 motifs mdscan requires at least 20 sequences with bind-
ing sites for convergence and so could be only run on four
data sets for all mdscan runs the maximum possible set
of sequences was used during the search step up to the max-
imum of 30 allowed by the tool for gibbs sampling based
tools three independent runs were performed to evaluate the
consistency of the results because the results generated by
these algorithms depend on the seed to the random number
generator in spite of the difficult nature of the problem the results
reflect a fair amount of consistency in the prediction quality of
all the tools the quality or the strength of prediction increased
with increase in both the number of sequences containing the
motif and the width of the motif mira found 14 of the 14 sites tested consensus
recognized 10 of the 14 sites the results generated by con-
sensus however whenever it found the site were in the top
two the results generated by gibbs sampling tools biopro-
spector and alignace showed a certain amount of variability
based on the random seed alignace located approximately 11
of 14 motifs bioprospector with whole genome background
model recognized 13 of the 14 sites for these tools and cer-
tain data sets different consensus sequences were obtained
for certain recognized motifs compared to other methods the consistency of gibbs sampling algorithms decreases more
1958 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from algorithm for regulatory element detection
table 4 results
motif r n a known consensus
versus
predicted consensus
by mira
mira consensus align-
ace
biopro-
spector whole
genome biopro-
spector inter-
genic mds can intergenic abf1 20 20 cgtnnnnnntga
cgtnnnnnntga 1 1 3 3 4 1 1 1 1 2 1 1
csre 7 7 tccattcatccg
tccattcatccg 1 1 2 2 2 1 1 1 1 1 1 na
ebf1 16 16 acaccca
acaccca 2 1 5 4 6 — — 1 3 1 na
fkh1 42 42 gtaaaca
gtaaaca 1 — — 1 1 1 5 2 1 1
fkh1 micro-array clb2 18 32 gtaaaca
gtaaaca 1 — — 2 3 2 6 5 10 na
gal4 10 10 cggnnnnnnnnnnnccg
cgggcgacagtgatccg 1 1 1 1 1 1 1 1 1 1 1 na
hap1 8 8 ccgntannnccg
accgata tttcc 1 — — 5 10 9 — — — na
mac1 6 6 tttgctca
tttgctca 1 1 1 1 — 10 2 2 na
mcm1 24 24 tttcccaaannggaaa
tttcctaattaggtaa 1 1 2 2 2 1 1 1 1 1 1 1
met-spr urs1 20 47 tagccgccga
tagccgccga 2 1 2 3 3 4 3 2 9 1 1 1 na
mig1 11 11 aaaaatctgggg
aaaatacgggga 2 1 3 3 3 1 1 1 1 2 1 na
pdr 11 11 tccgcgga
tccgcgga 2 2 10 8 2 1 1 1 2 4 1 na
rap1p 44 44 tacacccatacatt
ttacacccatacat 1 1 1 1 1 1 1 1 1 1 1 1
reb1 13 13 ttacccg
ttacccg 1 1 6 9 1 1 1 1 2 1 na
table shows the various data sets along with the known binding sites and the comparison of ranks of various algorithms an ‘—’ indicates that the motif was not ranked in the top
ten putative sites by that algorithm for gibbs sampling based algorithms 3 or more runs were performed to check the consistency of the results an is the total number of sequences in the data set r is the estimate of number of sequences likely to have the binding site dramatically as the strength of the motif decreases for
example these methods miss the hap1 and reb1 sites in
the gibbs sampling tools the score obtained for the same
motif from two different runs can be different liu 2001b for example bioprospector with motifs hap1 and mac1 and alignace with ebf1 mdscan which can only be run
on sets of 20 or more sequences found 4 of 4 sites in these
data sets another interesting point about table 4 is the result
obtained for the gal4 data set according to scpd zhu and zhang 1999 gal4 binds to a consensus of
cggnnnnnnnnnnnccg hence according to table 1 mira should not be able to find the binding site because
the number of mismatches allowed for a 17 mer is 7 and
the mismatch required for gal4 is 11 the latest findings
reported from analyzing 23 sequences obtained from a chip-
array experiment ren et al 2000 shows the consensus to
be cggassaststsstccg liu et al 2002 this explains
the positive mira result for gal4 mira has been tested on
many other yeast data sets besides the ones shown in table 4 and virtually all of these met the mismatch condition discussion
regulatory element detection tools typically have two com-
ponents the first component is the alignment searching
algorithm that collects and aligns similar motifs together the second is a scoring scheme for such sets of the aligned
motifs for n sequences of length l an o ln algorithm is
required to explore all possible alignments and this is compu-
tationally intractable stormo and hartzell 1989 therefore
binding site detection tools make approximations in the search
process for example consensus achieves efficiency by
only keeping a limited number of motif matrices as machine
memory will allow gibbs samplers sample the search space
to converge to binding sites the current algorithm mira 1959 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from c narasimhan et al is a semi-heuristic semi-iterative algorithm with a predefined
neighbor range defined in table 1 unlike the other methods
in table 2 which are fast linear time algorithms mira is
an iterative search over all sequences within a running time
of o wnl2 the sensitivity of mira is partly derived
from the extensive searching scheme built into it unlike
mdscan mira is intended to work on microarray data
sets and does not require a ‘confidence set’ for search-
ing this is yet another reason why the search needs to be
performed over all sequences mira does require a root
sequence that has the binding motif it may fail to detect
the binding site if the root sequence does not contain the
motif researchers can get around this problem by choos-
ing the upstream region of a specific gene of their interest
as the root sequence or perhaps testing several sequences as
the root the second component the scoring scheme seems to be
a significant drawback of many regulatory element detection
tools the llc score used by mira is similar to the one
presented in liu et al 2002 and has two components ic and
br to provide a means to examine the individual effects of
these two components in separating signal from noise mira
generates ic and br values for all word sets after the first
phase of iteration the binding sites generally have both high
ic and br scores and generally fall in the top right-hand
corner of the plots of ic versus br figs 4 6 the initial
premise that binding sites in upstream regions are rare else-
where in the genome seems to be justified and the br factor
significantly helps the scoring in mira compared to the use
of ic alone the br scores of binding sites perhaps reflect the
selective pressure against occurrences of these motifs in other
regions of a genome figure 4 shows ic versus br plots for the hap1 data set and for both first order and fourth order markov chain back-
ground models for the fourth order plot the points are more
dispersed in space compared to the first order plot even though
the location of binding site remains unchanged this indicates
that the fourth order model is more effective at discriminating
real from false signals these plots also demonstrate several other points from
the nature of the cluster one can get a sense of how closely
related the input sequences are evolutionarily close sequences
would land near each other on the plot they also provide
a means for visualizing the relative strength of a prediction
and how well true motifs contrast against false positives and
background noise by comparing these results with biopro-
spector and alignace both of which are gibbs sampling based
algorithms one can say that the markov chain based scor-
ing scheme improves the results presented by these previous
works liu et al 2001a thijs et al 2001 hughes et al 2000 state that the ‘drawback of map scores
is that some motifs occurring ubiquitously in a genome a rich
regions are scored very highly but are not likely to be relev-
ant to specific set of genes in question’ figure 5 illustrates an a b fig 4 plot of ic versus br for all the word sets for the hap1 data
set using first order a and fourth order b markov backgrounds one can see that the points in the graph are more dispersed in b than in a the binding site stands out more in figure b example of such a case and shows ic versus br plots for the
met-spr data set in table 2 genes active in the metabolic
phase of sporulation in yeast there are clearly many word
sets representing common words having low br scores and
which have higher ic than that of the binding site but a scor-
ing scheme that includes both of these factors can overcome
the ‘common word’ problem while both ic and br are important for scoring potential
binding site sets ic is less adequate for shorter binding sites fig 6 this behavior is to be expected since with a shorter
word size there is a higher chance of ‘accidental consensus’
among upstream sequences and therefore a higher likelihood
1960 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from algorithm for regulatory element detection
fig 5 ic versus br plot for sporulation data set it is interesting
to note that there are four word sets which have higher ic than the
motif yet another interesting feature is the shape of the cluster this
probably reflects the composition of sequences in the data set fig 6 ic versus br plot for ebf1 data set note that there are a lot
of word sets with higher ic than the binding site of false positives if ranking is based wholly on information
content scores in the current study over-represented words
occurring ubiquitously in the upstream regions are penalized
using their br score the br model also helps longer motifs
appear over-represented in upstream regions results in table 4 show that a simple mononucleotide back-
ground model gives reasonable results consensus and
alignace use this model it is interesting to note the data sets
for which the mononucleotide background model fails for
instance the fkh 1 binding site consensus is gtaaaca
or tgtttac reverse complement thus 5 out of 7 posi-
tions are either an a or a t in a at rich genome like yeast 62 the mononucleotide background scoring for a sequence
like gtaaaca is likely to estimate this as being too com-
mon to represent a real binding site however algorithms like
bioprospector mdscan and mira that use more complex
background models find this site the basic premise on which the br model is built that an
occurrence of transcription factor binding sites in the genome
is tightly controlled by natural selection may not be valid in
all situations if a transcription factor usually requires other
factors to facilitate transcription cooperative binding factors the selective pressure against a single motif occurring through-
out the genome could be diminished compared to situations
where a factor acts more independently regulatory element detection continues to be a challenging
problem for computational biologists for many reasons some
of which were described above despite the limitations the
comparison across existing tools represented in table 2 indic-
ates the general usefulness of these methods in predicting
regulatory elements further improvements in the perform-
ance will most likely require more complex background model
treatment as well as an incorporation of potential interactions
among sites and factors this will be especially true as more
complex eukaryotes with even more complicated regulatory
networks are examined acknowledgements
the author would like to thank dr xiaole shirley liu
at stanford medical informatics stanford university for
providing yeast intergenic data and some references the
author would also like to thank colleagues dr ying xu dr victor olman dr loren hauser dr frank larimer
and dr jay snoddy at ornl this research was sponsored
by the office of biological and environmental research at
oak ridge national laboratory managed by ut-battelle llc for the us department of energy under contract no
de-ac05 00or33735 references
bailey t l and elkan c 1995 unsupervised learning of multiple
motifs in biopolymers using expectation maximization mach learn 21 51 80 bailey t l and elkan c p 1994 fitting a mixture model by expect-
ation maximization to discover motifs in biopolymers intelligent
syst mol biol 2 28 36 bucher p 1999 regulatory elements and expression profiles curr opin struct biol 9 400 407 burge c and karlin s 1997 prediction of complete gene structures
in human genomic dna j mol biol 268 78 94 1961 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from c narasimhan et al chu s derisi j eisen m b mulholland j botstein d brown p o and herskowitz i 1998 the transcriptional program
of sporulation in budding yeast science 282 699 705 deckert j perini r balasubramanian b and zitomer r s 1995 multiple elements and auto-repression regulate rox1 a repressor
of hypoxic genes in saccharomyces cerevisiae genetics 139 1149 1158 derisi j l iyer v r and brown p o 1997 exploring the meta-
bolic and genetic control of gene expression on a genomic scale science 278 680 686 fickett j w and tung c s 1992 assessment of protein coding
measures nucleic acids res 20 6441 6450 gross c kelleher m iyer v r brown p o and winge d r 2000 identification of the copper regulon in saccharomyces
cerevisiae by dna microarrays j biol chem 275 32310 32316 hertz g z hartzell g w iii and stormo g d 1990 identification
of consensus patterns in unaligned dna sequences known to be
functionally related comput appl biosci 6 81 92 hertz g z and stormo g d 1999 identifying dna and pro-
tein patterns with statistically significant alignments of multiple
sequences bioinformatics 15 563 577 hughes j d estep preston w tavazoie s and church g m 2000 computational identification of cis-regulatory elements
associated with groups of functionally related genes in saccharo-
myces cerevisiae j mol biol 296 1205 1214 iyer v r horak c e scafe c s botstein d snyder m and
brown p o 2001 genomic binding sites of the yeast cell-cycle
transcription factors sbf and mbf nature 409 533 538 jensen l j and knudsen s 2000 automatic discovery of reg-
ulatory patterns in promoter regions based on whole cell
expression data and functional annotation bioinformatics 16 326 333 keich u and pevzner p a 2002 finding motifs in the twilight zone bioinformatics 18 1374 1381 kennedy m a barbuch r and bard m 1999 transcriptional
regulation of the squalene synthase gene erg9 in the yeast
saccharomyces cerevisiae biochim biophys acta 1445 110 122 klein c j olsson l and nielsen j 1998 glucose control in sac-
charomyces cerevisiae the role of mig1 in metabolic functions microbiology 144 13 24 kratzer s and schuller h j 1997 transcriptional control of the
yeast acetyl-coa synthetase gene acs1 by the positive regulat-
ors cat8 and adr1 and the pleiotropic repressor ume6 mol microbiol 26 631 641 kullback s and leibler r a 1951 on information and sufficiency ann math stat 22 79 86 lascaris r f mager w h and planta r j 1999 dna-binding
requirements of the yeast protein rap1p as selected in silico from
ribosomal protein gene promoter sequences bioinformatics 15 267 277 lawrence c e and reilly a a 1990 an expectation maximiza-
tion em algorithm for the identification and characterization of
common sites in unaligned biopolymer sequences proteins 7 41 51 lieb j d liu x botstein d and brown p o 2001 promoter-
specific binding of rap1 revealed by genome-wide maps of
protein dna association nat genet 28 327 334 liu x brutlag d l and liu j s 2001a bioprospector discov-
ering conserved dna motifs in upstream regulatory regions of
co-expressed genes pac symp biocomput 6 127 138 liu j s 2001b monte carlo strategies in scientific computing springer new york liu x s brutlag d l and liu j s 2002 an algorithm for find-
ing protein dna binding sites with applications to chromatin-
immunoprecipitation microarray experiments nat biotechnol 20 835 839 lockhart d j dong h byrne m c follettie m t gallo m v chee m s mittmann m wang c kobayashi m horton h and brown e l 1996 expression monitoring by hybridiza-
tion to high-density oligonucleotide arrays nat biotech 14 1675 1680 lodi t alberti a guiard b and ferrero i 1999 regulation of
the saccharomyces cerevisiae dld1 gene encoding the mitochon-
drial protein d-lactate ferricytochrome c oxidoreductase by hap1
and hap2 3 4 5 mol gen genet 262 623 632 pavlidis p furey t s liberto m haussler d and grundy w n 2001 promoter region-based classification of genes pac symp biocomput 151 163 pic a lim f l ross s j veal e a johnson a l sultan m r west a g johnston l h sharrocks a d and morgan b a 2000 the forkhead protein fkh2 is a component of the yeast
cell cycle transcription factor sff embo j 19 3750 3761 ren b robert f wyrick j j aparicio o jennings e g simon i zeitlinger j schreiber j hannett n kanin e et al 2000 genome-wide location and function of dna binding
proteins science 290 2306 2309 roth s and schuller h j 2001 cat8 and sip4 mediate regulated
transcriptional activation of the yeast malate dehydrogenase gene
mdh2 by three carbon source responsive promoter elements yeast 18 151 162 salzberg s delcher a kasif s and white o 1998 microbial
gene identification using interpolated markov models nucleic
acids res 26 544 548 schneider t d stormo g d gold l and ehrenfeucht a 1986 information content of binding sites on nucleotide sequences j mol biol 188 415 431 simon i barnett j hannett n harbison c t rinaldi n j volkert t l wyrick j j zeitlinger j gifford d k jaakkola t s et al 2001 serial regulation of transcriptional
regulators in the yeast cell cycle cell 106 697 708 spellman p t sherlock g zhang m q iyer v r anders k eisen m b brown p o botstein d and futcher b 1998 com-
prehensive identification of cell cycle-regulated genes of the yeast
saccharomyces cerevisiae by microarray hybridisation mol biol cell 9 3273 3297 stormo g d and hartzell g w iii 1989 identifying protein-
binding sites from unaligned dna fragments proc natl acad sci usa 86 1183 1187 thijs g lescot m marchal k rombauts s de moor b rouze p and moreau y 2001 a higher-order background model
improves the detection of promoter regulatory elements by gibbs
sampling bioinformatics 17 1113 1122 van helden j andr e b and collado-vides l 1998 extracting
regulatory sites from upstream region of yeast genes by computa-
tional analysis of oligonucleotide frequencies j mol biol 281 827 842 1962 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from algorithm for regulatory element detection
van helden j rios a f and collado-vides j 2000 dis-
covering regulatory elements in non-coding sequences
by analysis of spaced dyads nucleic acids res 28 1808 1818 vanet a marsan l labigne a and sagot m f 2000 inferring
regulatory elements from a whole genome an analysis of helico-
bacter pylori sigma 80 family of promoter signals j mol biol 297 335 353 walther k and schuller h j 2001 adr1 and cat8 synergisti-
cally activate the glucose regulated alcohol dehydrogenase gene
adh2 of the yeast saccharomyces cerevisiae microbiology 147 2037 2044 wolfsberg t g gabrielian a e campbell m j cho r j spouge j l and landsman d 1999 candidate regulatory
sequence elements for cell cycle-dependent transcription in
saccharomyces cerevisiae genome res 9 775 792 xu y einstein j r mural r j shah m and uberbacher e c 1994 an improved system for exon recognition and gene mod-
eling in human dna sequences proc int conf intell syst mol biol 2 376 384 zhang m q 1999 large-scale gene expression data analysis a new
challenge to computational biologist genome res 9 681 688 zhu g spellman p t volpe t brown p o botstein d davis t n and futcher b 2000 two yeast forkhead genes
regulate the cell cycle and pseudohyphal growth nature 406 90 94 zhu j and zhang m q 1999 a promoter database of yeast
saccharomyces cerevisiae bioinformatics 15 607 611 1963 at m
asaryk university on february 22 2012
http bioinform
atics oxfordjournals org d
ow
nloaded from