Home | About | Journals | Submit | Contact Us | Français |

**|**BMC Bioinformatics**|**v.7; 2006**|**PMC1534070

Formats

Article sections

Authors

Related links

BMC Bioinformatics. 2006; 7: 279.

Published online 2006 June 2. doi: 10.1186/1471-2105-7-279

PMCID: PMC1534070

Weichun Huang: huang6/at/niehs.nih.gov; David M Umbach: umbach/at/niehs.nih.gov; Uwe Ohler: uwe.ohler/at/duke.edu; Leping Li: li3/at/niehs.nih.gov

Received 2006 May 8; Accepted 2006 June 2.

Copyright © 2006 Huang et al; licensee 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.

This article has been cited by other articles in PMC.

Identifying functional elements, such as transcriptional factor binding sites, is a fundamental step in reconstructing gene regulatory networks and remains a challenging issue, largely due to limited availability of training samples.

We introduce a novel and flexible model, the Optimized Mixture Markov model (OMiMa), and related methods to allow adjustment of model complexity for different motifs. In comparison with other leading methods, OMiMa can incorporate more than the NNSplice's pairwise dependencies; OMiMa avoids model over-fitting better than the Permuted Variable Length Markov Model (PVLMM); and OMiMa requires smaller training samples than the Maximum Entropy Model (MEM). Testing on both simulated and actual data (regulatory *cis*-elements and splice sites), we found OMiMa's performance superior to the other leading methods in terms of prediction accuracy, required size of training data or computational time. Our OMiMa system, to our knowledge, is the only motif finding tool that incorporates automatic selection of the best model. OMiMa is freely available at [1].

Our optimized mixture of Markov models represents an alternative to the existing methods for modeling dependent structures within a biological motif. Our model is conceptually simple and effective, and can improve prediction accuracy and/or computational speed over other leading methods.

Biological sequences, including DNA, RNA and proteins, contain functionally important motifs, such as transcription factor binding sites (TFBS), RNA splice sites, and protein domains. With the increasing-availability of genome sequences, identification of such functional motifs not only plays important roles in gene finding and function prediction but also is a fundamental step in reconstructing gene regulatory networks and in revealing gene evolutionary mechanisms [2-6].

A commonly used model for motif identification is the Weight Matrix Model (WMM) proposed by Staden [7], also called the Position Weight Matrix (PWM) or Mononucleotide Weight Matrix (MWM). A PWM is usually generated from a set of aligned instances of known motif sequences, using the observed position-specific base frequencies and/or prior information. Stormo and Fields [8] showed that the PWM score of a motif is proportional to the total binding energy contributed by individual bases. PWM has been used by many motif identification programs, *e.g*., Matlnspector [9] and Match [10], and performs reasonably well for motif identification. While a PWM can capture both nucleotide preferences at each position and different levels of position specificity, it does not account for functional dependencies between positions. Recent studies [11-15] indicate that there are often important interactions between positions, adjacent as well as non-adjacent, within a motif. The inability of the PWM to capture such dependencies is a limitation as the PWM model often produces a large number of false positives in a genome-wide scan [16].

Many models have been developed to incorporate position dependencies. Motif models, such as the Dinucleotide Weight Matrix Model (DWMM) [17] and the Weight Array Model (WAM) [18], can incorporate dependencies between adjacent positions. To incorporate further dependencies of non-adjacent positions, Ponomarenko *et al. *[19] extended DWMM by introducing the Oligonucleotide Weight Matrix model, which includes a comprehensive set of oliogonucleotide matrices classified into 5 biological function categories. A WAM could also be extended to a high order WAM in principle, *e.g*., windowed *2*^{nd }order WAM [2]. However, the exponentially increased number of parameters of these models makes them impractical due to insufficient training data. To address the weaknesses of WAM in incorporating long-range interactions, Burge and Karlin [2] proposed the Maximal Dependence Decomposition (MDD) model, which has a binary tree structure formed by a set of conditional WAMs. While the MDD model can capture non-adjacent dependencies through the conditional WAM models, it still requires a rather large number of training sequences, which are partitioned into smaller subsets to train all conditional WAMs. To alleviate the requirement of a large training set, Cai *et al. *[20] developed a Bayesian tree to model dependencies within RNA splice sites; Ellrott *et al. *[21] suggested a position order optimized Markov chain model, which reorders motif positions to bring distant but dependent positions into near neighbors. More recently, several other models have been developed, including Bayesian networks for modeling protein-DNA binding sites [22], Maximum Entropy Model (MEM) for splice site identification [23], Permuted Variable Length Markov Model (PVLMM) for finding transcription factor binding sites and splice sites [24]. For a biological motif with position dependencies, these models can show improvement in prediction accuracy over the models that assume independence. Incorporating position dependencies can also improve the accuracy of *de novo *motif discovery [25].

In this paper, we present a new and flexible motif model, the OMiMa, to incorporate position dependencies within a motif. OMiMa can not only adjust model complexity according to motif dependency structures but also minimize model complexity without compromising prediction accuracy. As an integrated part of OMiMa, we also introduce the Directed Neighbor-Joining (DNJ) method to optimally rearrange positions to minimize Markov order. We then describe and discuss the methods for selecting the best model. We implement our model into the OMiMa system that is freely available to the public.

Let *X*_{i }be the discrete random variable associated with position *i *in a biological motif **X **of length *w*. For DNA sequences, *X*_{i }takes values from set B = {*A, C, G, T*}; and for protein sequences, *X*_{i }takes values from 20 different amino acids. *X*_{i }follows a multinomial distribution. Let = *X*_{i-k}...*X*_{i-1 }and = *x*_{i-k}...*x*_{i-1}, where *k *= 0,..., *w *- 1; upper case **X **(*X*_{i}) is a random variable and lower case **x** (*x*_{i}) is a particular value. The denotes an empty sequence and Pr( = ) = 1. Additionally, let *X*_{-j }= *X*_{w-j}, *x*_{-j }= *x*_{w-j}, where *j *= 0...*w *- 1. If one uses the *k*^{th }order Markov model (*M*_{k}), the probability of observing a motif sequence x is just the product of conditional/transition probabilities. Let be a *k*^{th }order Markov model of a linear chain, and be a *k*^{th }order Markov model of a circular chain. The probability of a motif sequence is given by equation (1) for a linear chain and equation (2) for a circular chain, respectively.

