Sequence features and chromatin structure around the genomic regions bound by 119 human transcription factors

Previously unannotated motifs

We identified 11 high-confidence motifs that did not match any annotated motifs in the JASPAR or TRANSFAC repositories (Fig. 1C; Supplemental Figs. S2, S6). Among these motifs, UA1–UA5 are likely the canonical motifs for four TFs, and UA9 is likely the canonical motif for a factor that functions in H1-hESC cells. Supplemental Figure S7 shows that the sites of the previously unannotated motifs tend to have high evolutionary conservation and show distinct DNase I footprints.

UA1 was detected as the primary motif of three TFs (ZBTB33, BRCA1, and CHD2), as well as a secondary motif for ETS1. Because ZBTB33 is a zinc finger protein that binds methylated CpG dinucleotides (Yoon et al. 2003) and the center of UA1 contains CGCG, UA1 most likely is the canonical motif of ZBTB33. BRCA1 and CHD2 do not have a DNA-binding protein domain, suggesting that they bind ZBTB33 to perform their functions in DNA repair and genome maintenance. Indeed, the 936 ZBTB33 peaks that contain UA1 sites and the 321 BRCA1 peaks that contain UA1 sites have 312 peaks in common. Similarly, the 936 ZBTB33 peaks that contain UA1 sites and the 1022 CHD2 peaks that contain UA1 sites have 719 peaks in common.

UA2 was the primary motif for the PBX3 data set in GM12878, with 44.3% of the 7431 peaks containing at least one UA2 site. We did not identify any previously published description of the sequence motif of PBX3.

UA3 was the primary motif for the ZBTB7A data set in K562, occupying 80.1% of the 19,942 peaks. CTCF was identified as a secondary motif but with a lower enrichment than UA3 (Supplemental Fig. S6). Genecard indicates that the consensus sequence of ZBTB7A binding sites is 5′-[GA][CA]GACCCCCCCCC-3′, which is similar to the reverse complement of the UA3 consensus 5′-CNGAGACCCCNNCCC-3′. Furthermore, the motif derived from the in vitro protein binding microarray method for ZBTB7B (a paralog of ZBTB7A) is very similar to UA3 (Badis et al. 2009).

UA4 and UA5 were discovered in the THAP1 data set in K562. UA4 is a gapped motif, and it is an extended version of the motif previously reported for the THAP family of TFs (Sabogal et al. 2010). UA5 shares the “GGGC” half of UA4 but further extends it. Thus both UA4 and UA5 are likely the canonical motifs for THAP1.

UA9 was discovered as the primary motif for NANOG (in H1-hESC cells) and BCL11A (in H1-hESC cells but not in GM12878 cells). It does not resemble the previously identified NANOG motif (Chen et al. 2008). We also discovered UA9 as a secondary motif for five other TFs in H1-hESC cells (see Supplemental Fig. S6). We, therefore, suspect that UA9 is the canonical motif of a yet uncharacterized TF that functions in H1-hESC cells.

We also identified two motifs that allow alternative spacing: The two GATA3 half sites, AGAT and ATCT, can be either 3 or 4 bp apart, and the two half sites of the AP-1 motif can be either 1 or 2 bp apart. The variant spacing of AP-1 was previously detected by the in vitro protein binding microarray method (Badis et al. 2009; Supplemental Fig. S5), reflecting intrinsic flexibility of the two leucine zippers of the heterodimeric AP-1 TF. The variant spacing of GATA3 has not been reported previously. We identified extensions of four annotated motifs—CREB, ZNF143, GATA1, and CTCF (Fig. 1D; Supplemental Fig. S6). ZNF143-ext (Myslinski et al. 2006) and CTCF-ext (Ohlsson et al. 2001; Rhee and Pugh 2011) have been documented before. GATA1-ext is the motif for the TAL–GATA1 complex (Xu et al. 2003). The extension for CREB has not been reported.

Comparison of bound vs. unbound motif sites

Although the ChIP-seq peaks are highly enriched in motifs, there are still many motif sites outside peaks (unbound motif sites). For example, there are, on average, 430 times more unbound motif sites (sites outside peaks; FIMO P-value < 1 × 10⁻⁴) than bound motif sites (sites within peaks) for the TFs with ChIP-seq data in K562 cells. We asked whether there were any sequence or chromatin features that could distinguish bound sites from unbound sites (see Supplemental Methods for details). Indeed, we found that the regions surrounding bound sites were more DNase I hypersensitive and enriched in TF motifs, compared with the regions surrounding unbound sites, as shown in Supplemental Figure S8 for the five cell lines with the most ChIP-seq data sets, one heat map per cell line. The histogram of log₂ (enrichment) has a heavier right-side tail in all cell lines, indicating an overall enrichment among all pairwise comparisons (Supplemental Fig. S8). As expected, regions around bound A-box sites are enriched in B-box sites and vice versa, consistent with these sites being the TFIIIC motifs in tRNA genes (Oettel et al. 1998). The bound regions of most motifs are enriched in sites of the same motif. Several motifs such as NRF1 are enriched in the bound sites of the majority of motifs across the cell lines.

Cobinding and tethered binding between different TFs

