PMS5: an efficient exact algorithm for the (ℓ, d)motif finding problem
 Hieu Dinh^{1},
 Sanguthevar Rajasekaran^{1}Email author and
 Vamsi K Kundeti^{1}
DOI: 10.1186/1471210512410
© Dinh et al; licensee BioMed Central Ltd. 2011
Received: 11 March 2011
Accepted: 24 October 2011
Published: 24 October 2011
Abstract
Background
Motifs are patterns found in biological sequences that are vital for understanding gene function, human disease, drug design, etc. They are helpful in finding transcriptional regulatory elements, transcription factor binding sites, and so on. As a result, the problem of identifying motifs is very crucial in biology.
Results
Many facets of the motif search problem have been identified in the literature. One of them is (ℓ, d)motif search (or Planted Motif Search (PMS)). The PMS problem has been well investigated and shown to be NPhard. Any algorithm for PMS that always finds all the (ℓ, d)motifs on a given input set is called an exact algorithm. In this paper we focus on exact algorithms only. All the known exact algorithms for PMS take exponential time in some of the underlying parameters in the worst case scenario. But it does not mean that we cannot design exact algorithms for solving practical instances within a reasonable amount of time. In this paper, we propose a fast algorithm that can solve the wellknown challenging instances of PMS: (21, 8) and (23, 9). No prior exact algorithm could solve these instances. In particular, our proposed algorithm takes about 10 hours on the challenging instance (21, 8) and about 54 hours on the challenging instance (23, 9). The algorithm has been run on a single 2.4GHz PC with 3GB RAM. The implementation of PMS5 is freely available on the web at http://www.pms.engr.uconn.edu/downloads/PMS5.zip.
Conclusions
We present an efficient algorithm PMS5 that uses some novel ideas and combines them with wellknown algorithm PMS1 and PMSPrune. PMS5 can tackle the large challenging instances (21, 8) and (23, 9). Therefore, we hope that PMS5 will help biologists discover longer motifs in the futures.
1 Background
The discovery of patterns in DNA, RNA, and protein sequences has led to the solution of many vital biological problems. For instance, the identification of patterns in nucleic acid sequences has resulted in the determination of open reading frames, identification of promoter elements of genes, identification of intron/exon splicing sites, identification of SH RNAs, location of RNA degradation signals, identification of alternative splicing sites, etc. In protein sequences, patterns have proven to be extremely helpful in domain identification, location of protease cleavage sites, identification of signal peptides, protein interactions, determination of protein degradation elements, identification of protein trafficking elements, discovery of short functional motifs, etc. Motifs are patterns found in biological sequences that are vital for understanding many biological subjects like gene function, human disease, drug design etc. As a result, the identification of motifs plays an important role in biological studies. The motif search problem has been attracting many researchers. In the literature, many versions of the motif search problem have been enumerated. Examples include Simple Motif Search (SMS), Planted Motif Search (PMS)  also known as (ℓ, d)motif search, and Editdistancebased Motif Search (EMS) (see e.g., [1]). In this paper, we will focus on the PMS problem (or PMS for short).
The definition of Planted Motif Search (PMS)
PMS is stated as follows. It takes as input n sequences, two integers ℓ and d. For simplicity, we assume that the length of each sequence is m. The problem is to identify all strings M of length ℓ such that M occurs in each of the n sequences with at most d mismatches. Formally, string M has to satisfy the following constraint: there exists a string M_{ i }of length l in sequence i, for every i (1 ≤ i ≤ n), such that the number of mismatches between M and M_{ i }is less than or equal to d. The number of mismatches between two strings of equal length is known as the Hamming distance between them. String M is called a motif. For example, if the input sequences are GCGCGAT, CAGGTGA, and CGATGCC; ℓ = 3; and d = 1, then GAT and GTG are some of the (l, d)motifs. PMS is a wellstudied problem and it has been shown to be NPhard. As a result, all known exact algorithms for PMS take exponential time in some of the underlying parameters in the worst case. Two kinds of algorithms have been proposed in the literature for PMS: exact and approximate. While an exact algorithm always finds all the motifs, an approximate algorithm may not always find all the motifs. Typically, approximate algorithms tend to be faster than exact algorithms. Some example approximate algorithms are due to Bailey and Elkan [2], Buhler and Tompa [3], Lawrence et al. [4], Pevzner and Sze [5], and Rocke and Tompa [6]. These algorithms are based on local search techniques such as Gibbs sampling, expectation optimization, etc. The WINNOWER algorithm of [5] is based on finding cliques in a graph. Some other approximate algorithms are: PROJECTION [3], MULTIPROFILER [7], PatternBranching [8], CONSENSUS [9], GibbsDNA [4], MEME [2], and ProfileBranching [8].
Although approximate algorithms are acceptable in some cases in practice, exact algorithms are preferable since they are guaranteed to report all the (l, d)motifs. For biologists, the motifs found by an algorithm could be much more important than its run time. As a result, we focus in this paper on efficient exact algorithms. Some exact algorithms known for PMS are: [10–18], and [19].
Buhler and Tompa [3] have employed PMS algorithms to find known transcriptional regulatory elements upstream of several eukaryotic genes. In particular, they have used orthologous sequences from different organisms upstream of four different genes: preproinsulin, dihydrofolate reductase (DHFR), metallothioneins, and cfos. These sequences are known to contain binding sites for specific transcription factors. The authors point out the differences between experimental data and synthetic data that PMS algorithms are typically tested with. For example, the background DNA in experimental data is not random. Their algorithm successfully identified the experimentally determined transcription factor binding sites. They have used the values of l = 20 and d = 2. The same sites have also been found using our PMS2 algorithm [11]. The algorithm of [3] is an approximation algorithm whereas PMS2 is an exact algorithm. Buhler and Tompa have also employed their algorithm to solve the ribosome binding site problem for various prokaryotes [3]. This problem is even more challenging since here the number of input sequences could be in the thousands.
Eskin and Pevzner [13] used PMS algorithms to find composite regulatory patterns. They point out that traditional pattern finding techniques (on unaligned DNA sequences) concentrate on identifying highscoring monads. A regulatory pattern could indeed be a combination of multiple and possibly weak monads. They employ MITRA (a PMS algorithm) to locate regulatory patterns of this kind. The algorithm is demonstrated to perform well on both synthetic and experimental data sets. For example, they have employed the upstream regions involved in purine metabolism from three Pyrococcus genomes. They have also tested their algorithm on four sets of S.cerevisiae genes which are regulated by two transcription factors such that the transcription factor binding sites occur near each other. Price and Pevzner [8] have employed their PatternBranching PMS technique on a sample containing CRP binding sites in E.coli, upstream regions of many organisms of the eukaryotic genes: preproinsulin, DHFR, metallothionein, & cfos, and a sample of promoter regions from yeast. They report finding numerous motifs in these sequences.
The performance of an exact algorithm is typically evaluated on random benchmark data generated as follows. Twenty input DNA sequences, each of length 600, are generated randomly from the alphabet Σ = {A, C, G, T}. A motif M of length ℓ is also generated randomly and planted in each of the input sequences within a Hamming distance of d to make sure that there always exists a motif in the input. Based on the values of ℓ and d, certain instances of PMS have been identified to be challenging. An instance is challenging if the expected number of the motifs that occur by random chance (in addition to the planted one) is one or more. For example, the following instances are challenging: (9, 2), (11, 3), (13, 4), (15, 5), (17, 6), (19, 7), (21,8), (23, 9), etc.
To compare the performance of exact algorithms, the challenging instances are commonly used. For example, the exact algorithm MITRA of [8] can solve the challenging instances (9, 2), (11, 3), and (13, 4). It takes either too much time or too much memory on the challenging instance (15, 5) or any larger instances. Both the exact algorithm Voting in [20] and the exact algorithm RISOTTO in [21] can solve the challenging instances up to (15, 5). In most of the cases, Voting runs faster than RISOTTO. The best uptodate exact algorithm is Pampa given in [10]. Pampa can solve the challenging instance (19, 7) within about 4.8 hours. The second best exact algorithm is PMSPrune [22] that can solve the challenging instance (19, 7) within about 10 hours.
In this paper we present an exact algorithm (named PMS5) that can solve the challenging instances (21, 8) and (23, 9). PMS5 takes about 10 hours on the challenging instance (21, 8) and about 54 hours on the challenging instance (23, 9). These run times are on a single 2.4GHz PC with 3GB of RAM. To the best of our knowledge, no other exact algorithm can solve these instances.
2 Methods
2.1 Notations and Definitions
In this section we introduce some notations and definitions that will help us describe our algorithm clearly.
Definition 2.1. A string x = x[1] ... x[ℓ] of length ℓ is called an ℓmer.
Definition 2.2. Given two strings x and y of equal length, we say the Hamming distance between x and y, denoted by d_{ H }(x, y). is the number of mismatches between them,
Definition 2.3. Given a string x = x[1] ... x[ℓ], we define the dneighborhood of x, B_{ d }(x), to be {y  d_{ H }(x, y) ≤ d}.
Note that ${B}_{d}(x)={\displaystyle {\sum}_{i=0}^{d}\left(\begin{array}{c}i\\ \mathrm{l}\end{array}\right)}{\left(\left\Sigma \right1\right)}^{i}$, where Σ is the alphabet of interest. Notice that B_{ d }(x) depends only on ℓ, d and Σ. For this reason, we define $\mathcal{N}(\ell ,d)={\displaystyle {\sum}_{i=0}^{d}\left(\begin{array}{c}i\\ \mathrm{l}\end{array}\right){\left(\left\Sigma \right1\right)}^{i}}$.
 1.