Compared to a linear Markov chain, a circular Markov chain incorporates additional dependencies that may contain subtle signals that allow the model to distinguish true motifs from false ones, especially when false motifs are similar to true motifs.

Suppose a motif **X **can be divided into *m *independent sub-motifs, that is **X **= **Y**_{1},... **Y**_{m }and each sub-motif is modeled as an independent Markov chain, that is , then the probability of the sequence (x) given the Markov models is:

These independent Markov models, each of which is position-optimized for its corresponding sub-motif, form an Optimized Mixture of Markov models (OMiMa). An example of OMiMa is illustrated in Figure Figure1.1. However, for short motifs such as transcription factor binding sites, we can use a simple mixture of Markov models consisting of only one 0^{th }order and one *k*^{th }order chain (Figure (Figure2).2). For convenience, we refer such a mixture model as *'a 0-k mixture model'*. Since the *k*^{th }order Markov chain of *'a 0-k mixture model' *can be either linear or circular, we also use terms *'a 0-k mixture linear model' *and *'a 0-k mixture circular model' *to distinguish them. In the following, we describe methods for the general mixture Markov model, while we use the simple 0-k mixture model for our testing.

Conceivably, the different parts of a motif could have distinct roles in the interaction with their partners. Motif positions involved in the same role can be highly dependent, whereas those involved unrelated roles are likely independent. A mixture of Markov models seems an ideal fit by modeling different signals with different sub-models. A 0^{th }order Markov chain can effectively model strong signals such as those embedded in highly conserved positions where the probability of a certain base occurring is almost one. In addition, positions where base composition contributes little or nothing to motif function need no more complex model than a 0^{th }order Markov model. On the other hand, a higher order Markov model is necessary for detecting subtle dependency signals that can be essential for distinguishing true motifs from false ones.

To apply the mixture of Markov models to a motif, the first step is to dissect the motif into several independent sub-motifs, each of which is modeled as a Markov chain. For a given set of sequences of a motif, we employ chi-square tests to find significant pairwise dependencies between positions within the motif (see also [21]). Based on pairwise dependencies, motif positions are grouped into independent sets, each forming a Markov chain. The outline of our procedure for grouping motif positions is described in the following steps.

1. Calculate base frequencies for each position, and find highly conserved positions where the observed frequency of a certain base (almost) equals 1. These conserved positions then are put into set *H *as defined below.

where *f*(*i, x*) is the frequency of base x at position i, and B is the set of bases.

2. Place remaining positions in the set *M*, and calculate pairwise chi-square values for every pair of positions in *M*.

where B_{i }and B_{j }are the sets of bases observed in positions *i *and *j*, respectively; *O*(*x*_{i}, *x*_{j}) and *E*(*x*_{i}, *x*_{j}) are the observed and expected counts of pair (*x*_{i}, *x*_{j}), respectively. *E*(*x*_{i}, *x*_{j}) is the product of observed base frequencies *x*_{i }and *x*_{j}. The degrees of freedom of this test is (|B_{i}| - 1) × (|B_{j}| - 1), where |B_{i}| and |B_{j}| are the number of different bases in sets B_{i }and B_{j}, respectively.

3. Based on the above *χ*^{2 }tests, find all positions that show little dependence on any other positions in *M*, and move them to the set *I*, as defined by

Here *p*_{i,j }is the p-value corresponding to , and *α *is the significance level, *e.g*., 0.05.

4. The remaining positions in *M *are further grouped into subsets by iterating the following rules:

(a) Set s = 1.

(b) Calculate *θ*_{i }= ∑_{jM,j≠i }*δ*(*p*_{i,j }<*α*) for each position *i *in *M*, where *δ *is a 0/1 indicator function. Find the largest *θ*_{i}, and move position *i *and all positions *j *that *p*_{i,j }<*α *from *M *into a new set *C*_{s}.

(c) For each remaining position, check if it significantly depends on any position in *C*_{s}. If it does, then move it from *M *into *C*_{s}.

(d) If *M *is not empty, update s = s+1 and go back to step (b).

Step 4 above essentially groups positions into independent subsets, each potentially forming a functional unit. For the special 0-k mixture model, we simply set *M *= *C*_{1 }at this step.

The next step is to arrange the positions in each subset into a Markov chain. Since the positions in sets *H *and *I *are independent of each other, they can be arranged in their natural order to form a 0^{th }order Markov chain. The positions in *H *can also be treated differently from those in set *I *in motif identification by requiring a perfect match for a true site. Sets *C*_{s }are different. The position arrangement for each set *C*_{s }needs to be optimized so that the Markov model can account for most dependencies while minimizing the Markov order. For a given set *C*_{s}, we use the median (*K*_{s}) of *θ *(*θ *= {*θ*_{j}, *j * *C*_{s}}) as the maximum order of its potential Markov model. We then optimize position arrangement for the *k*^{th }order Markov chain (*k *= 0,..., *K*_{s}) by the Directed Neighbor-Joining (DNJ) method described below.

The neighbor-joining (NJ) method proposed by Saitou and Nei [26] is a well-known distance method for phylogenetic tree reconstruction. The principle of the NJ method is to find pairs of operational taxonomic units that minimize the total branch length at each stage of clustering. Our DNJ method is based on the exactly same principle. The only major difference is that DNJ needs to consider the direction in joining two nearest neighbors to form a new node while NJ does not. Instead of producing a phylogenetic tree as the NJ method does, DNJ method creates a chain structure, which arranges closely dependent positions as the nearest neighbors. The DNJ method for constructing a *k*^{th }order Markov chain from a given subset (*C*_{s}) of motif positions is described in the following steps (see Figure Figure33 for an example).

1. For a given set *C*_{s}, put each position in the set into a different vector. Here a vector is represented by a letter, an arrow at the top of the letter may be used to indicate the direction of a vector, *e.g., a *stands for either or . If = (1, 2, 3), then = (3, 2, 1), = (1, 2, 3, 3, 2, 1), and = (1, 2, 3, 1, 2, 3). Initially, each vector has only one position.

2. Create an initial distance matrix (*d*) whose elements are *d*(*u, v*) = *p*_{i,j}, where *i *is the position in vector *u, j *is the position in vector *v*, and *p*_{i,j }is the p-value of chi-square test described above.