Many eukaryotic genes are coregulated by multiple TFs in a cell-type–specific manner (Maston et al. 2006). For 70 of the 87 sequence-specific TFs, we discovered the canonical motifs as well as significant secondary motifs that were distinct from the canonical motifs of the TFs in question and that correspond to the canonical motifs of other TFs. Two scenarios may result in secondary motifs: Two TFs bind to neighboring sites (cobinding), or one TF protein binds to another that, in turn, binds to DNA (tethered binding). To distinguish between these scenarios, we computed the percentages of peaks in a ChIP-seq data set that contain sites for the canonical TF only, a noncanonical TF only, or both, and then we sorted the data sets by the percentages of peaks with only noncanonical motif sites (Fig. 2A; see Supplemental Table S3 for the underlying data). We reasoned that if sites of a noncanonical motif were frequently found to be in the same ChIP-seq peaks as canonical motif sites (hence, adjacent to them), the two TFs are likely to interact at the protein level and influence each other in binding to their DNA sites. Conversely, if the majority of the peaks contain only sites for noncanonical motifs, then tethered binding is a more plausible model. In this fashion, we identified 151 potential tethered binding and 104 cobinding sequence-specific TF pairs (255 in total). We then compared the pairs we discovered with experimentally detected pairs reported in a mammalian two-hybrid study (Ravasi et al. 2010) and in the BIOGRID database (Stark et al. 2006) and found evidence for physical interaction for 27 (10.6%) of the pairs. Eighteen of the 151 tethered binding predictions were validated in the mammalian two-hybrid data. We randomly picked 151 TF pairs for 5000 trials, and on average, 4.19 pairs were validated in the mammalian two-hybrid experiments (maximum 13 pairs), indicating that our predicted TF pairs were highly significant (P-value < 2 × 10⁻⁴). Thus, our results both recapitulated previously reported observations and revealed novel potential interactions that can be tested by experimentation (see Supplemental Table S3 for summary of all pairs).

Open in a separate window

Figure 2.

Interactions between TFs. (A) Different modes of interaction between TFs are shown. Each bar indicates the canonical TF and one noncanonical TF whose motifs were identified in the same ChIP-seq data set, and the red, blue, and black segments of the bar indicate percentage of peaks in the ChIP-seq data set that contain only canonical motif sites, only noncanonical motif sites, or both. Cartoons depict examples of different models for TF-TF interactions. (B) Circos plot (Krzywinski et al. 2009) on the left depicts pairs of motifs (connected by an arch) with significant distance preferences between their sites. The thickness of a connection is proportional to the normalized frequency of the pair. A connection is depicted as blue, black, or red when the motif pair is discovered in different data sets, the same data set, or both, respectively. The heat map on the right shows the distributions of distances between motif pairs. Each row is a motif pair in a particular ChIP-seq data set, and each column represents an edge-to-edge distance (from 0 bp to 99 bp). (C) Similar to B except showing motif pairs discovered in repetitive regions.

SP1 (or SP2) and NF-Y (heterodimer of NF-YA and NF-YB) constitute an example of cobinding. (SP1 and SP2 mostly bind to common sites in the genome, and their motifs are indistinguishable.) The SP2 ChIP-seq data set in K562 cells contains 3025 peaks, of which 2496 peaks contain SP1/2 sites and 1711 peaks contain NF-Y sites, sharing 1512 peaks. Furthermore, 1562 of the 1711 NF-Y site-containing peaks overlap with the peaks in the NF-Y ChIP-seq data set in K562 cells, confirming that these NF-Y sites are, indeed, bound by the NF-Y protein. Thus, SP2 and NF-Y prefer to cobind neighboring sites in the genome. Similarly, the SP1 ChIP-seq data sets in K562 and GM12878 cells indicate cobinding between SP1 and NF-Y. In the next section, we show that NF-Y and SP1/2 binding sites tend to be within 30 bp of one another. Indeed, Roder et al. showed that NF-Y and SP1 proteins bind to each other and that they cobind to the same promoter (Roder et al. 1999). In another example, YY1 was shown to interact with MYC, and cooperatively, the two TFs regulate the expression of the ITGA3 gene of the α3β1-integrin complex in human osteosarcoma cells (de Nigris et al. 2007). We found that the cobinding of YY1 and MYC was not limited to the ITGA3 promoter. In K562 cells, 484 MYC peaks (11.9%) contain both YY1 and MYC motifs, and these peaks are bound by YY1 as well. Novel predictions of cobinding TF pairs include ESRRA and HNF4 in HepG2 cells and NFKB and SPI1 in GM12878 cells (Supplemental Table S3).

Tethered binding facilitates combinatorial regulation by additional sequence-specific TFs and is necessary for recruiting chromatin remodelers and other regulatory proteins that do not bind DNA directly. We describe examples of tethered binding between sequence-specific TFs here, and discuss the tethering of non-sequence-specific TFs onto sequence-specific TFs in a subsequent section. ATF3 contains a basic leucine zipper DNA-binding domain whose canonical motif is the CREB motif. However, 51.2%, 48.6%, 46.9%, and 49.5% of the peaks in the four ATF3 ChIP-seq data sets (in GM12878, H1-hESC, HepG2, and K562 cells, respectively) contain USF sites but not CREB sites, and 98.35%–98.95% of these peaks overlap with peaks of USF1 or USF2 in the respective cell lines, indicating that they are bound by USF (USF is the heterodimer of USF1 and USF2). USF belongs to the basic helix-loop-helix leucine zipper family, and its motif does not resemble the CREB motif (Supplemental Fig. S2). Thus, the above analysis suggests that ATF3 tethers to USF, which binds DNA directly. Although ATF3 and USF have been reported to coregulate the same promoters (Runkel et al. 1991), it has yet to be shown that these TFs physically interact with each other independent of target DNA. As another example, in K562 cells, NFE2 binds 75.98% of its ChIP-seq peaks directly and tethers through USF for 15.59% of the peaks. Interestingly, the peaks that involve tethered binding show the strongest ChIP signal for this data set (Supplemental Fig. S6), perhaps because the NFE2 proteins that bind to DNA directly adopt a different conformation from those involved in tethered binding, and because the antibody used to perform the ChIP experiments has a lower affinity for NFE2 in the direct DNA-binding conformation. This is an exception because, for most TFs, the peaks that contain sites for the canonical motif show the strongest ChIP-seq signal (Supplemental Fig. S6). Other examples of tethered binding include SP1 tethering to HNF4 in HepG2 cells and STAT3 tethering to CEBPB in HeLa cells, both with literature support (Kardassis et al. 2002; Zhou et al. 2010). Novel predictions of tethering include TCF12 to FOXA and HNF4 in HepG2 cells, IRF1 to NF-Y in K562 cells, SREBF1 to RFX5 in HepG2 cells, and SIX5 to ZNF143 in GM12878 cells (Supplemental Table S3).