We use the notation x ∈_{ℓ} s if x is a substring of s (equivalently, s = αxβ for some strings α and β). We also say that x is an ℓmer of s.
 2.
We define ${\stackrel{\u0304}{d}}_{H}\left(x,s\right)={min}_{r{\in}_{\ell}s}{d}_{H}\left(x,r\right)$.
Definition 2.5. Given a string x = x[1] ... x[ℓ] and a set of strings $\mathcal{S}=\left\{{s}_{1},...,{s}_{n}\right\}$ with s_{ i } = m for i = 1, ..., n and ℓ < m, we define ${\stackrel{\u0304}{d}}_{H}\left(x,\mathcal{S}\right)=\underset{1\le i\le n}{max}{\stackrel{\u0304}{d}}_{H}\left(x,{s}_{i}\right)$.
It is easy to see that x is an (ℓ, d)motif of $\mathcal{S}$ if and only if ${\stackrel{\u0304}{d}}_{H}\left(x,\mathcal{S}\right)\le d$.
Definition 2.6. Given a set of strings $\mathcal{S}=\left\{{s}_{1},...,{s}_{n}\right\}$, we define ${M}_{\ell ,d}\left(\mathcal{S}\right)$ to be the set of (l, d) motifs of $\mathcal{S}$.
The goal of PMS is to compute ${M}_{\ell ,d}\left(\mathcal{S}\right)$, given ℓ, d and $\mathcal{S}$.
2.2 PMS5  A fast algorithm
The idea of our algorithm is based on the following observations about ${M}_{\ell ,d}\left(\mathcal{S}\right)$.
Observation 2.1. Let $\mathcal{S}$, ${\mathcal{S}}^{\prime}$ and ${\mathcal{S}}^{\u2033}$ be three sets of strings such that $\mathcal{S}={\mathcal{S}}^{\prime}\cup {\mathcal{S}}^{\u2033}$ and ${S}^{\prime}\cap {S}^{\u2033}=\varnothing $. It is easy to observe that ${M}_{\ell ,d}\left(\mathcal{S}\right)={M}_{\ell ,d}\left({\mathcal{S}}^{\prime}\right)\cap {M}_{\ell ,d}\left({\mathcal{S}}^{\u2033}\right)$.
Observation 2.2. For any string s, ${M}_{\ell ,d}\left(\left\{s\right\}\right)={\bigcup}_{x{\in}_{\ell}s}{B}_{d}\left(x\right)$.
From Observation 2.1 and Observation 2.2, we can obtain the following observation.
Observation 2.3. Let ${\mathcal{S}}^{*}=\mathcal{S}\backslash \left\{{s}_{1}\right\}=\left\{{s}_{2},...,{s}_{n}\right\}$. We have ${M}_{\ell ,d}\left(\mathcal{S}\right)={\bigcup}_{x{\in}_{\ell}{s}_{1}}\left[{B}_{d}\left(x\right)\cap {M}_{\ell ,d}\left({\mathcal{S}}^{*}\right)\right]$.
Observation 2.3 tells us that ${M}_{\ell ,d}\left(\mathcal{S}\right)$ can be computed from ${B}_{d}\left(x\right)\cap {M}_{\ell ,d}\left({\mathcal{S}}^{*}\right)$.
Without loss of generality, we can assume that the size of ${\mathcal{S}}^{*}$ is even, i.e., $\left{\mathcal{S}}^{*}\right=n1=2p$, for some integer p. Otherwise we can add a copy of s_{ n }into ${\mathcal{S}}^{*}$, in which case ${M}_{\ell ,d}\left({\mathcal{S}}^{*}\right)$ will remain the same. Next, we partition ${\mathcal{S}}^{*}$ into pairs of strings ${\mathcal{S}}_{1},...,{\mathcal{S}}_{p}$, where ${\mathcal{S}}_{k}=\left\{{s}_{2k},{s}_{2k+1}\right\}$ for k = 1 ... p. From Observations 2.1 and 2.2, we can make the following observations.
Based on the above observations, we note that the process of computing ${M}_{\ell ,d}\left(\mathcal{S}\right)$ can be reduced to computing B_{ d }(x, y, z) = B_{ d }(x) ∩ B_{ d }(y) ∩ B_{ d }(z) repeatedly. We will discuss how to compute B_{ d }(x, y, z) efficiently in Section 2.2.2. The pseudocode of our algorithm PMS5 is given below.
Algorithm PMS5
1: for each x ∈_{ℓ} s_{1} do
2: for k = 1 to $p=\u230a\frac{n1}{2}\u230b$ do
3: $Q\leftarrow \varnothing $.
4: for each y ∈_{ℓ} s_{2k}and z ∈_{ℓ} s_{2k+1}do
5: Compute B_{ d }(x, y, z) = B_{ d }(x) ∩ B_{ d }(y) ∩ B_{ d }(z).
6: Q ← Q ∪ B_{ d }(x, y, z).
7: end for
8: if k = 1 then
9: Q' ← Q.
10: else
11: Q' ← Q' ∩ Q.
12: end if
13: if Q' is small enough then
14: break the for loop.
15: end if
16: end for
17: for each r ∈ Q' do
18: if ${\stackrel{\u0304}{d}}_{H}\left(r,\mathcal{S}\right)\le d$ then
19: Output r as an (ℓ, d) motif.
20: end if
21: end for
22: end for
In the pseudo code, the process of computing ${B}_{d}\left(x\right)\cap {M}_{\ell ,d}\left({\mathcal{S}}_{k}\right)$ for each k is from line 3 to line 7. After line 7, Q is actually ${B}_{d}\left(x\right)\cap {M}_{\ell ,d}\left({\mathcal{S}}_{k}\right)$. Within the loop at line 2, Q' is ${B}_{d}\left(x\right)\cap {M}_{\ell ,d}\left({\mathcal{S}}_{1}\right)\cap \cdots \cap {M}_{\ell ,d}\left({\mathcal{S}}_{k}\right)$ for each k after line 12. At line 13, if Q' is less than a certain threshold, the algorithm simply exits the loop and will not try other values of k. In practice, we set the threshold to be between 5000 and 10000. From line 17 to line 21, the algorithm checks if each string r ∈ Q' is actually an (ℓ, d)motif or not. To check if ${\stackrel{\u0304}{d}}_{H}\left(r,\mathcal{S}\right)\le d$ for any r, we only have to use the remaining sequences (s_{2k+2}, s_{2k+3}, ..., s_{ n }).
2.2.1 Analysis
PMS5 is correct
From the observations, it is straightforward to see that PMS5 outputs ${M}_{\ell ,d}\left(\mathcal{S}\right)$. Therefore, PMS5 is correct.
The worstcase run time of PMS5
Theorem 2.1. The worstcase run time of PMS5 is $O\left(n{m}^{3}d\mathcal{N}\left(\ell ,d\right)\right)$. Recall that $\mathcal{N}(\ell ,d)={B}_{d}(x)={\displaystyle {\sum}_{i=0}^{d}\left(\begin{array}{c}i\\ \mathrm{l}\end{array}\right)}{\left(\left\Sigma \right1\right)}^{i}$.
Proof. It is easy to see that the run time of PMS5 is dominated by the computation time of B_{ d }(x, y, z) in line 5. In Section 2.2.2, we will show that B_{ d }(x, y, z) can be computed in O(ℓ + dB_{ d }(x, y, z)) time. In the extreme case in which $x=y=z,{B}_{d}(x,y,z))={B}_{d}(x)=\mathcal{N}(\ell ,d)$. Since $\mathcal{N}\left(\ell ,d\right)$ is much larger than ℓ, the computation time of B_{ d }(x, y, z) is $O\left(d\mathcal{N}\left(\ell ,d\right)\right)$. Also, it is straightforward to see that the number of times B_{ d }(x, y, z) is computed is at most $\frac{n}{2}{\left(m\ell +1\right)}^{3}$. Hence, the run time of PMS5 is $O\left(n{m}^{3}d\mathcal{N}\left(\ell ,d\right)\right)$.
The expected run time of PMS5
We can compute the expected run time of of PMS5 by computing the expected value of B_{ d }(x, y, z). Let x, y, and z be three random ℓmers. How many ℓmers are there that are at a distance of ≤ d from each of x, y, and z? Let u be a random ℓmer. Prob.$[{d}_{H}(x,u)\le d]={p}_{\ell ,d}={\displaystyle {\sum}_{i=0}^{d}\left(\begin{array}{c}i\\ \mathrm{l}\end{array}\right)}{(3/4)}^{i}{(1/4)}^{\ell i}$. This means that Prob.$\left[{d}_{H}\left(x,u\right)\le d\&{d}_{H}\left(y,u\right)\le d\&{d}_{H}\left(z,u\right)\le d\right]={p}_{l,d}^{3}$. Therefore, the expected number of u's such that u is at a distance of ≤ d from each of x, y, and z, E[B_{ d }(x, y, z)], is ${4}^{\ell}{p}_{\ell ,d}^{3}$.
As a result, the expected run time of PMS5 is $O\left(n{m}^{3}d{4}^{\ell}{p}_{\ell ,d}^{3}\right)$, where ${p}_{\ell ,d}={\displaystyle {\sum}_{i=0}^{d}\left(\begin{array}{c}i\\ \mathrm{l}\end{array}\right){(3/4)}^{i}{(1/4)}^{\ell i}}$.
A comparison between $\mathcal{N}\left(\ell ,d\right)$ and E[B_{ d }(x, y, z)] for different values of ℓ and d
ℓ  d  $\mathcal{N}\left(\ell ,d\right)$  E[B_{ d }(x, y, z)] 

