DREME: motif discovery in transcription factor ChIP-seq data

2.2.1 Running motif discovery algorithms

We run DREME using its default settings. When only one set of sequences is input, DREME creates a dinucleotide-shuffled version of the input sequence set as the negative sequence set. For discriminative motif discovery, we provide two sets of input sequences to DREME. In both cases, the significance threshold for motifs is E = 0.05.

We compare the speed and accuracy of several popular motif discovery algorithms with that of DREME using the 13 mESC ChIP-seq datasets. Two of these algorithms [MEME (Bailey and Elkan, 1995) and nestedMICA (Down and Hubbard, 2005)] are extremely slow when run on more than 500 sequences, so we randomly choose 500 sequences from each ChIP-seq dataset to create datasets of reduced size. We run the other three algorithms [WEEDER (Pavesi et al., 2004), Trawler (Ettwiller et al., 2007) and Amadeus (Linhart et al., 2008)] on the same datasets that were used above with DREME. Where possible, we parameterize the programs to search for motifs of 4–7 bp. The exceptions are WEEDER up to 8 bp; Trawler minimum 4 bp; Amadeus exactly 7 bp. The number of motifs to search for is a required parameter for nestedMICA, and we set it to 10 to keep running times reasonably short. (It does not finish running after 134 h on the full-size E2f1 dataset, which contains 20 699 sequences.) We set the maximum number of motifs to find to 20 for MEME; WEEDER, Trawler and Amadeus have no such parameter. With MEME, we use a 0-order background model and with nestedMICA we use a 0-order four class background model. As background sequences for DREME, Amadeus and Trawler, we use 1, 5 and 10 dinucleotide-shuffled copies, respectively, of the foreground sequences. With WEEDER, we use its pre-computed background model for mouse, which contains the frequencies of all 8mers in regions 1000 bp upstream of genes. Exact command lines and more details are given in the Supplementary Material.

2.2.2 ChIP-seq datasets

We use 13 mESC ChIP-seq datasets (Chen et al., 2008), 3 mouse erythrocyte datasets (Cheng et al., 2009; Kassouf et al., 2010; Tallack et al., 2010) and 1 human lymphoblastoid cell line dataset (Ramagopalan et al., 2010). The mESC datasets are for 12 TFs that are key to the maintenance of pluripotency, plus CTCF. The mouse erythrocyte datasets are for Gata1, Klf1 and Tal1, key regulators of erythropoiesis. The human lymphoblastoid cell line data is for the vitamin D receptor (VDR) and includes separate ChIP-seq data for cells before and after stimulation with calcitriol. To prepare the datasets for use with motif discovery algorithms, we map the (centers of the) ChIP-seq peaks declared by the respective authors to the genome, and extract the 100 bp of genomic sequence centered on each peak. [The one exception is Gata1, where we used the peaks declared by Tallack et al. (2010) in the Cheng et al. (2009) data.]

2.2.3 Identifying discovered motifs

To determine the success of motif discovery, we compare the motifs to a database of known TF motifs using TOMTOM (Gupta et al., 2007). The motif database comprises all vertebrate motifs in the JASPAR database (Sandelin et al., 2004) (146 motifs) and all mouse motifs in the UniProbe database (Berger and Bulyk, 2009) (386 motifs). Because the two motif reference datasets lack a motif for the key pluripotency TF Nanog, we also include an in vitro motif for Nanog taken from the literature (Jauch et al., 2008), for a total of 533 motifs. We only consider statistically significant TOMTOM predictions (E ≤ 0.05, corrected for 533 tests).