As a special case of cobinding, TFs that belong to the same protein family usually share identical or similar motifs and may compete for sites that match both motifs. MYC/MAX and USF (the USF1/USF2 heterodimer) both contain basic helix-loop-helix leucine zipper DNA-binding domains, but they do not cross-dimerize (Sawadogo et al. 1999). Their motifs share the CACGTG core, but the USF motif contains two additional nucleotides (gtCACGTG) that have a moderate sequence preference (Supplemental Fig. S2). The USF motif was discovered as a secondary motif in all five MAX ChIP-seq data sets, and 77.37%–92.75% of USF sites identified in the MAX data sets overlap with peaks in the USF1 or USF2 ChIP-seq data sets in the same cell line. These results suggest that USF and MYC/MAX compete for these sites. It was reported that both USF and MYC/MAX can bind an E-box motif in the promoter of the hamster cad gene, but only the binding of MYC/MAX is required for the transcription of cad (Boyd and Farnham 1997).

Distance and orientation preferences between the sites of cobinding TFs

Cobinding TFs bind to neighboring sites in the genome. For some TFs, multiple molecules of the same TF also can occupy neighboring sites. We asked whether these neighboring sites prefer to be on the same strand or opposite strands and whether they prefer to be in a specific range of distances. In addition to the analysis presented in the previous section, which compared the canonical motif with each noncanonical motif discovered in the same data set, we also compared motifs discovered in different data sets collected using the same cell line. In Figure 2B,C, we summarize the heterotypic and homotypic TF pairs that show statistically significant orientation or distance preferences separately in nonrepetitive and repetitive regions of the genome (the underlying data are in Supplemental Table S4). Out of the 78 motifs discovered from ChIP-seq data sets, 36 motifs (92 pairs; 62 heterotypic pairs and 30 homotypic pairs) are included in Figure 2B, suggesting that preferred arrangements of nearby TF binding sites are a common phenomenon. The neighboring sites for many heterotypic TF pairs (e.g., CTCF–NF-Y, ELF1–GABP, and FOXA–HNF4) as well as the neighboring homotypic sites of many TFs (e.g., AP-1, CTCF, and USF) show a strong preference for an edge-to-edge distance of <30 bp and varying degrees of preference for one orientation over the other. For example, neighboring NF-Y sites prefer to be in the same orientation. NF-Y also prefers one orientation to the other when cobinding with SP1, PBX3 (its motif is UA2), and USF. We hypothesized that these 92 TF pairs are more likely to represent protein–protein interactions than the TF pairs we identified in the previous section without testing for position or orientation preferences. Indeed, 14 heterotypic pairs and 17 homotypic pairs (33.7%) were detected in the aforementioned mammalian two-hybrid study (Ravasi et al. 2010) or in the BIOGRID database (Stark et al. 2006).

TFs tend to bind gene-rich regions of the genome due to their role in regulating target gene expression (Carroll et al. 2006). Nonetheless, repetitive elements are known to harbor functional TF binding sites, especially when such elements occur near genes. We systematically compared our compilation of TF binding sites with all repeats annotated in the human genome, and the results are summarized in Figure 3A. We confirmed the previously reported enrichment of STAT1, NF-Y, and CTCF binding sites in various repetitive elements (Bourque et al. 2008; Schmid and Bucher 2010), and we uncovered many more TFs whose binding sites are enriched in certain repetitive elements, e.g., UA1 sites in THE1B and THE1D retrotransposons. It was shown that a long terminal repeat (LTR) region of the THE1D retrotransposon was recruited as an alternative promoter for the human IL2RB gene and that the activity of this alternative promoter is regulated by DNA methylation (Cohen et al. 2011). The UA1 motif we identified in ZBTB33 peaks contains a prominent CGCG center (Fig. 1C) and ZBTB33 is known to bind methylated CpG dinucleotides (Yoon et al. 2003), raising the interesting possibility that the THE1B/D retrotransposons spread ZBTB33 binding sites across the genome and that the regulation of the newly recruited target genes can be modulated by the DNA methylation mechanism. Figures 2C and ​and3B3B summarize all motif pairs that show statistically significant distance or orientation preference in repetitive regions of the genome. The NF-Y–USF site pairs that typically have an end-to-end distance of 5–6 bp are nearly all located in the MLT1 family of retrotransposons. Similarly, the NF-Y–NF-Y site pairs at a 9-bp distance are found most often in LTR12 retrotransposons. There are 181 copies of the MLT1J transposon in the genome that contain sites for the NF-Y, USF, and ZNF143 motifs simultaneously, bound directly by NF-Y, USF, and ZNF143 TFs, respectively. The relative distance among the sites are nearly invariant (Fig. 2C), indicating recent duplications of MLT1J. Our results suggest a mechanism whereby retrotransposons amplify functional TF site pairs across the genome through transposition, potentially bringing new genes under the regulation of those TFs.

Open in a separate window

Figure 3.

Binding sites of certain TFs or TF pairs are enriched in repeats. (A) Enrichment of TF binding sites in repetitive elements. The redness of each grid point is proportional to the negative logarithm of enrichment P-value. Repetitive elements are color-coded by family. (B) Enrichment of motif pairs that strongly prefer a narrow distance range in various repetitive elements (Fig. 2C).

Cell-type–specific binding of sequence-specific TFs