9  2  352  6.35 × 10^{4} 
11  3  4,984  7.04 × 10^{3} 
13  4  66,379  6.49 × 10^{2} 
15  5  853,570  5.39 × 10^{1} 
17  6  1.07 × 10^{7}  4.20 
19  7  1.33 × 10^{8}  3.12 × 10 
21  8  1.63 × 10^{9}  2.26 × 10^{2} 
23  9  1.99 × 10^{10}  1.60 × 10^{3} 
2.2.2 Computing the intersection of the dneighborhoods
In this section, we consider the problem of computing the intersection of the dneighborhoods B_{ d }(x, y, z). Given three ℓmers x, y, z and integer number d, we would like to list all of the ℓmers in B_{ d }(x, y, z). In this section we offer an algorithm FULLPRUNE for this task that runs in O(ℓ + dB_{ d }(x, y, z)) time.
FULLPRUNE is the heart of algorithm PMS5. The idea of FULLPRUNE is as follows. We first represent B_{ d }(x) as a tree ${\mathcal{T}}_{d}\left(x\right)$ in which each node is an ℓmer in B_{ d }(x) and its root is the ℓmer x. The depth of ${\mathcal{T}}_{d}\left(x\right)$ is d. We will describe ${\mathcal{T}}_{d}\left(x\right)$ in detail later. We traverse ${\mathcal{T}}_{d}\left(x\right)$ in a depthfirst manner. At each node t during the traversal, we output t if t is in B_{ d }(y) ∩ B_{ d }(z). We also check if there is a descendent t' of t such that t' is in B_{ d }(y) ∩ B_{ d }(z). If there is no such descendent, we prune the subtree rooted at node t. We will show that checking the existence of such a descendent can be done quickly in O(1) time, later. Formally, ${\mathcal{T}}_{d}\left(x\right)$ is constructed from the following rules.
 1.
