Transcription factor binding sites (TFBSs) are crucial in the regulation of gene transcription. Recently, chromatin immunoprecipitation followed by cDNA microarray hybridization (ChIP-chip array) has been used to identify potential regulatory sequences, but the procedure can only map the probable protein-DNA interaction loci within 1–2 kb resolution. To find out the exact binding motifs, it is necessary to build a computational method to examine the ChIP-chip array binding sequences and search for possible motifs representing the transcription factor binding sites.
We developed a program to find out accurate motif sites from a set of unaligned DNA sequences in the yeast genome. Compared with MDscan, the prediction results suggest that, overall, our algorithm outperforms MDscan since the predicted motifs are more consistent with previously known specificities reported in the literature and have better prediction ranks. Our program also outperforms the constraint-less Cosmo program, especially in the elimination of false positives.
In this study, an improved sampling algorithm is proposed to incorporate the binomial probability model to build significant initial candidate motif sets. By investigating the statistical dependence between base positions in TFBSs, the method of dependency graphs and their expanded Bayesian networks is combined. The results show that our program satisfactorily extract transcription factor binding sites from unaligned gene sequences.
Understanding transcription is central to understanding genetic regulatory mechanisms. The transcription of a gene is generally dependent on the presence of specific signals located at upstream regions of the core-promoter. These specific signals derive from their use as binding sites by transcription factors (TFs), and are therefore termed transcription factor binding sites (TFBSs). Recently, chromatin immunoprecipitation followed by cDNA microarray hybridization (ChIP-chip array) has been used to identify potential regulatory sequences, but the procedure can only map the probable protein-DNA interaction loci within 1–2 kilobases resolution . To find out the exact binding motifs, it is necessary to build a computational method to examine the ChIP-chip array binding sequences and search for possible motifs representing the TFBSs (motif discovery).
There are many computational TFBS motif finding tools available [2-4]. The traditional approach for finding TFBSs is to collect and align a set of promoter sequences of co-regulated genes from either the literature or systematic experiments. Numerous computational tools, such as CONSENSUS , EM , MEME  and the Gibbs sampler , have utilized the approach to identify short DNA sequence motifs which are statistically over-represented in the promoter sequences.
Other than the alignment-based motif finding algorithms in above, many approaches have tried to extend to the use of evolutionary conservation information such as phylogenetic footprinting or the detection of combinations of binding sites (termed as cis-regulatory modules; CRMs) [2,3]. Phylogenetic footprinting methods [9-11] is an approach that seeks to identify conserved regulatory elements by comparing genomic sequences between related species. However, due to the statistical nature of the approach, e.g., a small amount of closely related species, not all transcription binding sites can be found by using phylogenetic footprinting. Hence, some algorithms have emerged to combine the alignment-based motif prediction with phylogenetic footprinting such as PhyloGibbs  and MY sampler . On the other hand, by the detection of CRMs due to the cooperative interactions between TFs, algorithms like those in [14-16] can produce predictions of substantially better specificity than those of isolated sites.
Recently, more effective motif finders, e.g., MDscan , ANN-Spec , DMOTIFS , DME  and Cosmo , have taken the advantage of a background set, serving as a negative control. The goal of these discriminant motif finders is to search only for motifs that are most discriminating, that is, only those enriched in the foreground set relative to the background set . Although these motif finders have improved the performance of TFBS prediction, it is still a trouble to have a satisfactory solution. How to find out accurate binding motifs may require much attention in the computational biology community. In this study, an improved sampling algorithm is proposed to incorporate the binomial probability model to build significant initial motif sets. By investigating the statistical dependence between base positions in TFBSs, it appears feasible to use statistical models to formulate the structural dependence of a motif in the identification of TFBSs. In light of this observation, the method of dependency graphs and their expanded Bayesian networks  is combined and prediction results show that our algorithm is able to find out motifs more consistent with previously known evidence.
Let TF be one of the transcription factors to be investigated. The binding dataset of the transcription factor TF, denoted as BTF, consists of the sequences with low binding p-value (< 0.001) to the TF in the ChIP-chip array data . A sliding window of size w is used to extract segments of length w when sliding through each of the sequences in BTF.
Let STF be the collection of all extracted segments from BTF, M the number of sequences in the binding dataset BTF, Li the length of the ith sequence in the binding dataset BTF, and TTF the total number of segments in STF. Then
To discover the binding motifs of the transcription factor TF, a number of initial candidate motif sets for TF is subsequently built from the collection STF of extracted segments. Note that the contents of segments, called patterns, in STF may not be distinct.
Most of early motif finding algorithms, such as Gibbs sampler  and MEME , have a weakness, where initial candidate motif sets are built by randomly extracting segments from sequences in the binding dataset BTF (i.e. randomly selecting segments from the set STF). To improve the deficiency, the binomial probability distribution model is firstly utilized in the establishment of a number of initial candidate motif sets in our algorithm.
Then in the process of iterative sampling in our algorithm to expand and/or trim each of the initial candidate motif sets, the method of dependency graphs and their expanded Bayesian networks  is used to develop a statistical model for the background motif set identified as the union S = ∪TFSTF of segments extracted from all transcription factor binding datasets.
The basic procedure to find the binding motifs of the transcription factor TF is as follows:
1. Build N initial candidate motif sets.
(a) Take N distinct patterns from the set STF with the most highest significance scores as the candidates by the binomial distribution model (see the Binomial probability distribution model subsection).
(b) Then for each of the N significant binding site candidates for the transcription factor TF, in view of evolution, collect all segments in STF whose patterns have no more than d Hamming distance matching to the candidate pattern to form an initial candidate motif set.
At this stage, N initial candidate motif sets for the transcription factor TF are built.
2. Iteratively sample through the binding dateset BTF to expand and/or trim each of the N initial candidate motif sets so that their approximate maximum a posteriori (AMAP) scores [1,24] can keep increasing until the N candidate motif sets are invariant in K consecutive iterations (see the Iterative sampling subsection).
(a) In the calculation of AMAP scores in this stage, the background model for the background motif set S = ∪TFSTF is established under the method of dependency graphs and their expanded Bayesian networks (see the Method of dependency graphs and their expanded Bayesian networks subsection).
3. Refine each of the N candidate motif sets by re-examining all the segments already included in the motif set. A segment is removed from the motif set if doing so increases the AMAP score.
A simple flowchart for our algorithm is shown in Figure 1. The following subsections will expatiate on each stage of our algorithm. As an illustration of the dynamics of the PWM and the rank of different candidate motif sets at different stages of our algorithm, a summary of the prediction process for the motif of the transcription factor CBF1 is given in Figure 2.
Figure 1. A flowchart of our algorithm.
Figure 2. An illustration of the dynamics of the PWM and the rank of five candidate motif sets at different stages of our algorithm in the motif prediction of the transcription factor CBF1.
Initial motif sets building
Our method begins by enumerating those patterns in STF that appear most often in the binding dataset BTF than in others. What we want to do first is to calculate the appearance probability of a pattern in STF, which is the probability that the pattern appears no less than n times in the binding dataset BTF. If a pattern b appears more often than other patterns in STF and its occurrence probability in a generic intergenic region is comparatively low, the calculated significance score of b would be relatively high. We will take patterns with the most highest significance scores as the candidates to build a number of initial candidate motif sets.
Binomial probability distribution model
The probability to observe exactly j occurrences of pattern b in the collection STF of segments extracted from the binding dataset BTF is estimated by the binomial distribution
where occ(b) is the occurrence times of pattern b in STF and f (b) is the probability that pattern b occurs in the intergenic region and is estimated as the relative frequency of pattern b in the union S = ∪TFSTF of segments extracted from all transcription factor binding datasets. The probability to observe n or more occurrences of the pattern b in STF is
We define the significance score sigTF (b) of a pattern b to TF as
The less probable pattern b in S appears more than n times in STF, the more probable will it be a binding site candidate for the transcription factor TF. We will take N distinct patterns with the most highest significance scores as the candidates.
For each of the N significant binding site candidates for the transcription factor TF, in view of evolution, collect all segments in STF whose patterns have no more than d Hamming distance matching to the candidate pattern to form an initial candidate motif set. Thus N initial candidate motif sets for the transcription factor TF are built at the end of this stage. As an example, the PWM and the rank of five initial candidate motif sets for the motif prediction of the transcription factor CBF1 are shown in Figure 2.
In this stage, a sampling method is used to expand and/or trim each of the N initial candidate motif sets M1, M2,.., MN. For our purpose, a false motif set MN+1 is created by randomly selecting e0 (e0 is equal to the maximum size of the N initial candidate motif sets) segments from the collection STF such that Mi ∩ MN+1 = ∅, for all i = 1, 2,..., N. In addition, let the collection S = ∪TFSTF of segments extracted from all transcription factor binding datasets represent the intergenic background and here be denoted as MBG.
Approximate maximum a posteriori (AMAP) measure
where ps, j is the frequency of nucleotide j at base position s in the candidate motif set Mi(which can be retrieved from the position specific scoring matrix (PSSM) of Mi), ni is the number of segments in Mi, and P(m|MBG) is the probability of the pattern of segment m in the motif set Mi under an expanded Bayesian network (EBN) model  developed from the background motif set MBG (EBN model will be discussed shortly).
The first part of the AMAP score is a negative entropy, which is higher if there are more similar patterns in the candidate motif set Mi. A motif set Mi with all identical patterns has the maximum negative entropy 0, whereas equal nucleotide frequencies at every position in the PSSM of Mi has the minimum negative entropy. And a segment m in the candidate motif set Mi which has a pattern much different from the background motif model built from MBG would have lower appearance probability P (m|MBG) and hence increases the score of the AMAP measure of Mi.
In each iteration, there are two steps for the sampler, the S-step and the M-step.
In the S-step, the sampler samples a site by randomly selecting a sequence from BTF and then randomly picking up a site in the selected sequence to extract a segment ms of length w. For 1 ≤ i ≤ N, if the sampled segment ms appears in Mi, segment ms will be removed from Mi if the AMAP score of the candidate motif set Mi increases after its removal; otherwise, segment ms will be kept in Mi. Note that the PSSM of the motif model Mi should be retrained if the sampled segment ms is removed from Mi.
where ni is the size of current motif set Mi, P(ni) equals , P(ms|Mi) and P(ms|MBG) are the probabilities of the content of the sampled segment ms under the PSSM model developed from the current motif set Mi and under an EBN model developed from the background MBG, respectively. The sampled segment ms will be considered to append into the motif set Mi with the highest appendant score . If is the highest score, then the sampled segment ms is appended into the false motif set MN+1 unless ms is already there and the current iteration stops here. If for some i, 1 ≤ i ≤ N, is the highest score, the sampled segment ms will be further checked in the M-step to see if we really want to append ms into Mi unless we have processed ms for Mi at the beginning of this S-step as in above and the current iteration stops here.
In the M-step, the sampler has to decide whether the newly sampled segment ms should be appended into the candidate motif set Mi or not. The AMAP measure again will be used to evaluate our decision. The sampled segment ms is appended into the candidate motif set Mi if and only if the score of the motif model Mi is increased once the sampled segment ms is appended to Mi. Note that the PSSM of the motif model Mi should be retrained after the sampled segment bs is appended to Mi. Now the M-step is done and the current iteration stops here.
The sampler will iteratively sample through the binding dataset BTF to expand and/or trim the N candidate motif sets M1, M2,.., MN so that their AMAP scores will keep increasing. The N candidate motif sets will tend to be invariant after a (larger) number of iterations. The stopping criterion of the sampling process is that all the N candidate motif sets are invariant in K consecutive iterations. The parameter K is usually set to be 1% of the size of STF.
Alternative sampling strategy
There is an alternative sampling strategy as follows.
In the S-step, the new sampler also randomly samples a site from a sequence in BTF to extract a segment msof length w. For 1 ≤ i ≤ N, if the pattern of the sampled segment ms appears in Mi, all the segments in Mi whose pattern is the same as that of ms will be removed if the AMAP score of the motif set Mi increases after their removal. Otherwise, these segments will be kept in Mi.
Also in the S-step, if is the highest among all , 1 ≤ i ≤ N + 1, then all segments in the set STF having the same pattern as that of the sampled segment ms will be appended into the false motif set MN+1 unless these segments are already there and the current iteration stops here. If is the highest for some i, 1 ≤ i ≤ N, the sampled segment ms will be further checked in the M-step to see if we really want to append those segments in the set STF having the same pattern as that of the sampled segment ms into Mi unless we have already processed those segments for Mi at the beginning of this S-step as in above and the current iteration stops here.
In the M-step, all the segments in the set STF having the same pattern as that of the sampled segment ms are decided to append to the candidate motif set Mi if and only if the AMAP score of Mi increases after these segaments are appended into Mi.
Method of dependency graphs and their expanded Bayesian networks
Considering the binding mechanism of transcription factors to specific DNA sites (motifs), there must be distinctive features for the specific motif regions from other intergenic regions which represent the background DNA sequence. Hence, it is conceivable that we can use a statistical model to capture the feature of a specific DAN site (motif) or a generic DNA intergenic region (background). Since the size of a candidate motif set Mi is often small, a PSSM model is commonly used for Mi instead of any other more sophisticated statistical model. However, the size of the background motif set MBG is usually large enough to be equipped with a more sophisticated one.
As reported in , a dependency graph model is used to fully capture the intrinsic interdependency between base positions in a motif or region. The establishment of dependency between two positions is based on a χ2-test from known sample data. An edge is established between two nodes (a node represents a base position) in the graph if the two corresponding base positions of the motif or region are dependent. After all dependent edges have being established completely, a dependency graph for the motif or region is constructed. An example of a dependency graph with 7 nodes is shown in Figure 3.
Figure 3. An example of dependency graph.
As reported in , although the dependency graph can fully capture the intrinsic interdependency between base positions in a motif or region, it is difficult, if not impossible, to perform statistical inference based on the dependency graph. To resolve the dilemma, the dependency graph is expanded to form a Bayesian network (which is a directed acyclic graph that facilitates statistical reasoning) by allowing a base position in the dependency graph to appear more than once in the Bayesian network as nominally distinct nodes. Figure 4 shows an example of an expanded Bayesian network of the dependency graph in Figure 3. For the detailed procedure of constructing an expanded Bayesian network (EBN) from a dependency graph, please see . In this paper, we use EBNs to model the background motif set MBG.
Figure 4. An example of expanded Bayesian network.
Continued with the same example of the motif prediction of the transcription factor CBF1, the PWM and the rank of the five candidate motif sets at the end of the iterative sampling stage are also shown in Figure 2, together with the final results at the end of the refinement stage.
Results and discussion
In order to search for the transcription factor binding sites that regulate gene expressions, we collected binding promotor sequences from the cDNA microarray hybridization (ChIp-chip array) of yeast genome . Each of the binding sequences may contain some unknown motifs that are implanted at unknown positions. These data represent the binding affinity of a target transcription factor to the promoter region of a gene in vivo. The experiment protocol assigns a binding p-value to each binding promoter sequence of the corresponding transcription factor. A sequence with binding p-value less than 0.001 is considered to be bound by the corresponding transcription factor. The threshold of 0.001 is set up to reduce the false positive identification in yeast genome-wide screening.
We collected the ChIp-chip array sequence data from the "Motif discovery results – Discovered motifs, version 24" at . For a transcription factor TF to be investigated, we collected all sequences with binding p-value less than the threshold 0.001 to TF into the binding dataset BTF. There are 65 binding datasets BTF being able to be collected from Harbison's website.
Accuracy measurement and comparison
To evaluate the performance of our program, we collected known specificities from many famous websites, such as YPD, SCPD, Transfac and from the literature with experimental evidence  to compare with the discovered specificities predicted by our program.
Among the 65 binding datasets BTF collected from Harbison's website, we chose 36 transcription factor binding datasets which have known specificities with experimental evidence to evaluate the performance of our program. The results of our program for the 36 transcription factor binding datasets are listed in Figure 5. It is deserved to be mentioned that the specificity reported for transcription factor PHO2 in Harbrison et al's website is "GTGCGsyGCG", while the predicted result of our program is "ATTATC". In this case, the newly found motif by our program is more consistent with the results reported by Barbaric et al  that PHO2 binds to an AT-rich region than the specificity reported in Harbrison et al's website.
Figure 5. Predicted results of our program compared with known evidence. The letter symbols used in the 'Known specificity' column have the following mapping: aA: a tT: t gG: g cC: c wW: at rR: ag mM: ac kK: tg yY: tc sS: gc dD: atg hH: atc vV: agc bB: tgc nN: atgc
In this study, we compared our program with two online programs, MDscan  and Cosmo . MDscan is a famous program that can be used to examine the ChIP-array selected sequences and search for DNA sequence motifs representing the protein-DNA interaction sites. It takes the advantage of combining two widely adopted motif search strategies, word enumeration and position-specific weight matrix updating, and incorporates the ChIP enrichment information to accelerate the search and enhance its success rate. The comparison of MDscan with our program is shown in Table 1. Also reported in Table 1 is the performance of our algorithm when the the PSSM model  instead of the EBN model  is used to model the background motif set MBG in the calculation of the AMAP scores of the N candidate motif sets and the appendant scores of the N + 1 motif sets. In Table 1, for each transcription factor, the number in each 'Rank' column indicates the rank of the predicted motif which is most consistent with the known evidence from the top ten predicted candidate motifs.
Table 1. Comparison of MDscan and our program.
As shown in Table 1, our approach with EBN background model outperforms the other two methods. Our approach with EBN background model gives 30 out of the 36 most predicted motifs for the corresponding 36 transcription factors with the 1st rank, while MDscan and our approach with PSSM background model give only 20 out of 36 and 15 out of 36 most predicted motifs with the 1st rank, respectively. Moreover, MDscan fails in discovering a motif for three transcription factor binding datasets, while our approach in this study is still able to predict a motif consistent with the known evidence.
Cosmo (constrained search for motifs) is a general purpose algorithm for conserved motif detection that allows the search to be supervised by specifying a set of constraints that the PWM of the unknown motif must satisfy. Such constraints may be formulated derived from prior biological knowledge about the structure of the transcription factor, such as the length of the motif intervals. Although Cosmo is based on the same two-component multinomial mixture model used in MEME, it employs the likelihood principle instead of the E-value criterion in MEME. In addition, three model types (OOPS, ZOOPS, or TCM) can data-adaptively be selected in Cosmo to achieve better performance. Since there is no prior knowledge used in our program, we compared it to the constraint-less version of the Cosmo program. On the other hand, since the Cosmo program reports only one motif PWM for a dataset, instead of a list of ranked candidate motif PWMs as in MDscan, we adopted only the rank 1 results of our program in this comparison. To evaluate the performance of both programs, we used the statistics proposed by Tompa et al. . For a (computational) tool at the site level, the performance statistics on a dataset are defined as follows:
where TP is the number of known sites overlapped by predicted sites, FN is the number of known sites not overlapped by predicted sites, and FP is the number of predicted sites not overlapped by known sites. To summarize the performance of a given tool over a collection of datasets, we compute the "combined" statistics as though C were one large dataset by adding TP, FP and FN respectively over the datasets in . Then the combined statistics of our program are Sn = 0.6698, PPV = 0.8206, and ASP = 0.7452, while those of Cosmo are Sn = 0.6573, PPV = 0.5134 and ASP = 0.5854. For the detailed Cosmo prediction results and the comparison of the two programs, please see Figure S1 (see Additional file 1). The comparison shows that our program can offer better performance than Cosmo, especially in the elimination of false positives.
The parameters used in our program include the sliding window size w used to extract segments from binding datasets BTF to form STF, the Hamming distance d used to collect segments from STF to establish initial candidate motif sets, the number of most dependent edges used to form a dependency graph for the background motif model and the number of parents used in the construction of an expanded Bayesian network from the dependency graph . The parameters used in our program to give the best predicted motifs for each of the 36 transcription factors are listed in Table S1 (see Additional file 2). Comparing the performance of the two sampling strategies discussed in the Method section, as shown in Figure S2 (see Additional file 3), we found that the alternative sampler is faster and has almost identical best predicted motifs with those by the primary sampler, except that transcription factors GCN4, HAP4 and PHO4 have the best predicted motifs one nucleotide position shift from those by the primary sampler. In addition, the alternative sampler is slightly better than the primary sampler in the sense that the best predicted motif for the transcription factor DIG1 promotes its rank from the 2nd place by the primary sampler to the 1st place by the alternative sampler.
Format: PDF Size: 30KB Download file
This file can be viewed with: Adobe Acrobat Reader
In this study, we employed the binomial probability model to establish a number of initial candidate motif sets, and used the method of dependence graphs and their expanded Bayesian networks to model the background motif set as a control to predict TFBSs (motifs) from a set of unaligned DNA sequences. The prediction results suggest that, overall, our algorithm outperforms MDscan since the predicted motifs are more consistent with previously known specificities reported in the literature and have better prediction ranks. And when compared with the constraint-less Cosmo program, our algorithm has a slightly higher combined sensitivity Sn, a much higher positive predictive value PPV and a higher average site performance ASP. However, the performance of our algorithm is not much better if the length of possible binding sites are too long (more than 12 bps). Further research is needed to discover long motifs.
Furthermore, variable spacing within binding sites is legitimate for some transcription factors while this study focuses on ungapped motif discovery. Programs such as BIPAD  and spaced dyad  have investigated into such a bipartitie sequence element discovery problem. Therefore another direction for our future research is to investigate into gapped motifs.
The authors declare that they have no competing interests.
CL and WY developed and implemented the method. All authors participated in discussions and writing of the paper.
We thank Yi-Sian Liao for helping the execution of our prediction program. We also thank the reviewers for their valuable suggestions which improve the presentation of this paper. This work was supported by the National Science Council, Taiwan, under Contracts NSC 93-3112-B-007-004 and NSC95-3114-P-002-005-Y.
This article has been published as part of BMC Bioinformatics Volume 9 Supplement 12, 2008: Asia Pacific Bioinformatics Network (APBioNet) Seventh International Conference on Bioinformatics (InCoB2008). The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/9?issue=S12.
Tompa M, Li N, Bailey TL, Church GM, Moor BD, Eskin E, Favorov AV, Frith MC, Fu Y, Kent WJ, Makeev VJ, Mironov AA, Noble WS, 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.
Comput Appl Biosci 1990, 6:81-92. PubMed Abstract
Bailey TL, Elkan C: The value of prior knowledge in discovering motifs with MEME. In Proceedings of the Third International Comference on Intelligent Systems for Molecular Biology. Menlo Park, CA: AAAI Press; 1995::21-29.
Tagle DA, Koop BF, Goodman M, Slightom JL, Hess DL, Jones RT: Embryonic ε and γ globin genes of a prosimian primate (Galago crassicaudatus). Nucleotide and amino acid sequences, developmental regulation and phylogenetic footprints.
Harbison CT, Gordon DB, Lee TI, Rinaldi NJ, Macisaac KD, Danford TW, Hannett NM, Tagne JB, Reynolds DB, Yoo J, Jennings EG, Zeitlinger J, Pokholok DK, Kellis M, Rolfe PA, KT T, Lander ES, Gifford DK, Fraenkel E, Young RA: Transcriptional regulatory code of a eukaryotic genome.
J Am Stat Assoc 1995, 90:1156-1170. Publisher Full Text
Motif discovery results – Discovered motifs, version 24 [http://fraenkel.mit.edu/Harbison/release_v24/final_set/Final_Motifs/] webcite