CisModule: De novo discovery of cisregulatory modules by hierarchical mixture modeling
See allHide authors and affiliations

Edited by Michael S. Waterman, University of Southern California, Los Angeles, CA (received for review April 23, 2004)
Abstract
The regulatory information for a eukaryotic gene is encoded in cisregulatory modules. The binding sites for a set of interacting transcription factors have the tendency to colocalize to the same modules. Current de novo motif discovery methods do not take advantage of this knowledge. We propose a hierarchical mixture approach to model the cisregulatory module structure. Based on the model, a new de novo motifmodule discovery algorithm, CisModule, is developed for the Bayesian inference of module locations and withinmodule motif sites. Dynamic programminglike recursions are developed to reduce the computational complexity from exponential to linear in sequence length. By using both simulated and real data sets, we demonstrate that CisModule is not only accurate in predicting modules but also more sensitive in detecting motif patterns and binding sites than standard motif discovery methods are.
Transcription factors (TFs) regulate genes by binding to their recognition sites. The common pattern of the binding sites for a TF is called a motif, usually modeled by a positionspecific weight matrix (PWM). Experimental methods such as DNase footprinting (1) and gelmobility shift assay (2, 3) have allowed the determination of some binding sites for selected TFs. Because these procedures are timeconsuming, several computational methods have been developed for de novo motif discovery, including progressive alignment (4, 5), the expectationmaximization algorithm (6, 7), the Gibbs sampler (8–12), word enumeration (13, 14), and the dictionary model (15, 16). The propagation model (17) and the recursive Gibbs motif sampler (18) have been developed for locating multiple motifs simultaneously. In addition, methods also exist that combine motif discovery with gene expression data (19–21) or phylogenetic footprinting (22, 23). These experimental and computational analyses have given us a good number of useful TF motifs. However, there are still many important TFs whose motifs remain to be characterized. What is more, molecular analyses have established that most eukaryotic genes are not controlled by a single site but by cisregulatory modules (CRMs), each consisting of multiple TFbinding sites (TFBSs) that act in combination (24–27). It can be argued that motif discovery is but an intermediate step toward the characterization of CRMs. Current approaches on module prediction such as those based on logistic regression (28, 29) or hidden Markov models (30, 31) depend on the availability of known motifs, i.e., PWMs for several TFs hypothesized to bind synergistically to regulatory modules. Clearly, we cannot apply these methods to the situations where no prior knowledge on the TFs is available, and in these cases we must resort to de novo motif discovery algorithms. We hypothesized that greater sensitivity and specificity can be achieved for motif discovery by considering the colocalization of different TFBSs and searched for modules and motifs simultaneously. It is clear that the task of module discovery and motif estimation is tightly coupled: on one hand, motif patterns and binding sites are essential for predicting regulatory modules; on the other hand, discovery of modules will greatly improve the performance of motif detection.
In this article, we propose a hierarchical mixture (HMx) model and develop a fully Bayesian approach for the simultaneous inference of modules, TFBSs, and motif patterns based on their joint posterior distribution. We test the approach by using both simulated and real data sets. Simulation studies show that, by capturing the combinatorial patterns of cooperating TFBSs, our algorithm detects modules accurately and is much more precise than standard motif discovery algorithms are in finding true binding sites. Similar improvement is observed when the method is tested on the known CRMs from a number of Drosophila developmental genes (26, 32, 33) and on the regulatory regions of a set of musclespecific genes (28). Our approach for de novo motifmodule discovery is of great current interest. Expression microarrays (34) and serial analysis of gene expression (35) have provided powerful means to identify clusters of genes tightly regulated during various cellular processes. Genes in the same clusters have a higher likelihood of sharing similar CRMs. Comparative analysis of multiple genomic sequences can further identify conserved regions enriched for such modules (36, 37). Finally, chromatin immunoprecipitation followed by microarray (ChIPonchip) is able to predict the binding locations of a TF in the whole genome with a resolution of 500–2,000 bp. These approaches are expected to provide sets of sequences enriched for CRMs involving an unknown or a partially unknown set of regulatory TFs. The identification of the CRMs within these sequences and the clarification of their structures, which are essential steps in understanding the regulatory networks, will depend on computational methods such as those proposed in this article.
Methods
HMx Model for CisRegulatory Modules. Our goal was to search for the binding sites for K different TFs within the CRMs of a given set of sequences S. We proposed a twolevel HMx model for CRMs. At the first level, the sequences can be viewed as a mixture of CRMs, each of length l, and pure background sequences outside the modules; at the second level, modules are modeled as a mixture of motifs and withinmodule background. Detailed specification of the HMx model is illustrated in Fig. 1. The background sequences, both the regions outside the modules and the nonsite segments within the modules, are modeled by a firstorder Markov chain θ_{0}.
It is helpful to think of the HMx model as a stochastic machinery that generates sequences. Suppose the width of the kth motif is w_{k} and its product multinomial model (PWM) is Θ _{k} (k = 1,..., K). Starting from the first sequence position, we made a series of random decisions of whether to initiate a module or generate a letter from the background model, with probabilities r and 1  r, respectively. If a module was started at position i, within the region of [i, i + l  1], we generated background letters or initiated the kth motif sites, with probabilities q _{0} and q_{k} (), respectively. If a site for the kth motif was initiated at position n, we generated w_{k} letters from its PWM Θ _{k} and placed them at [n, n + w_{k}  1]. After we reached the end of the current module at position i + l  1, the decision at the next position was reverted back to the choice between sampling from the background or initiating a new module. Let M denote the module indicators and A _{k} denote the indicators for the binding sites for the kth motif. We used S(M) to denote the CRMs and S(M^{c} ) to denote the background outside the modules. To simplify the notation, we let A = {A _{0}, A _{1},..., A _{K} }, where A _{0} indicates the nonsite background sequences in the modules, Θ = {θ_{0}, Θ _{1},..., Θ _{K} }, q = {q _{0}, q _{1},..., q_{K} }, and W = {w _{1},..., w_{K} }. The notations for the model are summarized in Table 1.
Under the HMx model, the complete sequence likelihood with M and A given is Combining Eq. 1 with the prior distributions for all the parameters gives rise to the joint posterior distribution: where conjugate prior distributions are prescribed, i.e., a product Dirichlet distribution with parameter β _{k} (a w_{k} × 4 matrix) for (Θ _{k}w_{k} ), a Dirichlet distribution with parameter α (a vector of length K + 1) for q, and Beta(a, b) for r. We put a Poisson(w _{0}) prior on w_{k} (k = 1,..., K).
Bayesian Inference. We regarded M and A as missing data and used the Gibbs sampler (38–40) to perform Bayesian inference. Gibbs sampling algorithms are widely used for motif finding (8, 9, 17), but our problem was much more complex than traditional motif discovery because of its hierarchical structure. With a random initiation, our algorithm (CisModule) iteratively cycles through the steps of parameter update and modulemotif detection (Fig. 2A ). (i) Given current modules and motif sites (M and A), we updated all the parameters Ψ = (Θ, q, W, r) by sampling from their conditional posterior distributions [ΨM, A, S] (see Appendix A). (ii) Given current values of the parameters, we sampled modules and motif sites from the conditional distribution [M, AΨ, S]. Without loss of generality, suppose the sequence data are S = {x _{1} x _{2},..., x_{L} } = x _{[1,L]}. The computational bottleneck is the step of modulemotif detection. Sampling modules and sites naively results in a computational complexity of O((K^{l}l) ^{L/l} ), which increases exponentially with the total sequence length L. By using stochastic recursions we reduced the complexity to O(KL). First, we performed “forward summation” to compute P(SΨ) using the recursion (Eq. 5 in Appendix B). Then “backward sampling” was used to generate the module indicators as follows. Starting from n = L, at position n, we decided whether (i) x_{n} was at the last position of a module or (ii) x_{n} was from the background. The probabilities of these two events are proportional to the terms A_{n} (Ψ) and B_{n} (Ψ) in Eq. 6 in Appendix B, which are already computed from the forward summation. Depending on choosing event i or event ii, we moved to position n  l or n  1 and repeated the binary decision process. In this way, we generated all the module indicators. Once modules were updated, we again used forward summation (see Eq. 7 in Appendix B) and backward sampling to update motif indicators within each module. Suppose we have sampled the motif indicators backward up to position m in the current module. The sequence segment x _{[} _{m} _{} _{wk+1,m} _{]} (k = 0,..., K) is drawn as a background letter (k = 0, w _{0} = 1) or a site for one of the K motifs with probability proportional to the K + 1 terms in Eq. 7 . Apparently, because sites are sampled for each module separately, the combinatorial site patterns in the individual modules can be different.
By using the samples from the joint posterior distribution (Eq. 2 ), we obtained marginal distributions of the width and number of sites for each motif by smoothing their sampling histograms by means of a moving average. Based on the marginal modes that can be found through enumeration, we estimated ŵ_{k} and n̂_{k} (k = 1,..., K). The top n̂_{k} ŵ_{k} mers that were most frequently sampled as sites for the kth motif were aligned as output sites. Furthermore, we inferred the modules by the marginal posterior probability of each sequence position being sampled as within modules. The positions where this probability is >0.5 were output as modules (Fig. 2B ).
Strategies on l and K. In the discussion above, module length l and TF number K were left as userinput parameters. We now discuss how to determine l and K in case we have no prior knowledge of them.
An extra conditional sampling by a Metropolis update can be performed to determine the most likely module length. Let l be the current module length. We propose a new one, l + δ (δ = ±10), and accept it with the Metropolis ratio, where the prior distribution π(l) is geometric with mean l _{0} (usually between 100 and 200).
It is often desirable to provide some information about the TF number K. This can be formulated as a Bayesian model selection problem. Let H_{K} (K = 1, 2,...) denote the hypothesis that there are K motifs (TFs) and H _{0} denote the null hypothesis that S is generated from pure background. With π(H_{K} ) ∝ (1/3) ^{K} as the prior, we calculate the posterior odds of H_{K} over H _{0}, where P (SH _{0}) is of known form and P (SH_{K} ) can be calculated by importance sampling (see Appendix C for details). Thus we can run CisModule with K = 1,..., K_{m} , where with K_{m} the algorithm stops detecting new motifs, and treat the K ^{*} ∈ {1,..., K_{m}  1} that maximizes the posterior odds (Eq. 4 ) as our estimated number of motif types.
Results
We tested CisModule on both simulated and real biological data sets. Data Sets 1–4 are published as supporting information on the PNAS web site.
Simulation Studies. It is known that E2F, YY1, and c_MYC are potential cooperating factors (41). Thus, in our simulation, motif sites were generated according to the weight matrices of these three TFs based on TRANSFAC (42) matrix accession numbers, V$E2F_03, V$YY1_02, and V$MYCMAX_02, respectively. The background sequences were generated by a firstorder Markov chain with parameters estimated by >2,000 upstream 1kb sequences from the ensembl genome database (www.ensembl.org). In the first simulation study, each module was 100 bp long and contained one E2F site, one YY1 site, and one c_MYC site, randomly placed in the module. One data set consisted of 40 sequences, each 500 bp in length, and 20 modules were randomly located in these sequences. In the second simulation study, each data set contained 30 sequences, each 800 bp in length. Twenty 200bplong modules of different site combinations were generated, where four of them contained only three E2F sites, eight of them contained one E2F site, two YY1 sites, and one c_MYC site, and the rest contained one E2F site, one YY1 site, and two c_MYC sites. This different site combination mimics the fact that one TF (E2F) may work with different partners. For each of the simulation studies above, 10 data sets were generated independently. We applied CisModule to these data sets and fixed the module length to be 100 and 200 bp, respectively. The number of motifs K was set as 3 in both studies.
We evaluated our prediction for modules by their total length and coverage of true sites. The total lengths of our predicted modules were 2,009 and 4,108 bp on average for the two simulation studies, corresponding to excess rates of 0.5% and 2.7% over the actual module lengths (2,000 and 4,000 bp), respectively. The average true site coverage rates of the predicted modules were 84.3% and 94.0%, which showed that our module prediction was very informative with a high coverage of true sites and a low excess in length. In terms of motif discovery, we compared our predictions with MEME (7) and BioProspector (BP) (11) on these data sets. We set these algorithms to run multiple times and output the top 20 motifs they found. From Table 2 we see that, for all of the cases, CisModule showed the greatest success rates of discovering the correct motif patterns and found more true sites with comparable numbers of false positives. The improvement over MEME and BP was especially significant for weakly conserved motifs (c_MYC). These results demonstrate that the HMx model captures the colocalization of TFBSs and CisModule is capable of using this information to improve de novo motif discovery.
We repeated the experiments with K = 4, and, for all of the data sets, CisModule did not predict any new (false) motifs. By using the posterior odds calculation, CisModule correctly estimated the true motif numbers (K ^{*} = 3) for 19 of the 20 data sets. We also tested our algorithm assuming l unknown. The most likely module lengths predicted by CisModule were within 30 bp of the true lengths for 18 data sets.
Homotypic Regulatory Modules in Drosophila. Analyses of experimental data from the early developmental Drosophila gene enhancers show that these regions are highly enriched of homotypic clusters, i.e., multiple binding sites for one TF are tightly clustered together (32, 33). More than 60 regulatory modules for 20 different genes were collected and the known regulatory interactions using published data were annotated (32). We built three sequence sets, each of which contained all the CRMs for one of the three most frequent binding motifs in their data sets, Bicoid (Bcd), Hunchback (Hb), and Krüppel (Kr). Thirtyfour experimentally reported sites are in our data sets: 12 Bcd sites in three sequences, 14 Hb sites in four sequences, and 8 Kr sites in two sequences. Because binding sites are not reported in the remaining sequences, we scanned the data sets for putative target sites based on the known PWMs for the three TFs (32). These scannedbased sites served as an alternative basis for our comparison.
We applied CisModule to the three data sets with K = 1 (because the modules are clusters of binding sites for one TF) and l = 100. By the modulesampling step, CisModule provides more information through the marginal posterior probability of each position being sampled as within modules (P_{m} in Fig. 2B ). Some examples of the predicted modules with this probability are illustrated in Fig. 3A . For each data set, we selected all the sequence positions with this probability greater than a given value x, denoted by S(x), and calculated the density of S(x), defined as the ratio of the number of highscore sites (those within the top 0.5% in scanning) to the size of S(x), for x varying from 0.1 to 0.8 (Fig. 3B ). When x was increased to 0.9, the sizes of S(x) were too small to calculate the densities. From the figure it is clear that for x ≤ 0.5 the densities increase with x, i.e., those sequence positions that are more likely to be sampled within modules have a higher density of top sites. The densities for x = 0.5 (corresponding to the broken horizontal lines in Fig. 3A ) for all of the three data sets were significantly higher than 0.5% with P values < 3E6. If we further increased x (≥0.6), all of the positions in S(x) were selected from module regions, and thus the densities were approximately the same for different x.
As a comparison, we also applied MEME and BP to the data sets to find the top 20 motifs. From Table 3 we see that CisModule not only successfully discovered correct motifs in all three data sets but also found many more experimentally reported sites than the other two methods did. In total it reached a sensitivity of 56% for these reported sites. The numbers of output sites by CisModule were slightly more than those of scannedbased sites, because some weakly conserved sites missed by scanning can be detected by CisModule if they are close enough to other sites. The logo plots (43) for the three motifs found by CisModule are shown in Fig. 4, which is published in supporting information on the PNAS web site, where we see that they are consistent with the known consensus sequences listed in the figure legend. Furthermore, with the jaspar database (44, 45), the known Hb motif ranked number 1 compared with our predicted Hb matrix with a similarity score of 97/100. (The known motifs for Bcd and Kr are not collected in the jaspar database, so we did not compare these two factors to the database.) We also repeated the experiments with K = 2. For the Bcd and Kr data sets, CisModule did not output any new motifs. For the Hb data set, a weak motif with consensus GCMGGNM showed cooccurrence, but the posterior model odds was maximized at K ^{*} = 1. These results agreed with the homotypic cluster phenomenon.
MuscleSpecific Regulatory Regions. Logistic regression was proposed as a predictive model for the regulatory regions for musclespecific expression (28), where five TFs (Mef2, Myf, Sp1, SRF, and TEF) known to control the expression were used as predictors. The positive training set for the logistic regression was composed of 29 regulatory sequences sufficient for skeletalmusclespecific expression that have been experimentally localized to within 200 bp. We annotated 25 experimentally reported binding sites, 10 for Mef2, 7 for TEF, and 8 for SRF. Besides, by using the weight matrices for the five TFs (figure 1A in ref. 28), we scanned the 29 sequences and detected 19, 12, 23, 13, and 20 putative sites for the five TFs above at a falsepositive error rate of 5E4, which provided estimates for the numbers of target sites. Two data sets were constructed by adding 10 and 40 upstream sequences (200 bp each) randomly extracted from the ensembl database to the 29 positive training sequences. We tested how resistant the algorithm was to the presence of noisy sequences (those random upstreams). CisModule was applied to these data sets with K = 5 and l = 150. We also applied MEME and BP to the same data sets to output the top 20 motifs they could find. The logo plots for the motifs found by CisModule are shown in Fig. 5, which is published as supporting information on the PNAS web site.
It turns out that all three algorithms successfully found the Sp1 motif (GC box). We focus our comparison on the other four factors. The results are summarized in Table 4, where we tabulate among all the predicted sites from each method the number of reported sites (n _{1}), the number of putative sites in positive sequences that do not overlap with reported sites (n _{2}), and the number of falsepositive sites in random sequences (n _{3}). The nature of putative sites (n _{2}) is ambiguous because they may be unreported binding sites or false positives. For Mef2 and TEF, CisModule found more reported sites and usually fewer falsepositive sites for different cases. Furthermore, CisModule was the only algorithm that discovered the SRF motif (with a phase shift of two bases). None of the methods found the motif for Myf. From the summary in Table 4 we see that the sensitivity of CisModule in discovering reported sites (n _{1}) is 88% (22 of 25) and 68% (17 of 25) for the data sets with 10 and 40 random sequences, respectively, which is much higher than the sensitivity of the other two methods. CisModule is also most resistant to the mixed random sequences with the fewest falsepositive predictions (n _{3}). These results confirm the notion that module sampling based on the combinatorial effects of several motifs is more stable than sampling each motif individually. Taking the data set with 40 random sequences as an example, we found that 54% of our predicted modules were from the 29 positive sequences, but only 34% of the output sites predicted by MEME were from the positive sets. The predicted modules that do not overlap with positive sequences are most likely false positives, but the possibility exists that some might be unreported modules.
Discussion
The HMx model assumes that TFBSs are located within some relatively short sequence segments, the CRMs. The benefit of this model is that it captures the spatial correlation between different binding sites. It is clear that the more tightly clustered the motif sites, the more information the HMx model gains. Based on the model, a Bayesian module sampler, CisModule, is developed to simultaneously infer the motif modules and the binding sites for a set of TFs by means of the Gibbs sampling approach. The module detection step utilizes the combination of several motifs, which significantly enhances the sensitivity of the method.
As is true for all de novo motif discovery algorithms, CisModule may sometimes be trapped in local modes. To reduce this possibility, multiple trials are often needed. If some prior information is available for a particular data set, we can use it to initiate CisModule. For example, if we know that the sequences are controlled by one TF, and we are interested in finding the binding sites for this TF and its cooperating TFs, the weight matrix for the known TF can be used to prescribe more specific prior distributions. This will lead to faster convergence to the correct motif patterns.
An interesting future work would be to incorporate the information from comparative genomics into CisModule. Greater prior probabilities for modules and sites can be assigned to the regions that are highly conserved across species of appropriate evolutionary distances. This will effectively reduce the falsepositive discovery and is especially important for higher organisms, whose upstream sequences are long and regulatory mechanisms are complex. Finally, the model presented here should be regarded as a first step to the development of realistic models for de novo motifmodule discovery. The HMx model captures the colocalization tendency of cooperating TFBSs but not their order or precise spacing. It is possible that additional refinements to the model may further enhance its utility.
Acknowledgments
This work was supported by a National Institute of General Medical Sciences grant (to W.H.W.).
Footnotes

