Binding sites detection

Binding sites from ChIP-seq were called by Model-based Analysis for ChIP-Seq (MACS) (27) with default parameters for ERα, FoxA1, GATA3, AP2γ, Pol2 and p300. Enrichment regions for histone modifications were called by broad peak options with parameter setting ‘–nomodel –nolambda’ as suggested by Feng et al. (28).

EE, PP and enhancer–promoter interactions classification

Distance between two ChIP-seq binding peaks was calculated as the genomic distance between their summits called by MACS (27). Distance between two genomic regions was defined as the genomic distance between their mid-points. Density estimation was called by the R density function with default parameters. Peak position of a uni-modal distribution was chosen as the point with the highest density value. For multi-modal distribution, we first fitted the data with a Mixture Gaussian Distribution by R package mixtools (29) and chose the estimated expectation of each component as peak position.

We used log₂-transformed read-count ratio between H3K4me3 and H3K4me1 ChIP-seq data in the Pol2 binding sites to classify Pol2 ChIA-PET interaction clusters into three types: EE interactions, PP interactions and enhancer–promoter (EP) interactions. Windows for computing read counts were selected as ±1.5-kb region around the peak summits. Binding sites whose summits were closer than 1.5 kb were merged and the mid-point of merged regions was chosen as the new peaks’ ‘pseudo summit’. The data were then fitted with Mixture Gaussian Model and a Bayesian posterior probability of 0.5 was set to determine whether a Pol2 binding site was promoter-like or enhancer-like. PET clusters with both ends classified as enhancer-like were called EE interactions, both promoter-like were called PP interactions and the rest were called EP interactions. For ERα ChIA-PET data, we applied this classification in those interactions that overlapped with H3K4me1 or H3K4me3 peaks at both ends. The subsequent steps were similar to those of Pol2 ChIA-PET data.

Histone modification, DNase-seq and GRO-seq profiles surrounding ChIP-seq binding peaks

Each tag was extended by 200 bp along the read direction. We selected peak summit as the center for signal alignment. Each region was divided into 25-bp sub-regions and its read coverage was computed as the read counts covering this sub-region. The difference of read coverage between ChIP data and input data was computed as the log₂-transformed read-coverage ratio at each sub-region. DNase-seq data were not extended and only the 5′ end was used as aggregation inputs. GRO-seq data were divided into two groups according to the strand of the reads and then processed in the same manner as histone modification data.

To classify whether an ERBS was an ERE-containing one, we used STORM (30) to scan ±200-bp region around the summit of ERBS for ERE sequence motif, with P-value 1e−4.

Model to predict laERBSs

ERBSs were extracted from ChIP-seq data. To predict laERBSs, three types of features were computed:

1. F_i,j: log_e-transformed read counts for ChIP-seq data of protein j (or DNase-seq data) in ERBS i, with a window size of 400 bp centered around ERBS peak summit;

2. D_i: log_e-distance between the neighboring ERBS of ERBS i and ERBS i;

3. H_i,j: the differences between log₂-transformed ratio of read-coverage against input in central region (±100-bp region relative to the peak summit) to average of read-coverage against input in the two flanking regions (−400 bp to −200 bp and +200 bp to +400 bp relative to peak summit) of ERBS i for histone modification j.

We selected the 1822 ERBSs that related to 903 interactions identified in both two ChIA-PET experimental replicates as foreground training set. An equal number of ERBSs that did not overlap with any interactions in either replicate were randomly selected as background training set. We used a logistic classifier to perform the classification, i.e.

in which E_i is the indicator whether ERBS i is a laERBS, P is a probability function and k₀,a_i,j,b_i,c_i,j are the model parameters.

After training, we chose the most significant features that showed improvement over ERα ChIP-seq read counts in ERBSs to fix the final classifier. Then ERBSs were filtered by setting the threshold of the logistic classifier to be 0.2 to find putative laERBSs. Summits of these laERBSs closer than 3 kb were merged and the mid-point was selected as the new ‘pseudo summit’. This resulted in ~15 000 candidate anchors.

Model to predict ERBS interaction clusters