The majority of the ENCODE ChIP-seq data was produced using five cell lines: K562, GM12878, HepG2, H1-hESC, and HeLa. Integrating ChIP-seq data with RNA-seq data for these five cell lines, we asked whether genes that are preferentially expressed in a given cell line (defined by the average expression level in one cell line being more than 10-fold higher than that in any of the remaining four cell lines) show enriched TF binding sites in the corresponding cell line. This is, indeed, the case for a large fraction of genes, and Figure 4A shows five examples, one per cell line. (1) FCER2 (the low-affinity receptor for IgE) is a key gene for B-cell function. It is highly and specifically expressed in GM12878. Its promoter region and gene body are bound by nine TFs in GM12878, including SPI1. (2) The G protein-coupled receptor GPRC5A plays a role in epithelial cell differentiation. It is highly and specifically expressed in HeLa cells, and accordingly, its promoter region and gene body are bound by seven TFs in HeLa cells. (3) The Abd-B homeobox family member HOXB9 is a sequence-specific transcription factor. It is highly and specifically expressed in K562 cells, and accordingly, its promoter regions and gene body are bound by seven TFs including GATA1-TAL1 in K562 cells. (4) SERPINA1 encodes a serine protease inhibitor, and defects in this gene can cause liver diseases. It is four orders of magnitude more highly expressed in HepG2 than in the other four cell lines. FOXA, HNF4, RXRA, TCF7L2, and eight other TFs bind near this gene in HepG2 but not in other cell lines. (5) {"type":"entrez-nucleotide","attrs":{"text":"AC104304","term_id":"21490239","term_text":"AC104304"}}AC104304 encodes for a putative teratocarcinoma-derived growth factor that plays an important function in embryonic development. It is highly expressed in H1-hESC and bound by eight TFs, including NANOG.

Open in a separate window

Figure 4.

Cell-type–specific binding of sequence-specific and non-sequence-specific TFs. (A) Abundant TF binding sites are observed near cell-line–specific transcripts. Binding sites are shown as vertical bars and colored by cell line (dark blue for K562, red for HepG2, brown for GM12878, green for H1-hESC, and cyan for HeLa-S3). (Bottom, right) Expression levels (in RPKM) for example cell-line–specific transcripts across the five cell lines with the most ChIP-seq data. (B) Secondary motifs identified in the ChIP-seq data sets of sequence-specific TFs and their enrichment in the ChIP-seq peaks of non-sequence-specific TFs. The five cell lines are indicated with color-coded squares, noncanonical motifs of sequence-specific TFs are shown in pink circles and a solid line connecting each motif to the respective cell line. The thickness of the solid line is proportional to the normalized frequency in which a noncanonical motif is discovered in a particular cell line. Non-sequence-specific TFs are shown in diamonds whose colors match the color of the cell line if there is a ChIP-seq data set of the TF in that cell line. Dashed lines connect non-sequence-specific TFs and noncanonical motifs, indicating that a noncanonical motif of a sequence-specific TF is enriched in the ChIP-seq peaks of the non-sequence-specific TF. (Four insets) Expression profiles of sequence-specific TFs whose canonical motifs are found to be specific to a cell type and the TF binding sites around the genes that encode these TFs in the appropriate cell line. The expression levels in each cell line are assigned a similar color as the cell line. For four cell lines, two biological replicates were available for RNA-seq data; hence, there are two bars for each of these cell lines. Only one biological replicate was available for H1-hESC.

We then asked whether the noncanonical motifs we discovered also reflect cell type specificity. Figure 4B plots the noncanonical motifs (circles) detected in the ChIP-seq data sets of sequence-specific TFs for each of the five cell lines (squares) with the most ENCODE ChIP-seq data sets. Cell-line–specific, noncanonical motifs are placed close to their respective cell lines in Figure 4B. We defined cell-line–specific motifs as those that were discovered three times more often in one cell line than in any other cell line. The remaining noncanonical motifs are placed in the center of the figure, and these motifs correspond to TFs that cooperate with other sequence-specific TFs across multiple cell lines. The thickness of the solid line connecting a noncanonical motif to a cell line indicates the proportion of data sets in that cell line that revealed the motif as a noncanonical motif.

We highlight several motifs that were frequently discovered as noncanonical motifs in a particular cell line. (1) PU.1 was most frequently discovered in GM12878 cells. Its corresponding TF SPI1, a member of the ETS family, activates gene expression during myeloid and B-lymphoid cell development. The SPI1 gene is expressed in both GM12878 and K562 cells (RPKM = 4 and 9, respectively), but not in the other three cell lines (RPKM < 0.1). On the other hand, another member of the ETS family, SPIB, is only expressed in GM12878 cells, and the SPIB gene shows extensive TF binding sites specifically in GM12878 cells (bottom inset in Fig. 4B). SPIB and SPI1 have the same canonical motif (PU.1) and are both essential for B cell development (Sokalski et al. 2011). (2) GATA1 was the most frequently discovered noncanonical motif in K562 cells. It is bound by the GATA family of TFs, which are essential for erythroid development by regulating the fetal-to-adult switch of hemoglobin production (Takahashi et al. 1997). The GATA1 gene is highly expressed in K562 cells but not in the other four cell lines and shows extensive binding sites only in the K562 cell line (top right inset in Fig. 4B). (3) FOXA and HNF4 are the most frequently identified noncanonical motifs in HepG2 cells. Their corresponding TFs (FOXA1 and HNF4) are activators of many liver-specific genes and are essential for hepatocyte function (Lemaigre and Zaret 2004). Both the FOXA1 and HNF4 genes are more than 10-fold more highly expressed and show more extensive TF binding sites in the HepG2 cell line than in the other four cell lines (FOXA1 is shown in the middle right inset in Fig. 4B). (4) The SOX2-OCT4 combined motif was the most frequently identified noncanonical motif in H1-hESC cells. OCT4 is the canonical motif of POU5F1, a POU homeodomain-containing TF required for embryonic stem cell pluripotency. Their corresponding TFs (POU5F1 and SOX2) form a protein–protein complex and are required for embryonic stem cell pluripotency (Kashyap et al. 2009). Both POU5F1 and SOX2 are exclusively expressed in H1-hESC cells and extensively regulated by a large number of TFs, including by themselves (POU5F1 is shown in the lower left inset of Fig. 4B).