Because many TFs have very similar DNA-binding affinities, many of the known motifs in this database are similar to each other, so some discovered motifs have significant matches to multiple motifs in the database. Therefore, in order to determine to which TF a discovered motif may correspond, we use prior knowledge from the literature and evidence of increased expression in the ChIP-ed tissue type. To resolve multiple matches, we take the most significant TOMTOM prediction among matching TFs that are either upregulated in, or known to be central regulators of, the ChIP-ed cell type. We consider a TF as upregulated in mouse ES cells if it was so marked for at least one experiment in cell type EF0_000462 in the Gene Expression Atlas (Kapushesky et al., 2010) (http://www.ebi.ac.uk/gxa/). For motifs discovered in mouse erythrocytes and human lymphoblastoid cell lines, we restrict matching candidate motifs to TFs previously shown to be important to transcriptional regulation in those cell types.

3.1.1 Discovering motifs in mouse ES cell ChIP-seq data

DREME reports between 10 and 33 significant motifs (P < 0.05) for each of the mESC datasets (Table 1, column 3), and the number of significant motifs it finds is highly correlated with the number of ChIP-seq peaks in the input dataset. It discovers the ChIP-ed (primary) TF motif in 10 of the 13 mESC ChIP-seq datasets (Table 1, column 4), and it is the most significant motif found by DREME in those 10 cases. In two of the remaining three cases (Nanog and E2f1 ChIP-seq datasets), DREME finds a highly significant motif that is similar to the known motif, but the TOMTOM p-value does not meet our 0.05 threshold due to the large number of multiple tests (533 reference motifs). However, the uncorrected p-values are quite low (Nanog, p = 0.0021; E2f1, p = 0.0016). Since DREME reports 24 motifs in the Nanog dataset, we could reject a slightly different null hypothesis that stated that the known Nanog motif was not found (Nanog, p = 0.05, corrected for 24 tests) and similarly for E2f1 (E2f1, p = 0.04, corrected for 25 tests). In the final case (Smad1 ChIP-seq dataset), DREME does not find a motif that matches the only available in vitro Smad-family motif (UniProbe Smad3_primary), nor the motif reported by Chen et al. (The complete DREME and TOMTOM output on the mESC datasets is given in the Supplementary Material.)

In addition to finding the primary motif, DREME discovers binding motifs for up to 12 cofactors in each of the mESC ChIP-seq datasets (Table 1, column 4). Many of the motifs found by DREME in each mES cell dataset correspond to 1 of 12 key pluripotency TFs (shown in bold font in Table 1), and these predictions correspond well to the actual overlap of ChIP-seq peak regions reported by 0 (). For example, in addition to finding the motif for the ChIP-ed factor in the Tcfcp2l1 dataset, DREME also discovers the motifs for pluripotency factors STAT3, Klf4, Esrrb, Sox2, Oct4 and c-Myc/n-Myc (Myc-family motifs are essentially indistinguishable). In the CTCF dataset, DREME identifies the CTCF motif as well an unknown motif that was recently found to frequently occur upstream of the CTCF motif (consensus GCTGCAGT) (Boyle et al., 2010) in human cells.

Many of the other motifs found by DREME can be assigned to TFs that are known to be upregulated in mES cells, but whose roles in the maintenance of pluripotency are unknown (shown in normal font in Table 1, column 4). Several of these cofactor TF motifs are found in more than one mESC ChIP-seq dataset. These cofactor motifs include Sp1, CREB/ATF and Zic3, which was recently suggested as being required for the maintenance of pluripotency (Lim et al., 2010). On the other hand, we are currently unable to assign many of the motifs discovered by DREME to any known motif. Although some of these unassigned motifs may be artifacts, others may simply correspond to TFs whose motifs are not contained in the two motif databases we searched.

Several of the motifs discovered by DREME in the mES cell datasets differ substantially from those reported by Chen et al. in their original data analysis (0, ). Of particular note, DREME discovers the monomer motifs for Oct4 and Sox2 in their respective ChIP-seq datasets, rather than the Oct-Sox heterodimer motif. The in vivo binding motifs discovered by DREME for Oct4 and Sox2 closely match the in vitro motifs derived using protein binding microarray (PBM) technology (Berger and Bulyk, 2009), suggesting that the PBM motifs accurately capture DNA binding of these factors in vivo (Fig. 1a and b). This illustrates a particular advantage of focusing on short, core motifs and we discuss this issue in detail below. We also explore the hypothesis (suggested by Chen et al.) that Oct4 and Sox2 bind DNA exclusively as a homodimer. (The fact that the most significant motif found by DREME in each of the two ChIP-seq datasets is the primary factor's suggests this hypothesis may be false.)

DREME finds a motif similar to the known E2f1 motif in the E2f1 ChIP-seq dataset. This is in contrast to Chen et al., who report not finding a motif for E2f1 using either the WEEDER (Pavesi et al., 2004) or nestedMICA (Down and Hubbard, 2005) motif discovery algorithms. A previous study of E2f1 ChIP-chip data also failed to find a convincing E2f1 motif, concluding that E2f1 is mainly recruited to promoters via a mechanism distinct from recognition and binding to the consensus site (Bieda et al., 2006). Nonetheless, in the E2f1 dataset, the second most significant motif discovered by DREME is the word ATGGCG, which occurs in 1274 of the 20 699 ChIP-seq peak regions (E = 10 ^{− 98}) and resembles the in vitro motif for E2f-family member E2f3 (Fig. 1d). The strong enrichment of this word in the E2f1 ChIP-seq data and its resemblance to the in vitro E2f3 motif suggests that it may represent at least a subset of the in vivo binding sites of E2f1 in mES cells. The DREME motif also is similar to the motif for YY1 (JASPAR MA0095.1, consensus ATGG), suggesting the hypothesis that the sites discovered by DREME might bind both (but not simultaneously) E2f1 or YY1 in mES cells.

3.1.2 Discovering motifs in mouse erythrocyte ChIP-seq data

When we apply DREME individually to three ChIP-seq datasets from mouse erythroid cells [TFs Gata1 (Cheng et al., 2009), Klf1 (Tallack et al., 2010) and Tal1 (Kassouf et al., 2010)], the top motif is the ChIP-ed motif in the first two cases, but not in the last case. With the Gata1 dataset, DREME finds the Gata1 motif (E = 10 ^{− 892}) followed by the motif of SPI1 (E = 10 ^{− 79}, a key myeloid developmental regulator (Zakrzewska et al., 2010). It also finds the Klf1 motif (E = 10 ^{− 37}) and and E-box motif (E = 10 ^{− 19}) that might be due to Tal1 binding. DREME also discovers the motif for Runx1 (E = 10 ^{− 2}), which may also be involved in erythropoeisis (Yokomizo et al., 2008). In all, DREME finds 21 significant motifs in the Gata1 dataset.

DREME finds a total of six significant motifs in the Klf1 dataset. It identifies the Klf1 motif as the most significant (E = 10 ^{− 49}), followed by Gata1 (E = 10 ^{− 42}). Interestingly, the third most significant motif is the ‘secondary’ Klf motif discovered from in vitro binding data (Berger and Bulyk, 2009) (E = 10 ^{− 11}).

In the Tal1 ChIP-seq dataset, the most significant motif found by DREME (among 10 motifs found) is Gata1 (E = 10 ^{− 262}), whereas the actual Tal1 binding motif (E-box) is the third most significantly found motif (E = 10 ^{− 26}). The lower enrichment of the Tal1 motif compared with Gata1 may be due to the fact that Tal1 often binds DNA as a complex with Gata, Lmo2, Ldb1 and E47 (Wadman et al., 1997). In this complex, Tal1 has reduced contact with the DNA, binding to a motif consisting of one-half of an E-box followed ~10 bp later by a Gata1 site (Kassouf et al., 2010). This ‘semi-direct’ binding would reduce the prevalence of the full Tal1 motif in the Tal1 ChIP-seq dataset. DREME also finds motifs for Klf1 (E = 10 ^{− 23}) and SPI1 (E = 10 ^{− 16}) in this dataset.

3.1.3 Human lymphoblastoid cell line motifs

We apply DREME to two vitamin D receptor (VDR) ChIP-seq datasets from lymphoblastoid cells before and after stimulation with calcitriol (Ramagopalan et al., 2010). DREME finds 19 significant motifs in the unstimulated cell data, and only 11 in the stimulated cell data (E < 0.05). VDR binds DNA as a heterodimer with retinoid X receptor (RXR) and the left and right halves of the binding motif both have the consensus sequence (A/G)(A/G)G(T/G)TCA. DREME finds a single motif matching the consensus in the ChIP-seq data from the calcitriol-stimulated cells (fourth most significant motif, E = 10 ^{− 11}). The top-scoring motif in both datasets is an ETS-factor motif, possibly due to binding of Ets-1, which has been shown to cause VDR to function as a constitutive activator (Tolón et al., 2000). The second most significant motif in the stimulated cell data is an E-box, which has been proposed as being involved in negative regulation by VDR (Kato et al., 2007). The third most significant motif found in stimulated cell data strongly resembles STAT1, which is known to interact with VDR (Vidal et al., 2002). The fifth motif found is another Ets motif, and the sixth motif is a Klf motif. The seventh and eighth motifs significantly resemble Runx2 and Ap1/Fos/Jundm2 (TOMTOM E < 0.05), respectively, both of which are known to interact with VDR (Marcellini et al., 2010). In all, six of the eight most significant motifs found in the stimulated cell data appear to represent the ChIP-ed factor (VDR) or likely cofactors.

3.2.1 Determining if binding is solely via a heterodimer

It has been suggested that Sox2 and Oct4 bind their targets exclusively as a heterodimer in mES cells (0, ). If this is so, we should not observe enrichment of either motif when we use the ChIP-seq datasets for these two TFs as the positive and negative inputs to DREME. We performed this experiment (using default values for all DREME parameters), and found that the Oct4 dataset is enriched in Oct4 binding sites compared with the Sox2 dataset, and vice versa.

Using Oct4 as the positive and Sox2 as the negative input, DREME discovers the Oct4 motif (Fig. 2a, E = 10 ^{− 29}), but not the Sox2 motif (E > 0.05). DREME reports that this motif occurs in 18% of the Oct4 ChIP-seq peaks, but in only 9% of the Sox2 peaks. The Oct4 motif is thus significantly enriched in Oct4 ChiP-seq peaks compared with Sox2 peaks. Reversing the roles of the two ChIP-seq datasets, DREME finds the Sox2 motif (Fig. 2a, E = 10 ^{− 158}), but not the Oct4 motif. The Sox2 motif found occurs in 55% of Sox2 peaks, but in only 28% of the Oct4 peaks. These two results strongly suggest that both Oct4 and Sox2 often bind DNA in mES cells alone or with other partners, rather than exclusively with each other as a heterodimer.

3.2.2 Finding context-dependent motifs

It was recently reported that Oct4 may preferentially bind to a different motif with consensus ATGCGCAT when not binding near Sox2 (Mason et al., 2010). They used Oct4 and Sox2 ChIP-seq data from Marson et al. (2008) to create two Oct4 sequence datasets. The first dataset, Oct4woSox2, contains only Oct4 regions where there is no Sox2 peak within 5000 bp. The second dataset, Oct4wSox2, contains Oct4 regions that are within 50 bp of a Sox2 peak. Using the first and second datasets as positive and negative inputs, respectively, to their discriminative motif discovery algorithm CMF, they discovered the ATGCGCAT ‘Oct-only’ motif. They also noted that the ChIP-seq peaks containing this motif are enriched for a different set of transcription factors (nMyc, cMyc, E2f1 and Zfx1) than the Oct4 peaks containing a nearby Sox2 binding site (Nanog, Sox2 and Smad1), suggesting that the preferred in vivo binding motif of Oct4 is context dependent.

We applied DREME in an analogous way and found that it also discovers the ‘Oct-only’ motif when applied discriminatively to the Oct4woSox2 and Oct4wSox2 datasets. On these datasets, the ‘Oct-only’ motif is the most significant one found by DREME (E = 1e-47), and DREME runs about 5.4 times faster than CMF (74 versus 378 s). Interestingly, the ‘Oct-only’ motif is also among the 17 motifs found by DREME in the Chen et al. Oct4 ChIP-seq dataset, raising the total number of identifiable motifs in that dataset to five (Table 1, column 5).

3.2.3 Refining the search for the ChIP-ed factor's motif

Sometimes the most significant motif found by DREME is not that of the ChIP-ed TF, and it may not even be clear whether its motif has been found at all. For example, the most significant motifs found by DREME in the Nanog dataset are Oct4 and Sox2. This suggests that Nanog might sometimes be binding indirectly via one or both of these TFs. We wondered if a discriminative approach—looking for motifs enriched in the Nanog dataset relative to the Oct4 or Sox2 dataset—might increase our confidence in the motif we proposed as the Nanog DNA-binding motif (Fig. 1c). We, therefore, ran DREME using the Nanog mESC ChIP-seq dataset as the positive input, and each of the other two TFs' datasets (individually) as the negative input. The most significant motif in each of these two cases (Nanog versus Oct4, Nanog versuss Sox2) is indeed highly similar to the DREME motif we suggest above as the putative Nanog motif (Fig. 2c), and both are more similar to the previously published in vitro motif for Nanog [(C/A)(C/A)ATTA] than the motif found by DREME in non-discriminative mode. When we use TOMTOM to compare these two new putative Nanog motifs to our two reference motif databases plus the TRANSFAC motif database (Matys et al., 2006), each is most similar to a motif for Pbx1 (data not shown). So, as noted above, determining whether this motif is indeed the binding motif for Nanog, Pbx1 or a combination thereof, may require ChIP-seq data for Pbx1 in mES cells. However, the fact that this is the most highly enriched motif in the Nanog dataset relative to both the Oct4 and Sox2 datasets suggests that the motif is not merely due to the presence of Pbx1 binding sites in the Nanog dataset.