Computational inference of H3K4me3 and H3K27ac domain length

Domain length prediction using Random Forest regression algorithm

Feature extraction: For each histone modification domain and data feature being studied, we used the fraction of overlapping basepairs (bps) in our models to eliminate the bias that could be introduced by differences in domain lengths. For example, if an H3K4me3 domain has length 100 base pairs and 15 of these bases also have CHD1 binding, fraction of overlap for that domain for the CHD1 feature is (15∕100 = 0.15). Note that a single H3k4me3 domain can simultaneously overlap with multiple instances of the same feature (e.g., multiple CHD1 binding sites within a domain). For these instances, we used the sum of all of such overlaps. For example, as depicted in the Fig. 1C domain B overlaps with two instances of domains A, therefore the final score for this domain for feature A is = sum(bp_A)/length(B).

Model building: Random Forest regression with 100 trees was performed on the constructed data sets. We employed these models to unfold the relationship between domain length and the constructed features. We used an implementation available in scikit-learn library (Pedregosa et al., 2011). Gini impurity was used as a quality measure for choosing best splits. At each level √m randomly chosen attributes were considered as candidates for split (where m is the number of all attributes).

Feature ranking: Feature rankings enable identifying and ranking important features for the regression problem at hand. When building a decision tree, the feature that leads to the greatest decrease in Gini impurity score is chosen at each split. Importance score for each feature is the mean decrease of impurity for all tree nodes, which is averaged over all trees in the ensemble. Calculating importance scores enable comparison of attributes, but does not state the significance of importance scores. To assess significance we employed a Monte-Carlo technique based on contrast attributes, which are random permutations of the original attributes. To quantify the significance of a feature importance score, we followed the below procedure: This procedure enabled us to identify a small set of significant features for our models.

1. For each of the original attributes its values were permuted at random and added as a new contrast attribute to the original dataset.

2. Random Forest model was trained on the dataset consisting of both original and contrast attributes.

3. Standard deviation (SD) of importance scores for contrast attributes was calculated. The value equal to 2 × SD was used as a cutoff for minimal significant difference between importance scores.

4. Feature importance scores were sorted in decreasing order. Differences between subsequent scores were calculated. We looked for the last pair of subsequent scores for which the difference was larger than the cutoff value. The larger of the scores from that pair constituted a threshold for dividing features into significant and not significant.

H3K4me3- and H3K27ac-domain lengths are informative for cell-type specificity

It has been observed that epigenetic marks can decorate loci with domains of varying length. For example, when we study the domain lengths for H3K4me3, we observe that the length distribution is not uniform and the length of these peaks vary between a few hundred and 20 kb as shown in the length distribution of H3K4me3 domains in H1 human embryonic stem cells (hESCs) (Fig. 1A). The span of epigenetic mark deposition has gained attention in recent years with several studies, including ours, showing that particularly long stretches of DNA marked with H3K4me3 or H3K27ac specifically coincide with cell-type specific promoters or enhancers, respectively (Hnisz et al., 2013; Chapuy et al., 2013; Parker et al., 2013; Benayoun et al., 2014; Bernstein, Meissner & Lander, 2007).

Given the association of large domains with functionally important DNA elements, we built a computational framework that can (i) assess whether—and to what extent—the length of an H3K4me3 (or H3K27ac) domain at a locus can be predicted from other genomic and epigenomic characteristics of that locus and (ii) quantify the ability of each genomic and epigenomic characteristic to predict the length of domains of these two marks. Our framework consisted of three stages (summarized in Fig. 1B): (1) extracting genomic and epigenomic characteristics of each domain; (2) building a regression model (based on Random Forest regression) for the length of the domains as a function of other data characteristics; and (3) delineating the predictive genomic signature for domain length by prioritizing and selecting important predictors. The signatures we obtained will help us to identify candidate molecular mechanisms for setting or maintaining domain length of this epigenetic mark.

H3K4me3 and H3K27ac deposition lengths can be predicted with high accuracy by integrating other genomic datasets