↵ ‡ To whom correspondence should be addressed. Email: wwong{at}stat.harvard.edu.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: TF, transcription factor; TFBS, TFbinding site; CRM, cisregulatory module; HMx, hierarchical mixture; PWM, positionspecific weight matrix; BP, BioProspector; Bcd, Bicoid; Hb, Hunchback; Kr, Krüppel.
 Copyright © 2004, The National Academy of Sciences
References

↵
Galas, D. J. & Schmitz, A. (1978) Nucleic Acids Res. 5, 31573170. pmid:212715

↵
Fried, M. & Crothers, D. M. (1981) Nucleic Acids Res. 9, 65056525. pmid:6275366

↵
Garner, M. M. & Revzin, A. (1981) Nucleic Acids Res. 9, 30473060. pmid:6269071

↵
Stormo, G. D. & Hartzell, G. W. (1989) Proc. Natl. Acad. Sci. USA 86, 11831187. pmid:2919167

↵
Hertz, G. Z. & Stormo, G. D. (1999) Bioinformatics 15, 563577. pmid:10487864
 ↵
 ↵

↵
Lawrence, C. E., Altschul, S. F., Boguski, M. S., Liu, J. S., Neuwald, A. N. & Wootton, J. (1993) Science 262, 208214. pmid:8211139
 ↵