Tethered binding of non-sequence-specific TFs

In Figure 4B, we also included all non-sequence-specific TFs (diamonds) for which there are ChIP-seq data in these cell lines. Dashed lines connect non-sequence-specific TFs to the motifs discovered in their ChIP-seq peaks. Two non-sequence-specific TFs show cell-line–specific enrichment in motifs: the enhancer-binding protein EP300 and the histone deacetylase HDAC2. There are seven data sets for EP300 in seven different cell lines and three data sets for HDAC2 in three different cell lines. Distinct motifs were found in different cell lines: SPI1 for EP300 in GM12878 cells; GATA1 (and GATA1-ext) for both EP300 and HDAC2 in K562 cells; FOXA and HNF4 for HDAC2, and FOXA and TCF7L2 for EP300 in HepG2 cells; SOX2-OCT4 and UA9 for HDAC2, and TEAD1 for EP300 in H1-hESC cells; and CEBPB, AP-1, and CREB for EP300 in HeLa cells. As described in the previous section, many of these motifs were most frequently and specifically observed as secondary motifs for sequence-specific TFs in the respective cell lines. Because non-sequence-specific TFs do not bind DNA directly, they tether onto sequence-specific TFs to bind target DNA. EP300 is known to interact with AP-1 and CEBPB (Mink et al. 1997; C-C Wang et al. 2007) and HDAC2 with TAL1-GATA (the motif is GATA1-ext) (Hu et al. 2009). Our results highlight that the interactions of EP300 and HDAC2 with sequence-specific TFs are highly cell type dependent.

We further analyzed the subsets of peaks of the sequence-specific TFs through which EP300 and HDAC2 might tether in a cell-type–specific fashion. For example, in GM12878, 4.9% (1745 out of 35,821) of the SPI1 peaks were associated with EP300 (i.e., overlapping with the EP300 peaks in the same cell line by at least 1 bp). We asked whether these peaks differed from the remaining 95.1% of SPI1 peaks by performing two types of analysis. First, we tested whether these EP300 or HDAC2-associated peaks were differentially enriched in any of our collection of 78 motifs (with enrichment defined in the Supplemental Methods). Second, we performed de novo motif finding on the top 500 of the subsets of peaks (ranked by ChIP-seq signal) using MEME-ChIP.

For EP300-associated peaks, we observed differential enrichment in a number of known motifs (Supplemental Fig. S8), and a subset of these motifs were also found by MEME-ChIP. For example, the 1745 SPI1 peaks that were associated with EP300 were more enriched in 18 motifs than all SPI1 peaks in GM12878 (cf. the PU.1–EP300 row and the PU.1 row in the GM12878 heat map in Supplemental Fig. S8), and MEME-ChIP identified AP-1, SP1, and RUNX1, in addition to the canonical motif of SPI1 (PU.1), which were among the 18 motifs. Similar results were observed for the GATA1 peaks that were associated with EP300 in K562 cells and the FOXA1 peaks that were associated with EP300 in HepG2 cells (see Supplemental Fig. S8 for results on enrichment; the MEME-ChIP results are not shown).