For each H3K4me3 or H3K27ac domain, we extracted other genomic and epigenomic characteristics from ENCODE datasets. We used the fraction of overlap (as explained in methods) to normalize feature overlaps based on domain length, so longer domains would not have an inflation of overlapping features. Next, in seven ENCODE cell lines (H1 and H9 hESCs, K562, HCT116, GM12878, A549, HeLa) as well as in human pancreatic islets, we built Random Forest Regression models to predict the length of H3K4me3 and H3K27ac domains. Our analyses showed that these models can predict the length of epigenetic marks with high accuracy (>0.9 correlation between observed and predicted lengths; Fig. 2).

Conserved genomic features are predictive for H3K4me3 and H3K27ac deposition lengths

Random Forest quantifies the importance of each feature by calculating its contribution to the overall classification accuracy (see methods for details). We calculated feature importance scores for each genomic and epigenomic data type we used in our models. Next, we assessed the significance of these scores by comparing them with scores produced by random permutations of the original features (i.e., contrast features). With this methodology, we identified the significant predictors of length for H3K4me3 and H3K27ac in each studied cell type. For example, we showed that in K562 CHD1 is one of the most important predictors for both H3K4me3 and H3K27ac domains (Figs. 3 and 4).

Our random forest regression models revealed that there are frequently repeating features that are predictive for both H3K4me3 and H3K27ac epigenetic marks, which point out to putative upstream or downstream associations between the length of epigenetic marks and these predictive features. Certain genomic and epigenomic datasets frequently rank as important for predicting the length of these two marks across cell types (Figs. 3C–3D). For H3K4me3, CHD1- and Pol2-binding and the H3K79me2 histone modification are the most important predictors as assessed by feature importance scores. For H3K27ac, TEAD1- and CHD1-binding and the presence of other acetylation marks (H3K4ac and H2AK5ac) are important genomic characteristics.

Elongated CHD1 binding might promote robust and increased expression patterns

Among the top predictors, CHD1 is an H3K4me3-associated chromatin remodeler and it is shown to be essential for the expression of key developmental genes in mouse embryonic stem (ES) cells as well as for reprogramming of fibroblasts to the pluripotent stem cells (Gaspar-Maia et al., 2009). Moreover, CHD1 binding is shown to be essential for the assembly of preinitiation complexes (PIC), which is necessary to start transcription. PIC includes co-activators like Mediator complex and the TFIID complex, along with Pol2 and the general transcription factors (Lin et al., 2011). Given the importance of CHD1 in regulating PIC assembly and ultimately gene expression, we have conducted further analyses to study the relationship between CHD1 binding and the length of H3K27ac and H3K4me3 domains. Our analyses showed that there is a positive correlation between the length of CHD1 and the lengths of both H3K27ac and H3K4me3 domains (Fig. 4), implying that elongated CHD1 binding takes place at regions where longer H3K4me3 and longer H3K27ac deposition is observed. Broad H3K4me3 domains have been associated with increased and robust expression of genes, similarly genes marked by longer H3K27ac domains (i.e., super/stretch enhancers) have been associated with increased expression. Our findings show that longer H3K4me3 and H3K27ac domains have longer deposition of CHD1 binding; which might imply a more robust assembly of mediator complex and PIC enabled by the elongated binding of CHD1 protein. This robust assembly of PIC at these loci might explain their robust and increased expression patterns.

Broader H3K4me3 and H3K27ac are bounded by CTCF mediated interactions

CTCF, is a TF that is important for the insulator activity. Insulators are boundary elements that restrict the regulatory interactions between enhancers and promoters within a boundary element. An intriguing hypothesis based on the high predictive value of CTCF for the length of H3K4me3 and H3K27ac domains is that the longer H3K27ac and H3K4me3 domains tend to be enriched at the boundary of chromatin domains marked by the CTCF binding. To test this, we analyzed CTCF mediated chromatin interactions obtained by the Chromatin Interaction Analysis by Paired-End Tag Sequencing (ChIA-PET) technology in GM12878 cell line (Tang et al., 2015). We split H3K4me3 and H3K27ac domains into five bins based on their domain length and for each of these bins, we studied the localization of CTCF mediated chromatin interactions with respect to the center of the domain overlapping the CTCF fragment. Our analyses revealed an interesting pattern, where we observed that with increasing length of H3K27ac and H3K4me3 domains, CTCF interactions tend to cluster at the beginning or at the end of the domains (Fig. 5), which implies that longer H3K4me3 and H3K27ac domains tend to be at the boundary of CTCF defined boundary domains.