↵
Liu, X., Brutlag, D. L. & Liu, J. S. (2001) Pac. Symp. Biocomput. 6, 127138.

↵
Zhou, Q. & Liu, J. S. (2004) Bioinformatics 20, 909916. pmid:14751969

↵
Sinha, S. & Tompa, M. (2002) Nucleic Acids Res. 30, 55495560. pmid:12490723

↵
Hampson, S., Kibler, D. & Baldi, P. (2002) Bioinformatics 18, 513528. pmid:12016049

↵
Bussemaker, H. J., Li, H. & Siggia, E. D. (2000) Proc. Natl. Acad. Sci. USA 97, 1009610100. pmid:10944202
 ↵
 ↵

↵
Thompson, W., Rouchka, E. C. & Lawrence, C. E. (2003) Nucleic Acids Res. 31, 35803585. pmid:12824370
 ↵
 ↵

↵
Conlon, E. M., Liu, X. S., Lieb, J. D. & Liu, J. S. (2003) Proc. Natl. Acad. Sci. USA 100, 33393344. pmid:12626739

↵
Wang, T. & Stormo, G. D. (2003) Bioinformatics 19, 23692380. pmid:14668220
 ↵

↵
Yuh, C. H., Bolouri, H. & Davidson, E. H. (1998) Science 279, 18961902. pmid:9506933