Each node in ${\mathcal{T}}_{d}\left(x\right)$ is a pair (t, p) where t = t[1] ... t[ℓ] is an ℓmer and p is an integer between 0 and ℓ such that t[p + 1] ... t[ℓ] = x[p + 1] ... x[ℓ]. We refer to a node (t, p) as ℓmer t if p is clear.
 2.
Let t = t[1] ... t[ℓ] and t' = t'[1] ... t'[ℓ]. A node (t, p) is the parent of a node (t', p') if and only if
 (a)
p' > p.
 (b)
t'[p'] ≠ t[p'] (From Rule 1, t[p'] = x[p']).
 (c)
t'[i] = t[i] for any i ≠ p'
 3.
The root of ${\mathcal{T}}_{d}\left(x\right)$ is (x, 0).
 4.
The depth of ${\mathcal{T}}_{d}\left(x\right)$ is d.
It is not hard to see that ${\mathcal{T}}_{d}\left(x\right)$ has the following properties.
 1.
If a node (t', p') is a child of a node (t, p), then d_{ H }(x, t')  d_{ H }(x, t) = d_{ H }(t, t') = 1. As a result, if a node (t, p) at level h, then d_{ H }(x, t) = h.
 2.
Consider two nodes (t, p) and (t', p') with t = t[1] ... t[ℓ] and t' = t'[1] ... t'[ℓ]. Then (t', p') is a descendent of (t, p) if and only if:
 (a)
p' > p.
 (b)
t'[1] ... t'[p] = t[1] ... t[p].
 (c)
d_{ H }(x, t') ≤ d.
 1.
d_{ H }(x, t') = d_{ H }(x_{1}, t_{1}) + d_{ H }(x_{2}, w) ≤ d.
 2.
d_{ H }(y, t') = d_{ H }(y_{1}, t_{1}) + d_{ H }(y_{2}, w) ≤ d.
 3.
d_{ H }(z, t') = d_{ H }(z_{1}, t_{1}) + d_{ H }(z_{2}, w) ≤ d.
We will show that the constraints can be expressed by an integer linear program of ten variables. Each location i in x_{2}, y_{2} and z_{2} is one of five types.

Type 1 (or Type aaa): x_{2}[i] = y_{2}[i] = z_{2}[i].

Type 2 (or Type aab): x_{2}[i] = y_{2}[i] ≠ z_{2}[i].

Type 3 (or Type aba): x_{2}[i] = z_{2}[i] ≠ y_{2}[i].

Type 4 (or Type baa): x_{2}[i] ≠ y_{2}[i] = z_{2}[i].

Type 5 (or Type abc): x_{2}[i] ≠ y_{2}[i], x_{2}[i] ≠ z_{2}[i], y_{2}[i] ≠ z_{2}[i].
Let n_{1} (resp. n_{2}, n_{3}, n_{3}, n_{4}, and n_{5}) be the number of locations of Type 1 (resp. Type 2, Type 3, Type 4, and Type 5). Given x_{2}, y_{2} and z_{2}, n_{ j }is determined for j = 1 ... 5. Notice that n_{1} + ··· + n_{5} = ℓ  p.
Consider any (ℓ  p)mer w = w[1] ... w[ℓ  p]. We define the following variables.

Let N_{1,a}be the number of locations i of Type 1 such that w[i] = x_{2}[i]. We should have N_{1,a}≤ n_{1}.

Let N_{2},_{ a }(resp. N_{2,b}) be the number of locations i of Type 2 such that w[i] = x_{2}[i] (resp. w[i] = z_{2}[i]). We should have N_{2,a}+ N_{2,b}≤ n_{2}.

Let N_{3,a}(resp. N_{3},_{ b }) be the number of locations i of Type 3 such that w[i] = x_{2}[i] (resp. w[i] = y_{2}[i]). We should have N_{3,a}+ N_{3,b}≤ n_{3}.

Let N_{4,a}(resp. N_{4},_{ b }) be the number of locations i of Type 4 such that w[i] = y_{2}[i] (resp. w[i] = x_{2}[i]). We should have N_{4,a}+ N_{4,b}≤ n_{4}.

Let N_{5,a}(resp. N_{5},_{ b }, N_{5},_{ c }) be the number of locations i of Type 5 such that w[i] = x_{2}[i] (resp. w[i] = y_{2}[i], w[i] = z_{2}[i]). We should have N_{5,a}+ N_{5,b}+ N_{5,c}≤ n_{4}.
Now it is straightforward to calculate d_{ H }(x_{2}, w) through the variables. The number of mismatches between x_{2} and w caused by the locations of Type 1 (resp. Type 2, Type 3, Type 4, and Type 5) is n_{1}  N_{1,a}, (resp. n_{2}  N_{2,a}, n_{3}  N_{3,a}, n_{4}  N_{4,b}, and n_{5}  N_{5,a}). Therefore, d_{ H }(x_{2}, w) = (n_{1}  N_{1,a}) + (n_{2}  N_{2,a}) + (n_{3}  N_{3,a}) + (n_{4}  N_{4,b}) + (n_{5}  N_{5,a}). Similarly, d_{ H }(y_{2}, w) = (n_{1}  N_{1,a}) + (n_{2}  N_{2,a}) + (n_{3}  N_{3,b}) + (n_{4}  N_{4,a}) + (n_{5}  N_{5,b}), and d_{ H }(z_{2}, w) = (n_{1}  N_{1,a}) + (n_{2}  N_{2,b}) + (n_{3}  N_{3,a}) + (n_{4}  N_{4,a}) + (n_{5}  N_{5,c}). So the following integer linear program (ILP) expresses the constraints.
 1.
(n_{1}  N_{1,a}) + (n_{2}  N_{2,a}) + (n_{3}  N_{3,a}) + (n_{4}  N_{4,b}) + (n_{5}  N_{5,a}) ≤ d  d_{ H }(x_{1}, t_{1}).
 2.
(n_{1}  N_{1},_{ a }) + (n_{2}  N_{2,a}) + (n_{3}  N_{3,b}) + (n_{4}  N_{4,a}) + (n_{5}  N_{5,b}) ≤ d  d_{ H }(y_{1}, t_{1}).
 3.
(n_{1}  N_{1,a}) + (n_{2}  N_{2,b}) + (n_{3}  N_{3},_{ a }) + (n_{4}  N_{4,a}) + (n_{5}  N_{5,c}) ≤ d  d_{ H }(z_{1}, t_{1}).
 4.
N_{1, a}≤ n_{1}.
 5.
N_{2, a}+ N_{2, b}≤ n_{2}.
 6.
N_{3, a}+ N_{3, b}≤ n_{3}.
 7.
N_{4, a}+ N_{4, b}≤ n_{4}.
 8.
N_{5, a}+ N_{5, b}+ N_{5},_{ c }≤ n_{5}.
 9.
All of the variables are nonnegative integers.
 1.
N_{1, a}+ N_{2, a}+ N_{3},_{ a }+ N_{4},_{ b }+ N_{5,a}≥ ℓ  p  d + d_{ H }(x_{1}, t_{1}).
 2.
N_{1, a}+ N_{2, a}+ N_{3},_{ b }+ N_{4},_{ a }+ N_{5},_{ b }≥ ℓ  p  d + d_{ H }(y_{1}, t_{1}).
 3.
N_{1, a}+ N_{2, b}+ N_{3,a}+ N_{4,a}+ N_{5,c}≥ℓ  p  d + d_{ H }(z_{1}, t_{1}).
Because the integer linear program has only ten variables, checking whether it has a solution can be done in O(1) time. Notice that the integer linear program only depends on eight parameters n_{1}, ... n_{5}, d  d_{ H }(x_{1}, t_{1}), d  d_{ H }(y_{1}, t_{1}), and d  d_{ H }(z_{1}, t_{1}). The first five parameters are in the range [0, ..., ℓ] and the other parameters are in the range [0, ... d]. Therefore, we will store the results of all possible integer linear programs in a 8dimensional table of size (ℓ + 1)^{5}(d + 1)^{3} to speedup the checking time for the integer linear programs during the traversal on the tree in FullPrune. Notice that we only need to compute the table once before FULLPRUNE is executed, and reuse it as many times as needed. The pseudocode of FULLPRUNE is given below.
 1.
Compute d_{ H }(x, y) and d_{ H }(x, z).
 2.
Compute n_{1}, n_{2}, n_{3}, n_{4} and n_{5} for each p = 0... (ℓ  1).
 3.
Traverse the tree ${\mathcal{T}}_{d}\left(x\right)$ in a depthfirst manner. At each node (t, p), do the following steps.
 (a)
Incrementally compute d_{ H }(x, t), d_{ H }(y, t), and d_{ H }(z, t) from its parent.
 (b)
Incrementally compute d_{ H }(x_{1}, t_{1}), d_{ H }(y_{1}, t_{1}), and d_{ H }(z_{1}, t_{1}) from its parent. (Notice that t_{1} = t[1] ... t[p], x_{1} = x[1] ... x[p], y_{1} = y[1] ... y[p] and x_{1} = z[1] ... z[p]).
 (c)
If d_{ H }(x, t) ≤ d, d_{ H }(y, t) ≤ d and d_{ H }(z, t) ≤ d, then output t.
 (d)
Solve the integer linear program (ILP) with parameters n_{1}, n_{2}, n_{3}, n_{4}, n_{5}, ℓ  p  d + d_{ H }(x_{1}, t_{1}), ℓ  p  d + d_{ H }(y_{1}, t_{1}), and ℓ  p  d + d_{ H }(z_{1}, t_{1}).
 (e)
If d_{ H }(x, t) ≥ d and/or the ILP does not have a solution, then prune the subtree rooted at node (t, p). Otherwise, explore its children.
Theorem 2.2. Given three ℓmers x, y and z, FULLPRUNE computes B_{ d }(x, y, z) in O(ℓ + dB_{ d }(x, y, z)) time.
Proof. From the discussion above, FULLPRUNE outputs all of the ℓmers in B_{ d }(x, y, z). Now let us analyze its run time. In the pseudocode of FullPrune, step 1 and step 2 take O(ℓ) time. We will show that step 3 takes at most O(dB_{ d }(x, y, z)) time, that will complete our proof. Since in ${\mathcal{T}}_{d}\left(x\right)$ a node and its parent differ at exactly one location, step 3a and step 3b take at most O(1) time. It is easy to see that the other steps inside step 3 (from step 3c to step 3e) also take O(1) time. Therefore, FULLPRUNE spends at most O(1) time at each node it visits. As a result, the run time of step 3 is proportional to the number of the visited nodes. We will argue that the number of visited nodes is no more than dB_{ d }(x, y, z). Consider the tree $\mathcal{T}$ consisting of all the nodes visited by FullPrune. Obviously, ${\mathcal{T}}_{d}\left(x\right)$ contains $\mathcal{T}$. Because of the property of the integer linear program, every leaf in $\mathcal{T}$ is an element in B_{ d }(x, y, z). Therefore, the number of leaves in $\mathcal{T}$ is at most B_{ d }(x, y, z). On the other hand, in any tree the number of nodes is no more than its depth times the number of its leaves. Since ${\mathcal{T}}_{d}\left(x\right)$ contains $\mathcal{T}$, the depth of $\mathcal{T}$ is less than or equal to the depth of ${\mathcal{T}}_{d}\left(x\right)$, which is equal to d. Hence, the number of nodes in $\mathcal{T}$, which is equal to the number of nodes visited by FullPrune, is at most dB_{ d }(x, y, z).
We conclude this section with a remark that our algorithm FULLPRUNE can be generalized as follows. Right now we use the computation of the common dneighborhood of three ℓmers as the basic step. This can be generalized so that the basic step is that of computing the common dneighborhood of k ℓmers (for any value of k ≤ n).
2.3 Extended PMS5 for Solving the qPMS Problem
In this section, we consider a generalized version of the PMS problem called the qPMS Problem (see e.g., [22]). In the qPMS problem, we relax the constraints on the motifs. An ℓmer x is a motif if there are at least q sequences s_{ i }in $\mathcal{S}$ such that d_{ H }(x, s_{ i }) ≤ d. Apparently, the qPMS problem becomes the PMS problem if q = n. In practice, the qPMS problem is a more realistic model of motifs since these motifs usually appear in some of the given sequences, instead of appearing in all of them.
We can extend the algorithm PMS5 to solve the qPMS problem by exploiting the heart of PMS5, i.e., the algorithm FULLPRUNE that computes B_{ d }(x,y, z). One simple and straightforward way to extend PMS5 for the qPMS problem is as follows. We consider every tuple of sequences (s_{ i }, s_{ j }, s_{ k }), 1 ≤ i < j < k ≤ n. For each tuple (s_{ i }, s_{ j }, s_{ k }), we compute B_{ d }(x, y, z) where x, y, and z are in s_{ i }, s_{ j }and s_{ k }, respectively. For each ℓmer t in B_{ d }(x, y, z), we check whether there are at least q3 sequences s_{ p }in $\mathcal{S}\backslash \left\{{s}_{i},{s}_{j},{s}_{k}\right\}$ such that d_{ H }(t, s_{ p }) ≤ d. If t satisfies this constraint, we output t as a motif. The psuedocode is provided below.
Extended Algorithm PMS5 for qPMS
1: for each tuple of sequences (s_{ i }, s_{ j }, s_{ k }), where 1 ≤ i < j < k ≤ n do
2: for each tuple (x, y, z) of ℓmers where x ∈_{ℓ} s_{ i },y ∈_{ℓ} s_{ j }, and z ∈ _{ℓ} s_{ k }do
3: Compute B_{ d }(x, y, z) using FULLPRUNE.
4: for each t ∈ B_{ d }(x, y, z) do
5: if there are at least q3 sequences ${s}_{p}\in \mathcal{S}\backslash \left\{{s}_{i},{s}_{j},{s}_{k}\right\}$ such that d_{ H }(t, s_{ p }) ≤ d then
6: output t.
7: end if
8: end for
9: end for
10: end for
The two following theorems are immediate:
Theorem 2.3. The worst run time of the above algorithm is $O\left({n}^{4}{m}^{3}d\mathcal{N}\left(\ell ,d\right)\right)$.
Theorem 2.4. The expected run time of the above algorithm is $O\left({n}^{4}{m}^{3}d{4}^{\ell}{p}_{\ell ,d}^{3}\right)$, where ${p}_{\ell ,d}={\displaystyle {\sum}_{i=0}^{d}\left(\begin{array}{c}l\\ \mathrm{i}\end{array}\right){(3/4)}^{i}{(1/4)}^{\ell i}}$
2.4 Challenging Instances for qPMS
The challenging instances for qPMS have to be defined appropriately. For every value of ℓ, we can define a corresponding challenging instance with a relevant value for d. We define the challenging instance corresponding to a given value of ℓ to be the pair (ℓ, d) if d is the smallest value for which the expected number of (ℓ, d)motifs that occur by random chance is at least 1. In fact the same definition is used for the PMS problem as well. However, the identification of such instances is slightly different. We could identify the challenging instances for qPMS as follows. Let S_{1}, S_{2}, ..., S_{ n }be the given input strings. Consider a random ℓmer w. Let S be any one of the input strings and x be an ℓmer in S.
Probability that the Hamming distance between w and x is ≤ d is $P={\displaystyle {\sum}_{i=0}^{d}\left(\begin{array}{c}l\\ \mathrm{i}\end{array}\right){\left({\scriptscriptstyle \frac{3}{4}}\right)}^{i}{\left({\scriptscriptstyle \frac{1}{4}}\right)}^{\ell i}}$. Probability that the Hamming distance between w and x is > d is (1  P). Probability that the Hamming distance between w and each ℓmer of S is > d is Q' = (1  P)^{ℓm+1}. Here we assume that the ℓmers of S are independent, which is clearly incorrect. A similar assumption has been made in prior analyses (see e.g., [3]) and in practice conclusions made using such analyses seem to hold nearly. Probability that S has at least one ℓmer x such that the Hamming distance between w and x is ≤ d is Q = 1  Q'. If the Hamming distance between w and x is ≤ d, call x as an instance of w.
The expected number of random motifs for qPMS corresponding to various values of ℓ and d with n = 20, m = 600 and q = 10
1  d  Expected Number of Random Motifs 

9  2  1.599 
9  1  0.159 
11  2  1.424 
11  1  8.643 × 10^{12} 
13  3  22.090 
13  2  1.530 × 10^{9} 
15  4  154 
15  3  7.150 × 10^{8} 
17  5  640 
17  4  5.277 × 10^{6} 
19  6  1883 
19  5  8.504 × 10^{6} 
3 Results and Discussion
3.1 Performance of PMS5 on the challenging instances
In this section, we show the performance of PMS5 on the challenging instances as described in Section 1. We also compare the performance of PMS5 with that of other wellknown exact algorithms such as Pampa [10], PMSPrune [22], Voting [20], and RISSOTO [21]. Algorithms for planted motif search are typically tested on random input datasets. Any such dataset will consist of 20 random strings each of length 600 (n = 20, m = 600). A random motif of length ℓ is planted at a random position in each of the strings, mutating it in exactly d places. We test the algorithms for varying ℓ and d values. In particular, we have employed the following challenging instances: (13, 4), (15, 5), (17, 6), (19, 7), (21, 8), and (23, 9).
To have a fair comparison, we have run all of the algorithms on the same machine. The configuration of the machine is Windows XP Operating System, Dual Core Pentium 2.4GHz CPU with 3GB RAM. PMS5 is written in C language. Pampa, PMSPrune and RISSOTO were also written in C language. Only Voting was written in C++. All of the algorithms were compiled using Microsoft Visual C++ 6.0 Compiler.
Time comparison on challenging instances
Algorithm  (13,4)  (15,5)  (17,6)  (19,7)  (21,8)  (23,9) 

PMS5  117s  4.8 m  21.7 m  1.7h  9.7h  54h 
Pampa  35s  6 m  40 m  4.8h     
PMSPrune  45s  10.2 m  78.7 m  15.2h     
Voting  104s  21.6 m         
RISOTTO  772s  106.4 m         
While comparing PMS5 and PMSPrune, we notice an interesting fact that as the challenging instance increases in size, the ratio between their run times increase exponentially. In particular, the ratio is roughly 2,4, and 8 on the challenging instances (15,5), (17,6), and (19,7), respectively. This fact perhaps explains why PMS5 can solve the challenging instances (21, 8) and (23, 9) but PMSPrune cannot. If this observation is true in general, PMSPrune will probably take about 16 × 9.7 = 155.2 hours on the instance (21, 8), and 32 × 54 = 1728 hours on the instance (23, 9).
Notice that the run time of PMS5 does not include the time for building the ILP table. It takes 1.5 hours and 500MB to build and store the ILP table for ℓ = 23 and d = 9.
3.2 PMS5 on real data: predicting transcript factorbinding sites
In this section, we discuss how to use algorithm PMS5 to find transcript factorbinding sites in the real data provided in [23]. The real data is broadly used to test many existing algorithms [23], [11], [22], [3]. Each dataset in the real data is comprised of DNA sequences with verified transcript factorbinding sites. The datastes are from many species including mouse, human and yeast.
We have used the algorithm PMS5 to find transcript factorbinding sites as follows. For any given dataset, we have run PMS5 with ℓ = 21, d = 8, and obtained a set of motifs. Some of these motifs could be spurious hits. Hence, we need a scoring scheme to rank the motifs. We have used the simple scoring function ∑_{1≤i≤n}d_{ H }(M, s_{ i }), where d_{ H }(M, s_{ i }) is the hamming distance between motif M and sequence s_{ i }. We take the motif with the lowest score and then predict transcription factorbinding sites based on it. Notice that we have only used one value for ℓ (namely, 21) because smaller values of ℓ have been tested in [22].
PMS5 on real datasets: predicting transcript factorbinding sites
Dataset  Best motif found by PMS5  Matched binding sites at: 

mus05r  AGAGGAAAAAAAAAAGGAGAG  seq 1: GGAAAAACAAAGGTAATG 
mus07r  GCTGCCCACCCTCTGCAACCC  seq 4: CCCAACACCTGCTGCCTGAGCC 
mus11r  AGGGCGGGGGGCGGAGCGGGG  seq 2: GCCGCCGGGGTGGGGCTGAG 
seq 3: GGGGGGGGGGGCGGGGC  
seq 4: GTGGGGGCGGGGCCTT  
seq 9: GAACAGGAAGTGAGGCGG  
hm03r  AAAAGAAAAAAAAATAAACAA  seq 1: CGGGTGTTATTCAAGCAAAAAAAATAAATAAATACCTATGCAATAC 
seq 2: GGATGTTACACAAGCAAACAAAATAAATATCTGTGCAATAT  
seq 3: TGGGTGTTATATGAGCAAACAAAATAAATACCTGTGCAACAT  
hm08r  CAGCGTGCAGTCCCCTTCATC  seq 10: TATGGTCATGACGTCTGACAGAGC 
hm19r  CCCCCTTCCACCACCCACAGA  seq 2: CACTTTTAGCTCCTCCCCCCA 
hm20r  CCTCCTTCCTCCCCCTCCCCC  seq 10: TCCTCCCCACCTTCCCCACCCTCCCCACCCTCCCCATAAGCGCCCCTCCCG 
seq 11: GCAAACTCCGCCTCCCCCAA  
seq 14: GTCCCTCCTCCTCCCGCCC  
hm22r  GGACACGGCAGAGCCTGGGGA  seq 4: GAGGCAGACCACGTGAGAGCCTGGCCAGGCCTTCC 
hm24r  CGCCTGCGCCCCGCCCCGCCC  seq 2: CCCCGCCCCGCGCTCCCC 
hm26r  CCCCCCGCCTCCCGCTCCCAG  seq 3: CCCCGCCTCAGGCTCCCGGGG 
seq 7: CTCAGCCTGCCCCTCCCAGGGATTAAG  
seq 8: GCGCCGAGGCGTCCCCGAGGCGC 
3.3 Performance of Extended PMS5 on the qPMS challenging instances
Run time of Extended PMS5 on qPMS challenging instances
qPMS challenging instance  Run time 

(9,1)  78s 
(11, 2)  7.48 m 
(13, 3)  18.43 m 
(15, 4)  1.39h 
(17, 5)  15.93h 
(19, 6)   
4 Conclusions
In this paper we have presented an efficient exact algorithm for the (ℓ, d)motif search problem. This algorithm is more efficient than any known exact algorithm for PMS. In particular, using this algorithm we can solve the challenging instances (21, 8) and (23, 9). No prior exact algorithms could solve these instances. Therefore, we hope that PMS5 will help biologists discover longer motifs in future. Our algorithm is based on some novel ideas that will be of independent interest to solve PMS and other variations of the motif search problem. One of the basic ideas we employ is that of computing the common dneighborhood of three ℓmers. This is done using an integer programming formulation. An open problem will be to exploit this idea to further improve the performance of our algorithm. One possible direction is to use a basic step where the dneighborhood of k ℓmers is computed (for some relevant value of k). We have extended our algorithm to solve the qPMS problem as well. Challenging instances for the qPMS problem have been defined and computed. Our extended algorithm can solve the following qPMS challenging instances: (9,1), (11, 2), (13, 3), (15, 4), and (17, 5). In comparison, the exact algorithms MITRA, RISOTTO, and Voting also can only solve challenging instances up to d = 5 (but for the version where the motifs occur in all the input strings).
Declarations
Acknowledgements
This work has been supported in part by the following grants: NSF 0829916 and NIH 1R01LM01010101A1.
Authors’ Affiliations
References
 Rajasekaran S: Computational techniques for motif search. Frontiers in Bioscience 2009, 14: 5052–5065. 10.2741/3586View Article
 Bailey T, Elkan C: Fitting a mixture model by expectation maximization to discover motifs in biopolymers. Proceedings of Second International Conference on Intelligent Systems for Molecular Biology 1994, 28–36.
 Buhler J, Tompa M: Finding motifs using random projections. Proceedings of Fifth Annual International Conference on Computational Molecular Biology (RECOMB) 2001.
 Lawrence CE, Altschul SF, Boguski MS, Liu JS, Neuwald AF, Wootton JC: Detecting subtle sequence signals: a Gibbs sampling strategy for multiple alignment. Science 1993, 262: 208–214. 10.1126/science.8211139View ArticlePubMed
 Pevzner P, Sze SH: Combinatorial approaches to finding subtle signals in DNA sequences. Proceedings of Eighth International Conference on Intelligent Systems for Molecular Biology 2000, 269–278.
 Rocke E, Tompa M: An algorithm for finding novel gapped motifs in DNA sequences. Proceedings of Second International Conference on Computational Molecular Biology (RECOMB) 1998, 228–233.
 Keich U, Pevzner P: Finding motifs in the twilight zone. Bioinformatics 2002, 18: 1374–1381. 10.1093/bioinformatics/18.10.1374View ArticlePubMed
 Price A, Ramabhadran S, Pevzner PA: Finding subtle motifs by branching from sample strings. Bioinformatics 2003, 1: 1–7.
 Hertz G, Stormo G: Identifying DNA and protein patterns with statistically significant alignments of multiple sequences. Bioinformatics 1999, 15: 563–577. 10.1093/bioinformatics/15.7.563View ArticlePubMed
 Davila J, Balla S, Rajasekaran S: Pampa: An Improved Branch and Bound Algorithm for Planted (l, d) Motif Search. Tech. rep 2007.
 Rajasekaran S, Balla S, Huang CH: Exact algorithms for planted motif challenge problems. Journal of Computational Biology 2005, 12(8):1117–1128. 10.1089/cmb.2005.12.1117View ArticlePubMed
 Blanchette M, Schwikowski B, Tompa M: An exact algorithm to identify motifs in orthologous sequences from multiple species. Proceedings of Eighth International Conference on Intelligent Systems for Molecular Biology 2000, 37–45.
 Eskin E, Pevzner P: Finding composite regulatory patterns in DNA sequences. Bioinformatics 2002, S1: 354–363.View Article
 Brazma A, Jonassen I, Vilo J, Ukkonen E: Predicting gene regulatory elements in silico on a genomic scale. Genome Research 1998, 15: 1202–1215.
 Galas DJ, Eggert M, Waterman MS: Rigorous patternrecognition methods for DNA sequences: Analysis of promoter sequences from Escherichia coli. Journal of Molecular Biology 1985, 186: 117–128. 10.1016/00222836(85)902621View ArticlePubMed
 Sinha S, Tompa M: A statistical method for finding transcription factor binding sites. Proceedings of Eighth International Conference on Intelligent Systems for Molecular Biology 2000, 344–354.
 Staden R: Methods for discovering novel motifs in nucleic acid sequences. Computer Applications in the Biosciences 1989, 5(4):293–298.PubMed
 Tompa M: An exact method for finding short motifs in sequences, with application to the ribosome binding site problem. Proc. Seventh International Conference on Intelligent Systems for Molecular Biology 1999, 262–271.
 van Helden J, André B, ColladoVides J: Extracting regulatory sites from the upstream region of yeast genes by computational analysis of oligonucleotide frequencies. Journal of Molecular Biology 1998, 281(5):827–842. 10.1006/jmbi.1998.1947View ArticlePubMed
 Chin F, Leung H: Algorithms for Discovering Long Motifs. Proceedings of the Third AsiaPacific Bioinformatics Conference (APBC2005), Singapore 2005, 261–271.View Article
 Pisanti N, Carvalho A, Marsan L, Sagot MF: RISOTTO: Fast extraction of motifs with mismatches. Proceedings of the 7th Latin American Theoretical Informatics Symposium 2006, 757–768.
 Davila J, Balla S, Rajasekaran S: Fast and practical algorithms for planted ( l,d ) motif search. IEEE/ACM Transactions on Computational Biology and Bioinformatics 2007, 544–552.
 Blanchette M: Algorithms for Phylogenetic Footprinting. Proceedings of Fifth International Conference Computational Biology (RECOMB 2001) 2001.
 Tompa M, Li N, Bailey T, Church G, Moor BD, Eskin E, Favorov A, Frith M, Fu Y, Kent W, Makeev V, Mironov A, Noble W, Pavesi G, Pesole G, Regnier M, Simonis N, Sinha S, Thijs G, van Helden J, Vandenbogaert M, Weng Z, Workman C, Ye C, Zhu Z: Assessing Computational Tools for the Discovery of Transcription Factor Binding Sites. Nature Biotechnology 2005, 23: 137–144. 10.1038/nbt1053View ArticlePubMed
 Sharma D, Rajasekaran S, Dinh H: An Experimental Comparison of PMSPrune and Other Algorithms for Motif Search. CoRR abs/1108.5217 2011.
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.