To predict ER-associated interactions, the candidate anchors were paired with each other to form candidate ERBS pairs. Those pairs that crossed topological domain boundaries in h1-ESC were excluded since domain boundaries were roughly invariant across cell lines, and long-range interactions were largely restricted within topological domains (31). To predict interactions, two types of features were computed for each candidate pair:

1. : the sum of log_e-transformed jth ChIP-seq (or DNase-seq) read counts for each candidate ERBS pair (i₁,i₂), with a window size of 3 kb for each end (±1.5-kb region relative to the peak summit);

2. : the log_e-distance and inverse distance between each ERBS pair (i₁,i₂), since previous study (32) suggested that probability of long-range interaction was not monotonically related to genomic distance.

Eight hundred ERBS interactions (out of 903 interactions identified in both two ChIA-PET experimental replicates) with both ends restricted within the same topological domain were selected as foreground training set. Equal number of candidate ERBS pairs that did not overlap with any ERα interactions identified in either experimental replicate was randomly selected background training set. We still used a logistic classifier to perform the classification, i.e.

in which is the indicator whether ERBS pair (i₁,i₂) formed an interaction, P is a probability function and are the model parameters.

Performance of the classifier was evaluated by a 5-fold cross-validation; AUC (area under curve) and ROC (receiver operating characteristic) were plotted as the average of each 5-fold cross-validation.

To predict novel ERBS–ERBS interactions, we trained the model parameters on the whole training set and applied it to the remaining candidate ERBS pairs. Predicted interactions that were not reported in either replicate were further analyzed.

laERBSs are more likely to be nucleosome depleted