Loots, G. G., Locksley, R. M., Blankespoor, C. M., Wang, Z. E., Miller, W., Rubin, E. M. & Franzer, K.A. (2000) Science 288, 136140. pmid:10753117

↵
Berman, B. P., Nibu, Y., Pfeiffer, B. D., Tomancak, P., Celniker, S. E., Levine, M., Rubin, G. M. & Eisen, M. B. (2002) Proc. Natl. Acad. Sci. USA 99, 757762. pmid:11805330

↵
Banerjee, N. & Zhang, M. Q. (2003) Nucleic Acids Res. 31, 70247031. pmid:14627835
 ↵

↵
Krivan, W. & Wasserman, W. W. (2001) Genome Res. 11, 15591566. pmid:11544200

↵
Frith, M. C., Hansen, U. & Weng, Z. (2001) Bioinformatics 17, 878889. pmid:11673232

↵
Sinha, S., van Nimwegan, E. & Siggia, E. D. (2003) Proc. Int. Conf. Intell. Syst. Mol. Biol. 11, 292301.

↵
Lifanov A. P., Makeev, V. J., Nazinna, A. G. & Papasenko, D. A. (2003) Genome Res. 13, 579588. pmid:12670999

↵
Makeev, V. J., Lifanov A. P., Nazinna, A. G. & Papasenko, D. A. (2003) Nucleic Acids Res. 31, 60166026. pmid:14530449