For HDAC2-associated peaks, we did not observe differential enrichment of any annotated motifs (for example, cf. the GATA1 row and the GATA1–HDAC2 row in the K562 heat map of Supplemental Fig. S8). For the HNF4 peaks associated with HDAC2 in HepG2 cells, MEME-ChIP only identified the HNF4 motif. Intriguingly, MEME-ChIP identified UA9 from the subset of POU5F1 peaks associated with HDAC2 in H1-hESC cells (in addition to POU5F1's canonical motif SOX2-OCT4). In comparison, MEME-ChIP did not identify UA9 from the overall top 500 POU5F1 peaks but identified UA9 from the overall top 500 HDAC2 peaks in H1-hESC cells (Supplemental Fig. S6), indicating that UA9 is specific to the regions bound by HDAC2. Indeed, 199 of the 473 HDAC2-associated POU5F1 peaks (42.1%) contain significant UA9 sites, compared with 23.1% in the POU5F1 peaks not associated with HDAC2 (Supplemental Fig. S9). Furthermore, MEME-ChIP identified another previously unannotated palindromic motif (Supplemental Fig. S9). This motif was significantly enriched in the HDAC2-associated GATA1 peaks in K562 cells compared to a set of randomly chosen genomic regions with matching GC% and lengths (P-value = 1.1 × 10⁻⁸, using only the 1078 peaks that were not used as input to MEME-ChIP). This motif was distinct from all of the 78 motifs that were identified in the overall top 500 ChIP-seq peaks; thus, we named it UA12 and added it to our collection (Supplemental Fig. S2). Indeed, 36.9% of the HDAC2-associated GATA1 peaks in K562 cells contained UA12 sites (Supplemental Fig. S9), which was higher than the percentage of GATA1 peaks not associated with HDAC2 (32.1%) and the flanking regions (26.1%).

The ChIP-seq peaks of several non-sequence-specific TFs were enriched in the same motifs, regardless of cell line. The best example is the chromodomain helicase DNA-binding protein CHD2, for which the UA1 motif (likely the canonical motif of ZBTB33) was prominently discovered in GM12878, HepG2, and K562 cells, suggesting that CHD2 functions by interacting with ZBTB33. REST is highly enriched in the ChIP-seq peaks of SIN3A in H1-hESC cells. Accordingly, SIN3A is known to associate with REST and repress neuronal genes in nonneuronal cells (Huang et al. 1999).

The ChIP-seq peaks of many TFs are flanked by positioned nucleosomes

After analyzing the motif content of TF binding peaks, we set out to investigate the chromatin structure around these peaks. It is well known that transcription has a profound impact on nucleosome occupancy: Active TSSs in all eukaryotes are flanked by an upstream nucleosome-depleted region and several well-positioned downstream nucleosomes (Radman-Livaja and Rando 2010). We previously reported that the binding sites of the insulator binding protein CTCF were flanked by an array of strongly positioned nucleosomes, shown as a periodic oscillatory pattern in the average nucleosome occupancy profile centered on CTCF binding sites (Fu et al. 2008). Another study showed that NRSF (also called REST) binding sites are flanked by positioned nucleosomes in CD4+ T cells, CD8+ T cells, and granulocytes (Valouev et al. 2011).

In order to investigate where TF binding peaks were located with respect to nucleosomes, we computed an average nucleosome occupancy profile centered on the peak summits of each TF with available ChIP-seq data in GM12878 or K562 cells (Fig. 5A,B for YY1 in GM12878 cells; Supplemental Fig. S10 for all data sets). We had ChIP-seq data for 51 TFs in GM12878 cells, 73 TFs in K562 cells, and 32 TFs in both cell lines. Some TFs were tested by multiple labs in the same cell line, and we included all these data sets. To account for the impact of transcription, we computed the average nucleosome profile anchored on TSS-proximal and TSS-distal peak summits separately. Nucleosome profiles anchored on TSS-proximal peaks were oriented such that the nearest transcript is downstream from the anchor. We further stratified peaks in each data set as top, middle, and bottom thirds according to the ChIP-seq signal, reflecting the extent to which a peak is bound by the TF (averaged over a population of cells). We distinguish nucleosome occupancy and nucleosome positioning, with occupancy defined as the area under the occupancy profile and positioning defined as the regularity of the oscillatory pattern in the occupancy profile. Thus, the regions around TSS-proximal summits tend to show lower nucleosome occupancy and lower nucleosome positioning than regions around TSS-distal summits (cf. Fig. 5A,B; similarly the proximal and distal panels in Supplemental Fig. S10; note the difference in the y-scale). This difference may reflect the effects of RNA polymerase on chromatin structure (Weiner et al. 2010). Within the proximal and distal categories, the top, middle, and bottom third peaks, which correspond to the peaks with strongest, medium, and weakest TF binding, tended to show the greatest, medium, and weakest nucleosome positioning (Fig. 5A,B; Supplemental Fig. S10). Thus regions that are more strongly bound by TFs are flanked by better-positioned nucleosomes.

Open in a separate window

Figure 5.

Chromatin structure and GC content around TF binding regions. (A,B) Nucleosome occupancy profiles anchored on the summits of TSS-proximal (A) and TSS-distal (B) peaks of YY1 grouped by ChIP-seq signal strength: top (green), middle (red), and bottom (blue) third peaks in terms of ChIP-seq signal. Nucleosome depletion for the top third peaks is shown as D in each panel. (C) Distribution of nucleosome depletion “D” across all tested TFs, with peaks stratified according to TSS proximity (proximal or distal) and ChIP-seq signal strength (top, middle, or bottom third). P-values for pairwise comparisons based on paired Wilcoxon rank-sum tests are: P1 = 8.2 × 10⁻¹⁷, P2 = 7.6 × 10⁻²¹, P3 = 3.8 × 10⁻²³, P4 = 8.8 × 10⁻¹⁰, P5 = 1.1 × 10⁻⁹, P6 = 1.1 × 10⁻¹¹, and P7 = 6.6 × 10⁻²². (D) TF binding is correlated with significantly more nucleosome depletion than TSS. Wilcoxon rank-sum test P-values are shown separately for GM12878 and K562 cells. For the box plots in C and D, only those subcategories with 200 or more peaks are included, and whiskers represent the 1.5 inter-quartile range. (E) Nucleosome occupancy genome-wide is correlated with GC%. The smoothed density scatter plot contains 40,000 data points; each data point is a randomly chosen 250-bp region of the human genome. (Black dots) Those regions that overlap with ChIP-seq peaks. (Black line) Least square fit. Pearson correlation coefficient = 0.62; P-value < 2.2 × 10⁻¹⁶. (F) Comparison of in vivo (green) and in vitro (black) nucleosome occupancy profiles around peak summits of YY1. GC% profile around the same summits is plotted in orange. Note elevated GC% at summit coincides with high in vitro nucleosome occupancy and low in vivo nucleosome occupancy.

The cohesin components SMC3 and RAD21 (Supplemental Fig. S10) show the most striking patterns of positioned flanking nucleosomes, similar to what we previously reported for CTCF (Fu et al. 2008), to which these factors bind (Parelho et al. 2008). Two other TFs—CTCFL (a paralog of CTCF) and ZNF143 (a zinc-finger protein with a long motif)—also show striking patterns of positioned flanking nucleosomes. The binding sites for ∼70% of the tested TFs are flanked by positioned nucleosomes (Supplemental Fig. S10), indicating that this is a general phenomenon for most TFs. To quantify the regularity of nucleosome positioning around TF binding sites, we applied Fourier transforms to the nucleosome occupancy profiles, yielding power spectra. The height of the power spectrum at the spatial frequency corresponding to the nucleosomal repeat length was used as an indicator of how periodically nucleosomes were positioned (an example power spectrum is shown in Supplemental Fig. S11B). The spectrum height correlated significantly with the extent of positioning of the –1 and +1 nucleosomes (measured as the maximum of the nucleosome occupancy profile) (Supplemental Fig. S11C). Thus, how well the –1 and +1 nucleosomes are positioned strongly predicts how periodically the flanking nucleosomes are positioned.

Most TFs bind at GC-rich, nucleosome-depleted, and DNase I-accessible regions

The nucleosome occupancy profile dips at the peak summits of most TFs (Fig. 5A,B; Supplemental Fig. S10), indicating that TFs prefer to bind nucleosome-depleted regions or that the binding of a TF excludes nucleosomes. In the vicinity of TSS-proximal summits, lower nucleosome occupancy is seen in the direction of transcription than upstream of transcription. We define nucleosome depletion as the amount that nucleosome occupancy dips at the peak summit, as compared to the nucleosome occupancy at 2 kb from the summit (considered as background). TSS-proximal summits show significantly greater nucleosome depletion than TSS-distal summits (Fig. 5C). It is well known that the binding of the transcriptional machinery to the TSS excludes nucleosomes to a considerable extent (Radman-Livaja and Rando 2010). Indeed, average nucleosome occupancy anchored on the TSS shows an overall loss of nucleosomes (Supplemental Fig. S12). Interestingly, we observed that TSS-proximal TF peak summits show a significantly greater depletion in nucleosome occupancy than do TSSs (Fig. 5D). The median nucleosome depletion at the summits of TSS-proximal peaks is 0.56 for GM12878 cells and 0.59 for K562 cells, significantly greater than the maximal nucleosome depletion around TSS (0.42 for GM12878 cells and 0.48 for K562 cells; Wilcoxon rank-sum test P-value = 7.1 × 10⁻²⁸ and 1.1 × 10⁻²², respectively). Within the proximal and distal categories, the top, middle, and bottom third peaks showed greatest, medium, and weakest nucleosome depletion, respectively (Fig. 5C). This result indicates that TFs and nucleosomes compete for the genomic DNA and that stronger TF binding is correlated with greater nucleosome depletion, above and beyond the effect of transcription.

The peaks of seven TFs (BRF2, HDAC8, TRIM28, SETDB1, WRNIP1, ZNF274, and ZZZ3) do not show nucleosome depletion, nor are these peaks flanked by well-positioned nucleosomes, indicating these TFs tend to bind nucleosomal DNA (Supplemental Fig. S10). Three of these TFs function with each other to repress transcription. SETDB1 is a histone methyltransferase that catalyzes H3K9me3, which signals for the silencing of euchromatic genes (Bilodeau et al. 2009). TRIM28 (commonly known as KAP1) represses transcription by recruiting SETDB1 (C Wang et al. 2007). ZNF274 is a zinc-finger containing TF that binds to the 3′ end of zinc-finger coding genes and recruits chromatin-modifying proteins such as SETDB1 and TRIM28, which leads to transcriptional repression (Frietze et al. 2010). HDAC8 is a histone deacetylase and a transcriptional repressor. We caution that the HDAC8 ChIP-seq data set had only 287 peaks. BRF2 is a component of the RNA Pol III machinery (Moqtaderi et al. 2010). WRNIP1 (commonly known as WHIP) regulates DNA synthesis. ZZZ3 is a component of the ATAC complex and a histone H3 acetyltransferase and has been shown to acetylate both free and nucleosomal H3 (Wang et al. 2008).

We next asked whether the intrinsic DNA sequence properties of ChIP-seq peaks contribute to nucleosome depletion. In an earlier study, we reported a strong correlation between GC-rich sequences and their potential to form nucleosomes (Peckham et al. 2007). In vitro data also indicate that GC-rich sequences promote nucleosome formation (Valouev et al. 2011). Indeed, there is positive correlation between nucleosome occupancy and GC content for randomly chosen 250-bp regions of the genome (r = 0.62 and P-value < 2.2 × 10⁻¹⁶) (Fig. 5E). Many of those regions that overlap ChIP-seq peaks (Fig. 5E, black dots) are located above and to the left of the best-fit line, indicating that they have high GC% and low nucleosome occupancy. Compared with the average GC content of 40% in the human genome, ChIP-seq peaks are considerably more GC rich (61 ± 5% for TSS-proximal peaks and 53 ± 6% for TSS-distal peaks across the TFs). The high GC content may be due to the GC-richness of some TF motifs, but the motif sites are much smaller than peaks (8–21 bp vs. ∼250 bp), and we found similar GC patterns around TF summits without a motif site (data not shown). We conclude that TFs tend to bind GC-rich regions in the genome, regardless of the distance from the TSS. These results are seemingly contradictory—GC content is highly predictive of sequences that promote nucleosome formation, yet the GC-rich sequences surrounding TF binding sites are nucleosome-depleted in vivo.

To determine whether TF binding sites are, indeed, favorable sites for nucleosome formation, we used recent data from in vitro reconstitution of human genomic DNA into nucleosomes (Valouev et al. 2011) to construct in vitro nucleosome occupancy profiles around ChIP-seq peaks, confirming that in vitro nucleosome occupancy is much higher on the peak compared with flanking regions for the vast majority of TFs (Fig. 5F for YY1; Supplemental Fig. S13 for all TFs). Thus, TFs or their cofactors (such as chromatin remodelers) prevent the formation of nucleosomes or evict nucleosomes at these GC-rich locations of the genome.

Chromatin structure around cell-line–specific TF binding regions

In order to further investigate the relationship between TF binding and chromatin structure, we examined two sets of cell-line–specific ChIP-seq peaks for each TF—the set of peaks detected in GM12878 but not in K562 and the set of peaks detected in K562 but not in GM12878. We computed nucleosome occupancy profiles and DNase I cleavage profiles anchored on the summits of these two sets of peaks separately in each cell line (Fig. 6A,B for YY1; Supplemental Fig. S14 for all TFs). Strikingly, the peaks that were occupied by a TF in GM12878 (or K562) but not occupied by the TF in K562 (or GM12878) tend to be occupied by a nucleosome in K562 (or GM12878), similar to the in vitro nucleosome profiles of these peaks (Supplemental Fig. S15). Accordingly, the increase in nucleosome occupancy is reflected in decreased DNase I cleavage in K562 (or GM12878). For many TFs, the ChIP-seq peaks that were occupied by a TF in GM12878 (or K562) but not occupied by the TF in K562 (or GM12878) were no longer flanked by positioned nucleosomes in K562 (or GM12878), yet positioned nucleosomes were observed for other TFs, albeit at a lesser extent of positioning than the nucleosomes flanking TF-occupied peaks (Supplemental Fig. S14). Thus, for the same set of genomic sequences in two cell lines, TF binding level deviates from thermodynamic preference for nucleosome formation—TF binding either was enabled by, or caused, cell-type–specific depletion of nucleosomes from intrinsically favorable genomic locations.

Open in a separate window

Figure 6.

Chromatin structure around YY1 ChIP-seq peaks occupied differentially between GM12878 and K562. (A) Nucleosome occupancy profiles (solid lines) and DNase I cleavage profiles (dashed lines) anchored on the summits of YY1 peaks in GM12878 but not in K562. Note the average nucleosome occupancy at these peaks (x = 0) is lower in GM12878 than in K562, while the average DNase I cleavage at these peaks is higher in GM12878 than in K562. (B) Same as A, but around the summits of YY1 peaks in K562 but not in GM12878. (C) Nucleosome occupancy profiles in K562 anchored on the summits of the ChIP-seq peaks occupied by YY1 in GM12878 but not in K562. These 11,079 peaks were divided into two groups: 6754 peaks were bound by one or more TFs in K562 (dashed line), and 4325 peaks were not bound by any TF for which we had ChIP-seq data in K562 (solid line). Note high nucleosome occupancy at the summits of the unoccupied peaks (x = 0) and the lack of positioned nucleosomes flanking the unbound peaks, in sharp contrast to the lack of nucleosome occupancy at the peak summits and well-positioned nucleosomes flanking the peaks bound by other TFs. (D) Same as C, but around the summits of the ChIP-seq peaks occupied by YY1 in K562 but not in GM12878.

We further partitioned the peaks occupied by a TF in GM12878 (or K562) but not occupied by the TF in K562 (or GM12878) into two subsets: group 1 peaks that overlapped with one or more ChIP-seq peaks of any other TF tested in K562 (or GM12878), and group 2 peaks that did not overlap any ChIP-seq peaks in K562 (or GM12878). For the vast majority of the TFs, nucleosome occupancy profiles for the group 2 peaks show high nucleosome occupancy on the peak and no positioned nucleosomes flanking the peak. In contrast, the group 1 peaks show nucleosome depletion on the peak and well-positioned nucleosomes flanking the peak (Fig. 6C,D; Supplemental Fig. S16). The ChIP-seq data we have only cover up to 10% of TFs in a particular cell line, thus group 2 peaks can still be bound by other TFs for which we had no data, which could account for any residual pattern of nucleosome positioning. These results further strengthen the correlation between TF-binding and flanking positioned nucleosomes and indicate that such correlation can be regulated in a cell-type–specific manner.

De novo sequence motif discovery in ChIP-seq peaks

We constructed a de novo motif discovery pipeline as illustrated in Supplemental Figure S1 and described in detail in Supplemental Methods. Briefly, we separated the peaks in each ChIP-seq data set into independent training and testing sets to ensure the quality of the motif discovery. We used the MEME-ChIP software suite (Machanick and Bailey 2011) to identify enriched sequence motifs in the [–50 bp, +50 bp] window around the summits of the top 500 peaks (the training set) for each ChIP-seq data set. We asked MEME to report up to five significant motifs per data set. Then, we performed further testing to ensure the quality of the motifs, as described in Supplemental Methods. In total, we identified 1092 motifs. We compared the discovered motifs with annotated motifs in the JASPAR (Bryne et al. 2008) and TRANSFAC (Matys et al. 2003) databases using TOMTOM (Gupta et al. 2007). We manually merged similar motifs and identified 79 distinct motifs, which included 11 previously unannotated motifs (UA1–UA11) that we felt most confident about because they were highly enriched in one or more ChIP-seq data sets or were supported by the literature. We further identified another motif, UA12, in the subset of HDAC2-associated GATA1 peaks in K562.

Distance and/or orientation preferences between motif sites

For all pairs of motifs, we computed the edge-to-edge distance and relative orientation and computed a P-value using the Kolmogorov-Smirnov test, as described in Supplemental Methods. The analysis was also applied to peaks in nonrepetitive regions and repetitive regions of the genome separately. For motif pairs in repetitive regions that had exactly the same distance as the mode distance of the observed distribution, we tested their enrichment in various repetitive elements using a hypergeometric test. All P-values were corrected for multiple testing, and an FDR cutoff of 0.025 was applied throughout.

We thank ENCODE data producers for the rich data. We thank consortium members and Mike Pazin for insightful discussions and critical comments on the paper. We thank John Stamatoyannopoulos for bringing to our attention the distinction between tethered binding and cobinding of transcription factors. This work was funded by NIH grants U01 HG004561, U01 HG004695, and R01 GM103544, and NSF grant DBI-0850008.

Author contributions: Z.W. designed and supervised the analysis. J.W., J.Z., and S.I. performed the majority of the analysis. XY.L. engineered factorbook.org. The remaining authors contributed to the analysis. Z.W., J.W., J.Z., and S.I. wrote the manuscript. Other authors contributed to editing the manuscript.