3. Convert the distance matrix *d *to the transformed distance matrix *D*, whose elements are *D*(*u, v*), by the following conversion function (see [26]):

Where *V *is the number of vectors under consideration, and its value decreases by 1 in each iteration.

4. Find the minimum *D*(*u, v*) in *D*. Then a new vector *x *is formed by joining vector *u *and *v *according to **Algorithm 1 **[see Additional file 1] for a *k*^{th }order Markov chain.

5. Update the matrix *d *by replacing *u *and *v *by *x*. The distance of *x *to other remaining vector *y *is defined by:

*d*(*x, y*) = (*d*(*u, y*) + *d*(*v, y*) - *d*(*u, v*))/2

6. Go back to step 3 if the number of vectors in *C*_{s }is larger than 2, otherwise join the last two vectors according to **Algorithm 1**.

The order of positions in the final vector is the optimized linear chain for Markov model. Joining the first position to the last position in the vector forms a circular chain. A linear chain could be further optimized by forming a circular chain first from the final vector, then breaking the circular chain between positions with the weakest dependency, *e.g*., between positions *i *and *j *where *p*_{i,j }is the largest or the log-likelihood of the corresponding linear chain model is maximized. DNJ not only optimizes position order for linear chain models but also improves circular chain models, particularly when the order of Markov model is low, *e.g*., 1^{st }or 2^{nd }order Markov models.

Many different mixtures of Markov models can be formed from the combination of different Markov chains. It is essential to choose the model that minimizes prediction error. In model selection, we first fit each model using maximum likelihood smoothed by a Dirichlet prior [see Additional file 1], then compute either the Akaike information criterion (AIC) [27] or the Bayesian information criterion (BIC) [28]. The model with the minimum value of AIC or BIC is selected as the potential best model. Minimizing AIC is the same as choosing the model with the minimum prediction error or loss, while minimizing BIC is equivalent to choosing the model with the largest posterior probability. Nonetheless, AIC and BIC have a similar form:

-2·*loglik *+ *λ*·*DF*

where *λ *= 2 for AIC and *λ *= log(*N*) for BIC (*N *is the number of sequences); *loglik *is the maximized log likelihood of data given the model; DF is the degrees of freedom (number of free parameters). We replace DF with the effective degrees of freedom (EDF) in calculating AIC or BIC of the mixture of Markov models, which enables an appropriate model to be selected (see sub-section *Effective degrees of freedom*). There is no clear better choice between AIC and BIC for model selection. AIC tends to choose a model too complex as *N *→ ∞, while BIC tends to choose a model too simple when *N *is small. In our test on 61 different TFBS datasets, whose sample sizes range from 20 to 130, we preferred AIC to BIC for picking an appropriate model.

Let B be the set of bases (|B| denotes the number of different bases in B), *e.g*., for DNA sequences B = {*A, C, G, T*} (|B| = 4). For a motif of length *w*, the DF for a *k *order Markov model is (|B|^{k }- 1) × (*w *- *k*) for a linear Markov chain; and (|B|^{k }- 1) × *w *for a circular chain model. That is, the DF increases exponentially as the order of Markov chain increases. As a result, AIC or BIC often pick a simpler mixture model than the best model, especially when |B| is large. Tested on 61 human regulatory motifs from the Transfac database (ver. 7.4) [29], we found that both AIC and BIC selected the 0^{th }order Markov models for all 61 DNA regulatory motifs when using the DF. To avoid picking overly simple models, we used the EDF described below to calculate AIC and BIC.

Generally, only a subset of bases from B appears in a particular position of a set of biological motifs. The more conserved a position, the fewer bases are in the subset. The EDF for a model is related to the observed bases in training samples. For example, suppose that one would like to estimate nucleotide frequencies occurring in a position in a set of DNA training motifs. If only base A is observed in the position, then one needs to estimate only the frequency of A, the remaining parameters, *i.e*., the frequencies of C, G, T, can be derived from any prior information. Therefore, the actual DF is one in this case. For our mixture of Markov models, the EDF is defined as the number of parameters that are direct estimates of the observed bases in a training motif set. Let *b*_{i }be the base set observed in a position *i *of a training set of motifs. Additionally, let *h*^{k }be the sequence of motif positions in the *k*^{th }order Markov chain, be the motif position in the *i*^{th }element of *h*^{k}, and ∑|*h*^{k}| = *w *( |*h*^{k}| is the number of positions in *h*^{k }), then we define the EDF for the *k*^{th }Markov chain as

where if *i *- *k *≤ 0; and are for linear and circular chains, respectively. The total EDF for a mixture Markov model is just the sum of EDFs of all individual chains. For example, the total EDF for the special 0-k mixture model equals to the EDF sum of the 0^{th }and the *k*^{th }order chains.