↵
Schena, M., Shalon, D., Davis, R. W. & Brown, P. O. (1995) Science 270, 467470. pmid:7569999

↵
Velculescu, V. E., Zhang, L., Vogelstein, B. & Kinzler, K. W. (1995) Science 270, 484487. pmid:7570003
 ↵

↵
Loots, G. G., Ovcharenko, I., Pachter, L., Dubchak, I. & Rubin, E. M. (2002) Genome Res. 12, 832839. pmid:11997350
 ↵
 ↵

↵
van Ginkel, P. R., Hsiao, K. M. & Farnham, P. J. (1997) J. Biol. Chem. 272, 1836718374. pmid:9218478

↵
Wingender, E., Chen, X., Hehl, R., Karas, H., Liebich, I., Matys, V., Meinhardt, T., Pruss, M., Reuter, I. & Schacherer, F. (2000) Nucleic Acids Res. 28, 316319. pmid:10592259

↵
Schneider, T. D. & Stephens, R. M. (1990) Nucleic Acids Res. 18, 60976100. pmid:2172928

↵
Sandelin, A., Alkema, W., Engström, P., Wasserman, W. & Lenhard, B. (2004) Nucleic Acids Res. 32, D91D94. pmid:14681366

↵
Lenhard, B. & Wasserman, W. (2002) Bioinformatics 18, 11351136. pmid:12176838