To identify genetic or epigenetic features that could discriminate ERBSs involved in loop formation from those solo ones, we carried out integrative analysis of multiple (TF and histone modification) ChIP-seq data sets from the same cell line (MCF7). Firstly, we analyzed histone modification patterns around laERBSs and non-loop-associated ERBSs (nlaERBSs). Here laERBSs were selected as ERBSs overlapping with PETs 2+ clusters identified in both two ChIA-PET experimental replicates (903 clusters). To align histone modification signals around ERBSs, we selected the strongest summit of ERBS as the center if more than one ERBS was contained in the same end of an interaction cluster, resulting 1203 ERBSs as the foreground set. We randomly selected an equal number of ERBSs that were not associated with any interaction clusters identified in either ChIA-PET experimental replicate as the background set. It was noticed that laERBSs showed subtle depletion of histone modification signal in the center compared with nlaERBSs (Supplementary Figure S3A–C). laERBSs and nlaERBSs were further divided into two groups by the presence or absence of an ERE. Surprisingly, we found that a significantly more proportion of laERBSs with an ERE shown depleted histone modification signals in the center (Figure 2A–C, Supplementary Figure S4 and Supplementary Table S1) than the other three groups. DNase-seq data confirmed that laERBSs were more frequently located in open chromatin regions (Supplementary Figure S3D), especially those that contained an ERE (Figure ​(Figure2D2D and Supplementary Figure S5). Since only about 10% laERBSs located in promoter regions, our observation could not be explained by nucleosome depletion in active promoter regions (Supplementary Figure S6). Therefore, although to the best of our knowledge, there were no direct genome-wide nucleosome positioning data available in E2-induced MCF7 cells, our integrative analysis of multiple histone modifications and DNase-seq data suggested that nucleosomes were more likely to be depleted in laERBS. In fact, a previous paper (22) reported that ERα binding to DNA was not sensitive to nucleosomes; while on the contrary, we found that those laERBSs with an ERE showed significant nucleosome depletion in the center. This indicated that there might be some other co-factors that can bind to these laERBSs cooperatively to facilitate long-range interactions and be responsible for nucleosome eviction.

Open in a separate window

Figure 2.

Histone modification and DNase-seq signal profile for four groups of ERBSs. Average H3K4me1 (A), H3K4me3 (B) and H3K9ac (C) read-coverage against input surrounding laERBSs with and without an ERE and nlaERBSs with and without an ERE. Average DNase-seq tag counts surrounding laERBSs with and without an ERE and nlaERBSs with and without an ERE (D).

To test this hypothesis, we conducted an enrichment analysis against the ChIP-seq data of three known ERα's co-factors, and the general co-activator p300. More than 80% laERBSs (either with an ERE or not) overlapped with FoxA1, GATA3 and AP2γ binding sites (Figure ​(Figure3).3). The percentages were significantly higher than those of nlaERBSs (P < 2.2e−16; Supplementary Table S2), indicating that these three TFs are likely in the complex mediating ERα-associated long-range interactions. Reported RNAi knock-down experiments also support the notion that these three co-factors are important for ERα loop formation (23,24). Thus they might be associated with nucleosome eviction in laERBSs. For p300, the percentage of overlapping with laERBSs was not as high as that of those three co-factors. However, interestingly, p300 was significantly enriched in laERBSs with an ERE in comparison with laERBSs without an ERE (P = 7.59e−4), but showed no difference between nlaERBSs with an ERE and nlaERBSs without an ERE (P = 0.72), while none of the three co-factors have a similar pattern. A previous study (36) showed that p300 is a component of an estrogen receptor co-activator complex and the formation of the p300 complex is associated with nucleosome eviction (37). This indicated that p300 might be more frequently recruited to laERBSs with an ERE for higher level of nucleosome depletion. In addition, with the GRO-seq data (GEO accession numbers GSM678539 and GSM678540) we found that laERBSs showed stronger bidirectional transcription of small RNAs (smRNAs) compared to nlaERBSs (Supplementary Figure S7), with the highest transcription level occurring in laERBSs with an ERE (Supplementary Figure S7A). This is consistent with the notion that the bidirectional smRNAs transcription may help maintain the open chromatin structure in these regions (38). Taken together, these observations suggest that laERBSs, especially those with canonical EREs, are often associated with an open chromatin structure, which is likely resulted by the formation of the looping-specific protein complex of ERα, its co-factors and p300, and associated with the bidirectional transcription of smRNAs.

Open in a separate window

Figure 3.

Proportion of FoxA1, GATA3, AP2γ and p300 binding peaks that overlapped with four groups of ERBSs. P-values smaller than 0.05 were listed above the bars, which were given by one-tailed Fisher test in R. ***P < 2.2e−16.

Clustered ERBSs are more likely to be associated with loops

Of all the 1299 regions that formed PETs 2+ interaction clusters identified in both two ChIA-PET experimental replicates, 30% contained more than one ERBS. Statistical tests showed the clustered ERBSs within 3 kb were more likely to associate with DNA loops (P < 2.2e−16, 45% versus 20%).

A logistic classifier can predict laERBSs and associated long-range interactions

We tested the prediction power of the above features extracted from previous sections by building a classifier to predict ERα-associated long-range interactions based on TF and histone modification ChIP-seq data. Given a pair of ERα ChIP-seq binding peaks, a powerful classifier should be able to judge how likely they are actually forming an interacting DNA loop. Our predictor was built in two steps. First, we used the above features to distinguish laERBSs and nlaERBSs. We used 903 PETs 2+ clusters confirmed in both two replicates as the positive training samples, since they were more reliable than those non-reproducible ones. We evaluated the performance of the classifier with 5-fold cross-validation and got the average true positive rate (TPR) at 74% and average false positive rate (FPR) at 21% at the Bayesian threshold (0.5) with only three features, i.e. ERα ChIP-seq signals, neighboring ERBS distances and DNase-seq signals. AUC for predictor with single feature and ROC comparisons for different feature combinations were shown in Figure ​Figure4A4A and ​andC,C, respectively. Although transcriptional rate of bidirectional smRNA and the binding intensity of p300 showed good predictive power, they are redundant with ERα ChIP-seq signals, neighboring ERBS distances and DNase-seq signals. And unexpectedly, histone modification features did not provide additional information over DNase-seq signals for the prediction. We finally set the Bayesian threshold 0.2 to contain as many laERBSs (in the training set) as possible (97%) while the relatively high FPR (61%) can be further dealt with in the next filtering step.

Open in a separate window

Figure 4.

AUC and ROC curves for the predictors. (A) AUC for each feature to predict laERBSs. Here AUC was computed as the average area under ROC curve for predictors with each single feature for 5-fold cross-validation. (B) AUC for each feature to predict interactions between predicted laERBSs. AUC was computed in the same manner as in (A). "dis" is distance between two ERBSs and "idis" is the inverse distance. (C) ROC trained for ERα alone, ERα+distance of neighboring ERBSs, ERα+distance of neighboring ERBSs+DNase-seq and all the features. (D) ROC trained for ERα alone, ERα+distance and ERα+distance+6 significant features and all the features.

In the second step, we tried to predict long-range interactions between the putative laERBSs obtained from the first step. Recently, Hi-C data analysis suggests that chromatins are organized as large topological domains and interactions between regulatory elements are largely restricted within these domains (31). As shown in (31), although strengths of chromatin interactions between and within these domains may vary across different cell lines, such domains are roughly invariant across different cell types. Thus, we further restricted the two ends of the same interaction cluster within the same Hi-C-determined topological domain. This criterion only excluded 69 (~8%) PETs clusters from our positive set, but reduced more than one half of ERBS pair-wise combinations, a large percentage of which are likely to be noise (there are 97 202 candidate ERBS pairs by using constrain of topological domain while there exist 239 016 and 6 377 706 possible intra-chromosome ERBS pairs with and without restriction of 1-Mb genomic distance, respectively). With the above restrictions, we revealed 800 PETs 2+ interactions between filtered laERBSs within the topological domains, which served as the foreground training set. Classifiers with different feature combinations on the training set were evaluated by 5-fold cross-validation. The best performer achieved a TPR of 93% and an FPR of 8% at the Bayesian threshold 0.5. The intensity of ERBSs and the distance between candidate pairs acted as the top two significant features to distinguish foreground and background (Figure ​(Figure4B).4B). The binding intensity of FoxA1 and AP2γ had a similar predictive power, but was somewhat redundant with that of ERα. Beyond these features, DNase-seq and histone modification signals that were able to describe chromatin accessibility could also improve the performance. Finally, eight features (Supplementary Tables S3 and S4) that showed significant improvement over the ERα+distance combination (Figure ​(Figure4D)4D) were selected to build the final classifier.

We applied the classifier to predict interactions among the 97 202 candidate ERBS pairs. Originally, Fullwood et al. (15) reported 3527 high confidential interactions that share ERα ChIP-seq binding peaks at both ends for the combined two replicates. Among them, 3113 were restricted in the topological domains, 76% (2356) of which could be recalled, while only 9% (8805) of the remaining 94 089 unreported pairs were predicted to be ERα-associated interaction clusters. Therefore we conclude that our classifier is highly effective in predicting known ChIA-PET interactions from multiple ChIP-seq data.

Many novo predicted interactions are supported by other data sources

Next, we compared our newly predicted 8805 interactions with other sources of data to confirm their biological significance. They are significantly overlapping with known Pol2 PETs clusters (33) (P < 2.2e−16) (Figure ​(Figure5A5A and Supplementary Table S5). And the overlapping proportion was comparable with that of 3113 high-confidential ERBS–ERBS pairs reported in (15) (Figure ​(Figure5A5A and Supplementary Table S5). This observation suggested that our predicted novel ERBS interactions might contain a considerate proportion of true positives that were not captured by the published ChIA-PET analysis (15), probably due to the fact that the original ERα experiments were not saturated and many interactions were filtered out to reduce false positives.

Open in a separate window

Figure 5.

Evaluation of novo predicted interactions by other independent data sets. (A) Proportion of reported, predicted and unreported and rest of the candidate (non-predicted and unreported) ERBS–ERBS pairs that overlapped with Pol2 PETs 2+ and PETs 3+ clusters, respectively. (B) Proportion of reported, predicted and unreported and rest of the candidate (non-predicted and unreported) ERBS–promoter pairs that overlapped with Pol2 PETs 2+ and PETs 3+ clusters, respectively.

We thank Dr Monica Sleumer for helpful discussion and checking the proofreading.