We test the effectiveness of our method on TFBS data and the donor splice sites, where training data for OMiMa are a set of sequences of a motif. For prediction results, we use the following abbreviations for empirical quantities: *TP *(# true positives), *TN *(# true negatives), *FP *(# false positives), *FN *(# false negatives), *Ac *(Accuracy), *Sn *(sensitivity), *Sp *(specificity), and *Mc *(Matthews correlation coefficient). *Sn, Sp, Ac*, and *Mc *are defined as:

Matthews correlation coefficient [30], also called Phi (correlation) coefficient, has a value between -1 and 1, with 1, 0, and < 0 indicating a perfect prediction, a random prediction, and a worse than random prediction, respectively.

OMiMa can use two ways to score a motif site **x**: log-likelihood and log-likelihood ratio, which are defined by

where *M*_{s }is the signal model trained by true motif sites, and *M*_{b }is the background model or false signal model trained by background sequences or false motif sites. A sequence **x** is predicted as a positive site if the score of **x** is larger than a certain threshold. We select a cutoff threshold using one of the following three criteria: balanced sensitivity and specificity, the maximum prediction accuracy, and the maximum Matthews correlation coefficient. Each potential threshold yields an estimated true positive rate and a false positive rate. The plot of true positive rates against false positive rates generates a Receiver Operating Characteristic (ROC) curve, which can be used for comparing and selecting the best model.

We used a three-symbol notation 'k-m-s' to distinguish different models, where 'k' stands for a 0-k mixture Markov model, 'm' is either 'L' or 'C' to indicate whether the *k*^{th }order chain is linear ('L') or circular ('C'), and 's' is either 0 or 1 to indicate whether log likelihood score (0) or log-likelihood ratio score (1) is used. For example, '1-L-l' stands for a 0–1 mixture of linear Markov models that uses log-likelihood ratio to score a motif site.

To assess the ability of our DNJ method for optimizing a Markov chain, we compared the DNJ method with random permutation method. In this evaluation, we used a 0-k mixture model (*k *= 0, 1, 2) (Figure (Figure2)2) to model transcription factor binding sites from the Transfac database. For each TFBS, we first fitted a 0-k mixture model (denoted as *M*_{DNJ}) with its *k*^{th }order Markov chain optimized by the DNJ method. We calculated the log-likelihood of the data given the model *M*_{DNJ }(log Pr(*data*|*M*_{DNJ})). Second, with the same data, we fitted a new 0-k mixture model (denoted as *M*_{R}), which is the same as *M*_{DNJ }except that the positions in its *k*^{th }order chain are ordered by random permutation, and calculated logPr(*data*|*M*_{R}). This step was repeated 1,000 times, so we have 1,000 log-likelihoods of the randomly permuted models . We then calculated the empirical *p_value *of the DNJ-optimized model as follows:

where *δ *is an indicator function with value 1 if condition is true, and 0 otherwise. The smaller the *p_value*, the better the DNJ optimization is; and *p_value *= 0 means the DNJ-optimized model performs better than any one of the 1,000 randomly permuted models. The *p_value *approximates the probability of observing log Pr(*date*|) larger than (log Pr(*data*|*M*_{DNJ}).

Fifty-three human transcription factors, whose binding sites contain at least four dependent positions by the *χ*^{2 }test given by equation (3), are selected for this evaluation (Table (Table1).1). The assessment was performed on four 0-k mixture models: 1^{st }order linear chain, 1^{st }order circular chain, 2^{nd }order linear chain, and 2^{nd }order circular chain.

The optimized 1^{st }order Markov chains for TFBS. The optimized arrangement of dependent positions within TFBS for the 1^{st }order Markov model. N and *N*_{D }are total number of motif positions and the number of positions significantly dependent, respectively. **...**

Results suggest that DNJ method performed remarkably well in optimizing the 1^{st }order linear Markov chains, that in 49 out 53 cases, the DNJ optimized models were the best or close to the best (Figure (Figure4a).4a). The optimization for the 2^{nd }order linear chains was slightly worse than that for the 1^{st }order linear chains, partially because the DNJ method relies only on the pairwise dependencies between two single positions. Nevertheless, most of the DNJ optimized models were still close to the best [see Additional file 1 Figure 1a]. Although our DNJ method was designed for optimizing linear Markov chains, it still worked well in optimizing the 1^{st }order circular Markov chains (Figure (Figure4b).4b). However, the DNJ method did not perform well in optimizing the 2^{nd }order circular Markov chains [see Additional file 1 Figure 1b].

We used AP1 (activating protein 1) transcription factor binding sites (Transfac ID V$AP1_Q4_01) as an example of how DNJ optimization can improve performance of a 0–1 or 0–2 mixture model. We plotted the histogram of the log-likelihood per instance given a model , log Pr(*data*|)/*N*, where *N *is the number of sequences in the data set, and *i *= 1,..., 1000 for 1,000 mixture models of randomly permuted Markov chains. The histogram represents a simulated null-distribution of log-likelihood per instance given a mixture model. Then we mapped the location of the likelihood per instance given DNJ optimized model, (log Pr(*data*|*M*_{DNJ})/*N*, in the histogram. For the transcription factor V$AP1_Q4_01, we found that for the 0–1 mixture model of either linear or circular structure, DNJ optimized models are better than any models from 1,000 random permutations (Figure (Figure55).

Theoretically, the optimal model can be found by exhaustively searching through all possible models. An exhaustive search is not always possible in practice, however, as the search space can be very large. The number of possible Markov chains is the factorial of the length of the Markov chain and dramatically increases as the length of chain increases. For example, the computational time for a motif of 15 bases (15! = 1.307674e + 12) can be practically unacceptable. Our DNJ method can deal with such long motifs because of its computational efficiency.

One interesting application of our mixture model is TFBS identification. In this assessment, we used a couple of examples to show how OMiMa can improve prediction accuracy when there are position dependencies within a TFBS. We first tested our method on simulated data where the exact dependency structure of a TFBS is known. We tested whether OMiMa can capture such dependency and optimize the Markov model accordingly. Next, we tested our method on real motif data for AP1. In both examples, we compared OMiMa performance to PWM, PVLMM, and the 1^{st }order Markov model (1stMM) with its motif positions in the natural order. PVLMM, run on Microsoft Windows, is based on the variable length Markov model (VLMM) [31,32]. Except for its order and depth parameters, PVLMM was run under its default settings in all comparisons.

Many TFBS are palindromic sites bound by heterodimers/homodimers (*e.g*.,: Jun-Fos, Myc-Max, Max-Max and p50-p50). The sequences in two half sites of a palindromic TFBS are usually not perfectly complementary, and the strong binding to one half site may compensate for weak binding to the other. We simulated two imperfect palindromic TFBS (named A and B) of length 12 bases. For each TFBS, the bases in each position were generated from the uniform distribution (the frequency of each base is 0.25). The base in one half site and its reverse complementary base in the other half were generated using the probabilities listed in Table Table2.2. Therefore, there are only pairwise position dependencies in the simulated TFBS. The position pair 0–11 has the strongest dependency, whereas the pair 5–6 has the weakest dependency (they are independent). Overall, motif A has stronger position dependencies than motif B. The false sites of TFBS were simulated from the uniform distribution of four nucleotides without any constraints of base pairing between the two half sites. The simulated data of each TFBS consist of a training set with 150 true sites, and a testing set with 150 true sites and 150 false sites. Using these simulated training sets, OMiMa found all true dependencies significant at level *α *= 0.05 (see Table Table2).2). For both TFBS, OMiMa was also able to arrange the positions of each dependent pair to be the nearest neighbors in their 0–1 mixture models: (See table table77)

Simulation of two palindromic TFBS. Simulation of two palindromic TFBS A and B. The first 2 columns are the complementary positions of the palindromic TFBS. The 3^{rd }and 4^{th }columns are simulation parameters, which specify the probabilities of forming **...**

In our simulation, the positions 5 and 6 were generated independently from all other positions, so they should be in the 0^{th }order chains. However, based on OMiMa's pairwise *χ*^{2 }tests for the training data, the position pair 5–8 (with p-value = 0.03) in TFBS A, and the position pair 6–10 (with p-value = 0.04) in TFBS B were declared dependent. That is why the positions 5 and 8 were arranged together in the model for TFBS A, and positions 6 and 10 were together for TFBS B. We compared the prediction results of OMiMa's 0–1 mixture model with those of PWM, 1stMM and the 1^{st }order PVLMM (with depth 1). Results (Table (Table3)3) showed that OMiMa outperformed all other models, and PVLMM performed better than IstMM and PWM. Additionally, we used smaller training sets to access the performance of these methods on the same testing set. Smaller training sets, in sizes ranging from 15 to 150, were independently sampled (without replacement) from the original 150 sites for training. Results suggested that OMiMa performed consistently better than the other methods, regardless the size of a training set (Figure (Figure6)6) [see Additional file 1 Figure 2].

We chose human AP1 TFBS (see Figure Figure7a)7a) for this evaluation because of its relatively large number of known sites. In total, we had 119 true sites and 5950 false sites. The true sites were extracted from Transfac database (Transfac ID V$AP1_Q4_01), and false sites were randomly sampled from the non-coding regions of the human genome. Our *χ*^{2 }tests on the 119 true sites suggested that all positions showed some level of dependency with the neighboring pairs 0–2, 4–5, 5–6, and 4–6 showing strong dependencies (p-value < 1.0e-6). Noticeably, the positions 4, 5 and 6 are also the most conserved positions, so we expect that PWM would be reasonable good model for the TFBS. We randomly split both the true sites and false sites into 10 roughly equal-sized parts, and used a 10-fold cross validation to compare the performance of OMiMa's 0–1 mixture model with the others. OMiMa had advantage over the other three models in predicting TFBS that do not have strong long-range dependencies (Table (Table4).4). Results also showed the first order PVLMM did not perform better than 1stMM and PWM. We found that the first order PVLMM arranged the position pair 5–6, which showed the strong dependency, differently from OMiMa and 1stMM. In only 3 out 10 times, PVLMM arranged positions 5 and 6 as direct neighbors, while OMiMa did in 9 out 10 times, and 1stMM did naturally all times. This is one possible reason why PVLMM performed slightly worse.

The transcription of most higher eukaryotic genes involves RNA splicing, in which primary transcripts become mature mRNA by removing introns. The donor or 5' splice sites and the acceptor or 3' splice sites on the boundaries of exons and introns provides critical signals for precise splicing. Therefore, splice site recognition has been widely used by gene finding tools such as GENESCAN [2] and GENIE [33] for gene prediction. The splicing process starts with U1 snRNP binding to the donor site via base-pairing of U1 snRNA. The base pairing between U1 snRNA and the donor site, however, need not be perfectly complementary in all positions [34,35]. Both experimental and computational evidence suggest that there are mutually dependent positions within the donor site: a mis-matched pair of U1 snRNA and the donor sites at one position can be compensated for by a matching base pair at another position, and vice versa [2,24,36,37]. Modeling such dependency structure within the donor site has been used to improve donor site prediction [2,23,24,33]. We used two independent datasets of human donor sites to assess the performance of OMiMa in comparison with leading competitors.

The test dataset of human donor splice sites (Reese data) was from [38]. This dataset has 6246 donor sites (1324 real and 4922 false) of length 15 bases from -7 to +8 around the conserved 'GT' dinucleotide. The dataset consists of a training set (containing of data) and a testing set (the remainder), which were previously used to assess the performance of NNSplice [33]. We used the same partitions for training and testing in the following comparisons.

First, we tested whether OMiMa which uses either AIC or BIC, can correctly pick the best model based on ROC analysis. We fitted a set of 0-k mixture models, in which the *k*^{th }order chains are either linear or circular and *k *ranges from 0 to 3, with the training data. We subsequently applied the fitted models to predict splice sites in testing data. The performances of different models were compared and evaluated by ROC analysis (Figure (Figure8).8). In addition, we compared the maximum accuracy (*Ac*) and the maximum Matthews correlation efficient (*Mc*) achieved by each model (data not shown). The best models were 0–1 mixture models (Figure (Figure8).8). Both the linear and circular models performed about the same. The best models picked by ROC analysis are consistent with those selected by OMiMa (AIC criterion). The selected models were further confirmed by a six-fold cross validation.

Using the best model selected above, we then compared OMiMa with NNSplice and PVLMM. NNSplice is based on a complex neural network model and is trained by both true sites and false sites. Since both OMiMa and NNSplice used the same training and testing data, their prediction results can be directly compared. We compared OMiMa's 1-L-1 and 1-C-1 models with the first order PVLMM (with depth 1) as all have similar model complexity. The results of NNSplice were reported at the NNSplice Web site [39]. We found that OMiMa had comparable prediction accuracy to NNSplice and PVLMM (Table (Table5).5). In addition, OMiMa is much more computationally efficient than NNSplice and PVLMM [see Additional file 1].

Given enough training data, we can use more complicated models than the 0–1 mixture model to improve prediction accuracy. In this evaluation, we test whether 0-k mixture models can compete with the MEM on a much larger dataset. This large donor site dataset (Yeo data), used to assess performance of the MEM, was from [40]. The dataset, extracted from 1821 non-redundant human transcripts, has 8,415 real and 179,438 decoy sites in the training set, and 4,208 real and 89,717 in decoy sites in the testing set. Each real site has length 9 bases from -3 to +6 around the conserved 'GT' of donor splice sites recognized by U-2 type spliceosome. The decoy sites are any other sequence segments in the exons and introns matching the pattern *N*_{3}*GTN*_{4}. So a decoy site can have the exactly same sequence as a real site. We applied this original training and testing sets to assess performance of OMiMa, where we used only log-likelihood ratio scoring. In addition, we ran a 3-fold cross-validation, in which the number of sites in new training and testing sets are roughly the same as those in the original ones [see Additional file 1 Table 2]. The top 4 models selected by AIC are 3-C-1, 3-L-1, 2-C-1 and 2-L-1, respectively, consistent with the ROC analysis. To find the top 4 sub-models of PVLMM by ROC analysis, we used a series of Markov orders and/or depths (1 ≤ *order *≤ 4 and *order *≥ *depth*) to predict the same data sets. For convenience, we use notation "P:*k*-*d*" to denote a PVLMM of order *k *and depth *d*. We adopt notation in [23] for sub-models of MEM.

Briefly, the notation has the form "me*K*s*D*" or "me*K*x*D*", where "me" stands for maximum entropy; "*K*" is a number for the marginal order or the maximum length of an oligomer in consideration; "*D*" is the skip number or maximum skip number determining which positions the bases of an oligomer are from; "s" stands for skip number and "x" for the maximum skip. For example, model "me5s0" considers all marginal distributions of *p*(*x*_{i}), *p*(*x*_{i}, *x*_{i+1}), *p*(*x*_{i}, *x*_{i+1}, *x*_{i+2}), *p*(*x*_{i}, *x*_{i+1}, *x*_{i+2}, *x*_{i+3}), *p*(*x*_{i}, *x*_{i+1}, *x*_{i+2}, *x*_{i+3}, *x*_{i+4}).

Comparison of the top 4 performers from each model class suggested that OMiMa performed comparably with MEM and better than PVLMM (Table (Table6)6) (all models' performances suffered on this data set because about 98% real sites appeared at least once in the decoy set). One advantage of OMiMa over MEM is that, for the models with similar performance, OMiMa's models generally have fewer parameters and thus require fewer training samples. Our test showed that OMiMa was able to retain similar performance of MEM even when trained by only 60% data of MEM's original training sets [see Additional file 1 Table 3].

To compare OMiMa's fitted donor site models to biological knowledge about dependencies among positions, we examined the best donor models for the first donor dataset (Reese data) and for the second donor dataset (Yeo data). For convenience, let us mark the invariant 'GT' nucleotides in the boundary of exon/intron as the positions 0 and 1 of the donor site, respectively (see Figure Figure7b).7b). First, based on 1,116 real donor sites in the Reese original training data, the 0–1 mixture model was selected as the best model with the following 1^{st }order chain.

-2 5 -1 3 4 -3 -7 -6 -5 -4 7 6 2

We found that this position arrangement is supported by the following biological evidence of base-pairing between U1 snRNA and the donor site: (a) 5'/3' compensation effect: a base pair at position -1 can prevent an aberrant splicing caused by a mis-matched pair at position 5 [37]; (b) Adjacent base-pair effect: a matching base pair at position 3 is rare in the absence of a matching base pair at position 4 [2,41]; (c) A matching base-pair at the non-conserved positions 6 and 7 can compensate for a mis-matched pair at position 2 [42]. Interestingly, the model also arranged non-conserved positions (-4, -5, -6, -7, 6, 7) together as it did for the other more conserved positions. Second, based on 8,415 real donor sites of the Yeo original training sites, the 0–3 mixture model was the best model. The optimized position order of its 3^{rd }order chain was:

2 5 -1 4 -2 3 -3

We can see that this model is consistent with the above evidence (a) and (b). In addition, it is well supported by experimentally verified position dependencies of position 4 on the positions -1, -2, 3 and 5 [37], and the computationally confirmed dependency of position -3 on position -2 due to the adjacent base-pair effect [2].

The prediction accuracy of a probabilistic model is largely determined by the effectiveness of the model in characterizing a biological motif. Since there is large variation of the signals embedded in biological motifs, an effective model can be as simple as a consensus sequence or as complex as a fully connected network model. In this paper, we described a mixture of Markov models to allow adjustment of model complexity for different motifs. Also, we extended the traditional linear chain Markov model to the circular chain Markov model, which can better represent position dependencies within a motif in some cases. We presented a novel method, DNJ, for efficiently optimizing position arrangement of a non 0^{th }order Markov chain to incorporate most dependencies. We described methods for calculating the EDF and for selecting the best mixture Markov model. We implemented these methods in our motif finding OMiMa system, which is freely available. Finally, we demonstrated from different aspects in several examples that OMiMa can improve motif prediction accuracy in biological sequences.

The interaction of biological macromolecules, such as transcription factors bound to DNA sites, usually involves several highly dependent positions functioning as a unit. Many methods including Markov chains, Bayesian trees, and neutral networks have been used to model dependency structures within a motif. The Markov model is the simplest yet can be very powerful when it is optimized. Our results showed that the optimized Markov models performed better than the neural network model and PVLMM, and comparably with MEM for splice site prediction. The optimized Markov model can incorporate both local and non-local dependencies into the model, which enables it to compete with tree or network models in predicting short biological motifs. We also showed that the optimized Markov model can be an excellent motif predictor. Moreover, it is also computationally efficient due to its simplicity.

Model complexity, measured by parameter number, is an important issue in motif modeling. The more complex a model, the more data are needed for adequate training. For many biological motifs, however, the number of known (experimentally determined) sites is small. This limits the usage of complex models, such as higher order Markov models, Bayesian trees, network models or MEM, even though these models in some cases can perform better than the simpler models given enough training data. For a standard Markov model, the number of its parameters increases exponentially as its Markov order increases. Without sufficient training data, it is difficult to accurately estimate all model parameters, even using more robust methods (*e.g*. interpolated Markov chains [43,44]). As a result, lack of sufficient training data often makes it impractical to train a higher order Markov model. On the other hand, a low order Markov model may perform poorly by failing to incorporate more distant dependencies. Several motif models and methods have been developed to address this issue. One of these models is the variable length Markov model (VLMM), whose Markov orders (also called context lengths) can vary among different positions. VLMM can effectively reduce Markov model complexity when the variation of actual context lengths is large. VLMM, however, is not the best choice to incorporate long-range dependencies. The position optimized Markov model (POMM) [21] is able to incorporate important distant dependencies without increasing Markov chain order. However, the effectiveness of this model largely depends on the optimization routine.

More recently, Zhao *et al. *[24] described the PVLMM in an attempt to combine advantages of both VLMM and POMM. The disadvantage of PVLMM is that the number of possible permutations is the factorial of motif length, which makes it more computationally expensive. In addition, the random permutation method used by PVLMM for optimization is more likely to overfit the model, *e.g*., incorporating non-significant dependencies into the PVLMM model that can reduce its prediction power. The optimized mixture of Markov models we presented here tries to inherit advantages of these existing models while avoiding their disadvantages. In OMiMa, we replace VLMM with a mixture of several lower order Markov models, which are subsequently optimized to account for long-range dependencies.

In comparison with other leading methods, OMiMa can incorporate more than the NNSplice's pairwise dependencies; OMiMa avoids model over-fitting better than the PVLMM; and OMiMa requires smaller training samples than the MEM. These are primarily reasons that OMiMa showed superior performance, in terms of prediction accuracy, required size of training data or computational time, over other leading-methods in our results.

With any model selection procedure, the possibility of choosing a model that drastically over- or underfits is a concern. OMiMa employs AIC and BIC, two standard criteria, that are widely used because they tend to avoid extreme over- or underfitting. Both have theoretical support [27,28]. In our application, neither criterion worked well when using the DF (results not shown); but both, particularly AIC, performed well when using EDF. We found that models selected by AIC using EDF were consistent with models selected by cross-validation and by ROC analysis.

Our OMiMa approach has two features that can be limitations when the size of the training data is small. First, the chi-square test that partitions motif positions into those with dependencies and those without dependencies will, like any statistical test, make mistakes, and its statistical power to detect dependencies will suffer with small training samples. Although the test will not always provide a correct partition, our approach should adapt to strong or weak dependencies overall and improve prediction when dependencies are strong. In addition, weakly dependent positions mistakenly placed in the set with no dependencies are often adequately modeled by a 0^{th }order chain, whereas independent positions mistakenly assigned to the set with dependencies will be placed by the DNJ algorithm in locations with the least impact on the *k*^{th }order chain. Second, the EDF that we used in model selection is an estimate based on the training data. For degenerate sites, the estimate should be accurate with even small training samples; whereas for conserved sites a larger training sample might reveal additional bases and change the EDF. Still, such additions should be minimal and would generally induce small changes in the EDF, so we expect little impact on model selection. Any methods that employ chi-square techniques to test for dependent sites face similar limitations. Nevertheless, OMiMa with its relatively small parameter space should adapt to small training datasets better than many competitors. Of course, any motif finding algorithm would do better with larger training samples.

OMiMa places no limit on the length of sequences that it can scan, and it could be used to find TFBS in any sequenced organism as long as a training motif set is available. The larger the genome evaluated, the more false positives are likely to be declared. Although OMiMa's prediction accuracy will help, other approaches to reducing false positives will be needed. Cross-species comparisons and relative location compared to transcription start sites have been used to reduce false positives and could be used with OMiMa too. Furthermore, OMiMa's ability to accurately and quickly identify splice sites should be easy to incorporate into probabilistic gene-prediction programs where correct prediction of splice sites is critical.

Our optimized mixture of Markov models represents an alternative to the existing methods for modeling dependent structures within a biological motif. Unlike existing methods, our model is conceptually simple and effective, which has advantages in a large scale motif prediction. In particular, with its ability to minimize model complexity, our method can work effectively even with limited training data. The optimized mixture of Markov models is implemented in our computational tool OMiMa, which can use a variety of mixture models for motif prediction. OMiMa, in which most parameters are configurable, is freely available to all users.

W. Huang provided the principal contributions to the conception and design of this study as well as to its analysis. D. M. Umbach and L. Li contributed to the design of the study and the interpretation of results. All authors contributed to writing and critically revising the manuscript.

The supplement includes the mathematical formulas for computing the probability of a motif site given a Markov model, the algorithmic pseudo-code for the DNJ method, and the description of the parameter estimation for our model. It also contains supplemental materials for the main results as well as other additional results, such as the application for protein domain identification, the comparison of computational time, and so on.

Click here for file^{(312K, pdf)}

We thank Drs Bruce Weir and Jeffrey Thorne for critically reading the manuscript, and Drs Clarice Weinberg and Joseph Nevins for helpful comments. This research was supported by Intramural Research Programs of the NIH, National Institute of Environmental Health Sciences.

- Weichun Huang's Research Domain http://BioMedEmpire.org
- Burge C, Karlin S. Prediction of complete gene structures in human genomic DNA. J Mol Biol. 1997;268:78–94. [PubMed]
- Wray GA, Hahn MW, Abouheif E, Balhoff JP, Pizer M, Rockman MV, Romano LA. The Evolution of Transcriptional Regulation in Eukaryotes. Mol Biol Evol. 2003;20:1377–1419. [PubMed]
- Kellis M, Patterson N, Endrizzi M, Birren B, Lander ES. Sequencing and comparison of yeast species to identify genes and regulatory elements. Nature. 2003;423:241–254. [PubMed]
- Negre B, Casillas S, Suzanne M, Sanchez-Herrero E, Akam M, Nefedov M, Barbadilla A, de Jong P, Ruiz A. Conservation of regulatory sequences and gene expression patterns in the disintegrating Drosophila Hox gene complex. Genome Res. 2005;15:692–700. [PubMed]
- Xie X, Lu J, Kulbokas EJ, Golub TR, Mootha V, Lindblad-Toh K, Lander ES, Kellis M. Systematic discovery of regulatory motifs in human promoters and 3' UTRs by comparison of several mammals. Nature. 2005;434:338–345. [PMC free article] [PubMed]
- Staden R. Computer methods to locate signals in nucleic acid sequences. Nucl Acids Res. 1984;12:505–519. [PMC free article] [PubMed]
- Stormo GD, Fields DS. Specificity, free energy and information content in protein-DNA interactions. Trends Biochem Sci. 1998;23:109–113. [PubMed]
- Quandt K, Freeh K, Karas H, Wingender E, Werner T. Matlnd and Matlnspector: new fast and versatile tools for detection of consensus matches in nucleotide sequence data. Nucl Acids Res. 1995;23:4878–4884. [PMC free article] [PubMed]
- Kel AE, Gössling E, Reuter I, Cheremushkin E, Kel-Margoulis OV, Wingender E. MATCH: A tool for searching transcription factor binding sites in DNA sequences. Nucleic Acids Res. 2003;31:3576–3579. [PMC free article] [PubMed]
- Agarwal P, Bafna V. RECOMB '98: Proceedings of the second annual international conference on Computational molecular biology. New York, NY, USA: ACM Press; 1998. Detecting non-adjoining correlations with signals in DNA; pp. 2–8.
- Man TK, Stormo GD. Non-independence of Mnt represser-operator interaction determined by a new quantitative multiple fluorescence relative affinity (QuMFRA) assay. Nucl Acids Res. 2001;29:2471–2478. [PMC free article] [PubMed]
- Benos PV, Lapedes AS, Fields DS, Stormo GD. SAMIE: statistical algorithm for modeling interaction energies. Pac Symp Biocomput. 2001:115–26. [PubMed]
- Bulyk ML, Johnson PLF, Church GM. Nucleotides of transcription factor binding sites exert interdependent effects on the binding affinities of transcription factors. Nucl Acids Res. 2002;30:1255–1261. [PMC free article] [PubMed]
- Roulet E, Busso S, Camargo AA, Simpson AJG, Mermod N, Bucher P. High-throughput SELEX SAGE method for quantitative modeling of transcription-factor binding sites. Nat Biotechnol. 2002;20:831–835. [PubMed]
- Krivan W, Wasserman WW. A Predictive Model for Regulatory Sequences Directing Liver-Specific Transcription. Genome Res. 2001. p. GR1806R.
- Schneider TD, Stormo GD, Gold L, Ehrenfeucht A. Information content of binding sites on nucleotide sequences. J Mol Biol. 1986;188:415–431. [PubMed]
- Zhang MQ, Marr TG. A weight array method for splicing signal analysis. Comput Appl Biosci. 1993;9:499–509. [PubMed]
- Ponomarenko MP, Ponomarenko JV, Frolov AS, Podkolodnaya OA, Vorobyev DG, Kolchanov NA, Overton GC. Oligonucleotide frequency matrices addressed to recognizing functional DNA sites. Bioinformatics. 1999;15:631–643. [PubMed]
- Cai D, Delcher A, Kao B, Kasif S. Modeling splice sites with Bayes networks. Bioinformatics. 2000;16:152–158. [PubMed]
- Ellrott K, Yang C, Sladek FM, Jiang T. Identifying transcription factor binding sites through Markov chain optimization. Bioinformatics. 2002;18:S100–S109. [PubMed]
- Barash Y, Elidan G, Friedman N, Kaplan T. RECOMB '03: Proceedings of the seventh annual international conference on Computational molecular biology. New York, NY, USA: ACM Press; 2003. Modeling dependencies in protein-DNA binding sites; pp. 28–37.
- Yeo G, Burge CB. Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals. J Comput Biol. 2004;11:377–394. [PubMed]
- Zhao X, Huang H, Speed TP. RECOMB '04: Proceedings of the eighth annual international conference on Computational molecular biology. New York, NY, USA: ACM Press; 2004. Finding short DNA motifs using permuted markov models; pp. 68–75.
- Zhou Q, Liu JS. Modeling within-motif dependence for transcription factor binding site predictions. Bioinformatics. 2004;20:909–916. [PubMed]
- Saitou N, Nei M. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987;4:406–425. [PubMed]
- Akaike H. A new look at the statistical model identification. IEEE Trans Automat Control. 1974;19:716–723.
- Schwarz G. Estimating the dimension of a model. Ann Stat. 1978;6:461–464.
- Wingender E, Chen X, Hehl R, Karas H, Liebich I, Matys V, Meinhardt T, Prüâ M, Reuter I, Schacherer F. TRANSFAC: an integrated system for gene expression regulation. Nucleic Acids Res. 2000;28:316–319. [PMC free article] [PubMed]
- Matthews B. Comparison of the predicted and observed secondary structure of T4 phage lysozyme. Biochim Biophys Acta. 1975;405:442–451. [PubMed]
- Rissanen J. Complexity of strings in the class of Markov sources. IEEE Trans Inform Theory. 1986;32:526–532.
- Bühlmann P, Wyner AJ. Variable length Markov chains. Ann Statist. 1999;27:480–513.
- Reese MG, Eeckman FH, Kulp D, Haussler D. Improved splice site detection in Genie. J Comput Biol. 1997;4:311–323. [PubMed]
- Ketterling RP, Drost JB, Scaringe WA, Liao DZ, Liu JZ, Kasper CK, Sommer SS. Reported in vivo splice-site mutations in the factor IX gene: severity of splicing defects and a hypothesis for predicting deleterious splice donor mutations. Hum Mutat. 1999;13:221–231. [PubMed]
- Staley JP, Guthrie C. An RNA switch at the 5' splice site requires ATP and the DEAD box protein Prp28p. Mol Cell. 1999;3:55–64. [PubMed]
- Thanaraj T, Robinson AJ. Prediction of exact boundaries of exons. Brief Bioinform. 2000;1:343–356. [PubMed]
- Carmel I, Tal S, Vig I, Ast G. Comparative analysis detects dependencies among the 5' splice-site positions. RNA. 2004;10:828–840. [PubMed]
- Berkeley Drosophila Genome Project http://www.fruitfly.org/seq_tools/datasets/Human/GENIE_96/splicesets
- BDGP: Splice Site Prediction by Neural Network http://www.fruitfly.org/seq_tools/splice.html
- Christopher Burge Lab http://genes.mit.edu/burgelab/maxent/ssdata
- Nelson K, Green M. Mechanism for Cryptic Splice Site Activation During Pre-mRNA Splicing. PNAS. 1990;87:6253–6257. [PubMed]
- Nandabalan K, Price L, Roeder GS. Mutations in U1 snRNA bypass the requirement for a cell type-specific RNA splicing factor. Cell. 1993;73:407–415. [PubMed]
- Salzberg S, Delcher A, Kasif S, White O. Microbial gene identification using interpolated Markov models. Nucl Acids Res. 1998;26:544–548. [PMC free article] [PubMed]
- Ohler U, Harbeck S, Niemann H, Noth E, Reese M. Interpolated markov chains for eukaryotic promoter recognition. Bioinformatics. 1999;15:362–369. [PubMed]
- Crooks GE, Hon G, Chandonia JM, Brenner SE. WebLogo: A Sequence Logo Generator. Genome Res. 2004;14:1188–1190. [PubMed]

Articles from BMC Bioinformatics are provided here courtesy of **BioMed Central**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's Canada Institute for Scientific and Technical Information in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |