Identification of promoters and enhancers in ESCs

We performed ChIP-seq experiments of H3K4me3, H3K27me3, H3K27ac, and H3K4me1 in mouse ESCs to identify the different types of regulatory regions. This set of histone modifications has been previously used to distinguish between AEs, PEs, APs, and BPs in mouse ESCs [21,22]. We next generated a 9-state chromatin segmentation model of ESCs using ChIP-seq data (Fig 1A and 1B). As expected, active states mark transcriptionally active regions, while repressed states denote transcriptionally repressed regions (Fig 1B; see also our previously published RNA-seq data [23]).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Identification of repressed and active functional regions in ESCs.

(A) State definition of the chromatin segmentation model in ESCs. The values represent the probability (from 0 to 1) of finding each histone modification (vertical) in genomic segments of states 1 to 9 (horizontal). The cells of the matrix are colored according to the value of probability they contain inside. Red: states with histone modifications associated to activation (active, 1–4); dark yellow, H3K4me1-only state (Intermediate, 5); grey, states in which H3K27me3 was present (repressed, 6–8); dark grey, poised states, in which H3K27me3 colocalized with H3K4me3 and/or H3K4me1 (states 6 and 7); light grey, H3K27me3-only regions (state 8); and white, unmarked state (9). (B) Example of a genomic region containing two expressed genes (Skap2 and Halr1), which are covered by active states (in red), and a cluster of repressed genes (HoxA), which are covered by repressed states (in grey). Active chromatin segments integrate the signal of H3K27ac, H3K4me3, and H3K4me1 and lack H3K27me3. Repressed chromatin segments integrate the signal of H3K27me3, H3K4me3, and H3K4me1 and lack H3K27ac. Expression of Skap2 and Halr1, and silencing of HoxA genes, were confirmed by the RNA-seq profiles [23]. Y-axis represents normalized count of reads by total reads. The screenshot was taken from the UCSC Genome Browser [62]. (C) Enrichment of state transitions (e.g., number of observed transitions divided by the number of expected transitions by chance) from the segments of one state (vertical) towards the segments of another state (horizontal) in the linear chromatin. The cells of the matrix are colored according to the value of enrichment they contain inside. (D) Expression of genes associated to active promoters (AP; 10,786 genes) or bivalent promoters (BP; 3,459 genes). The dotted line represents 1 FPKM. (E) Top GO biological process (2018 categories) for each list of genes in D.

https://doi.org/10.1371/journal.pcbi.1009368.g001

To understand the resulting map of states, we calculated the matrix of transition enrichments between all states in the model. The transition value between two different states, x and y, is defined as the number of times that a segment of state y is found after a segment of state x, as measured from left to right in the linear genome. The enrichment score is defined as the ratio between the observed number of transitions and the expected number of transitions by chance. Two groups of states (active, 1–4, and repressed, 6–8) clearly emerged from the global picture (Fig 1C) (note that state 5 represents a distinct category; see Discussion). The high enrichment of transitions between states belonging to the same category suggests that they might mark the same functional regulatory regions, rather than be caused by different functional regions separated by the unmarked stated. Visual inspection confirmed that states 1 to 4 marked the same active regulatory regions, revealing that differences in the state definition are due to differences in the shape of the ChIP-seq peaks (S1A Fig). We also observed that the repressed states 6 and 7 decorated poised or bivalent regulatory regions (e.g., marked with H3K27me3, in combination with H3K4me3 and/or H3K4me1), whereas state 8 was generated by the tail-end of broad peaks of H3K27me3 (S1B Fig). Similarly, state 5 was associated with the tail-end of broad peaks of H3K4me1 in active regions, but it also associated with single peaks of H3K4me1 near active regions (S1A Fig).

Based on these results, we decided to merge the contiguous segments of states 1 to 4 as a list of potential active regulatory regions, and the segments of states 6 and 7 as a list of potential poised or bivalent regulatory regions. We reasoned that functional regions should have a minimum length of 600 bp and thus discarded shorter regions, which we considered as background signal. We also discarded those cases in which an active region and a poised region were contiguous, as this was ambiguous. Promoters were defined as those regions that overlapped by at least 1 bp to a region ± 500 bp around a TSS according to RefSeq [24]. Enhancers were defined as regions that were not classified as promoters and overlapped by at least 1 bp with a peak of the enhancer mark p300 [21]. As H3K4me3 can be present in enhancers [25–28], we did not discard enhancers containing H3K4me3, although this histone modification has been traditionally only associated with promoters. In total, we found 9,421 APs, 3,344 BPs, 16,904 AEs, and 2,699 PEs (S1–S4 Tables).

Next, we matched our set of promoters to their target genes, using the same parameters as before (overlap by at least 1 bp to a region ± 500 bp around a TSS according to RefSeq [24]). Using RNA-seq data [23], we confirmed that genes associated with APs are expressed, while genes associated with BPs are not (Fig 1D). Gene ontology (GO) term enrichment analysis performed with Enrichr [29] confirmed that genes with APs are involved in housekeeping roles, while genes with BPs are mostly related to development and differentiation (Fig 1E), as is expected in mouse ESCs. On the other hand, we used available high-throughput chromosome conformation capture (3C) data of Hi-C [30], to link AEs and PEs with target genes. As interacting enhancers and promoters have been shown to match their chromatin state [31], we associated an enhancer with the target gene of a promoter when both enhancer and promoter are in the same category (e.g., both active, or both repressed) and each one overlaps with one of the two sides of the same Hi-C significant interaction (total of 43,892,155 significant Hi-C interactions). We confirmed that PEs were significantly enriched in interactions with BPs over APs (p < 2.2e-16 Exact Binomial Test, observed probability: 0.29, expected probability by chance: 0.24), whereas AEs were significantly enriched in interactions with APs over BPs (p < 2.2e-16 Exact Binomial Test, observed probability: 0.80, expected probability by chance: 0.76). In total, we found 10,786 genes associated to APs, 3,459 genes to BPs, 10,206 genes to AEs, and 2,526 genes to PEs (S1, S2, S5 and S6 Tables). Likewise, 15,841 AEs and 2,466 PEs were associated to at least one gene, and 8,931 APs and 2,443 BPs, to at least one enhancer (S5 and S6 Tables).

Development of a predictive model of gene expression using histone modifications at enhancers

After identifying the set of enhancers and promoters in ESCs, we built the gene expression predictive models. We first performed additional ChIP-seq experiments for other histone modifications, in order to have additional variables to predict gene expression. Specifically, we performed ChIP-seq experiments for trimethylation of histone H3 at lysine 36 (H3K36me3), ubiquitination of histone H2B (H2Bub), monomethylation of histone H3 at lysine 27 (H3K27me1), dimethylation of histone H3 at lysine 27 (H3K27me2), trimethylation of histone H4 at lysine 20 (H4K20me3), and dimethylation of histone H3 at lysine 79 (H3K79me2). The input of our predictive models consisted of the ChIP-seq data of these six histone marks as well as the four histone marks previously used to define promoters and enhancers (H3K4me3, H3K4me1, H3K27ac and H3K27me3), together with our previously-published RNA-seq expression data [23].

Initial studies on gene expression prediction revealed that using two to three histone modifications (rather than larger set) are sufficient to accurately predict gene expression, and do not find substantial improvements with the addition of other histone modifications into the models [9–11]. Indeed, follow up publications directly utilize three to four histone modifications [12, 13]. However, our objective is not only to predict gene expression, but also to assess histone modification contribution in enhancer predictive models in comparison to promoter predictive models. Therefore, we used a set of ten histone modifications, each one with different properties (associated to activation, associated to repression, broad marks, sharp marks, etc.). A total of 11,387 protein coding genes previously associated to a promoter (active or bivalent) and at least one enhancer (active or poised) entered the modelling. We divided the set of genes into two subsets: training and test. The Pearson’s correlation coefficient (r) between the measured expression in the test subset and the predicted one was used to assess performance.

As one of our aims is to compare enhancer predictive models to promoter predictive models, we needed to build them using the same approach to identify enhancers and promoters so they are comparable. Therefore, we first generated a predictive model for the promoters identified using our approach (named Hi-C–all promoter model) to confirm that histone modifications at these elements are predictive of gene expression. The predictive capacity of promoters has been previously shown in several cell types from different model organisms, including ESCs, where promoters were defined as a pre-set distance from a TSS [9–15]. As a control, we repeated the predictive model learning in a training subset in which expression values were randomized. We obtained an r-value of 0.81 for the promoter model, and an r-value of –0.07 for the random promoter model (S2A Fig). The low performance of the random promoter model strongly indicated that the high predictive power of the Hi-C–all promoter model was not due to random structures in the data. Importantly, the performance of our promoter model was comparable to previously described predictive models, in which r-values around 0.8 were reported [9–15]. Indeed, Karlic and colleagues obtained an r-value of 0.77 in CD4+ T-cells in the seminal paper on gene expression prediction from histone modification levels [9]. Coefficients and p-values of the predictors for all the predictive models generated in this study can be found in S7 Table.

Next, we trained a second predictive model of gene expression (the Hi-C–all enhancer model) using the levels of histone modifications at the previously-defined enhancers as predictors. We obtained a performance in the test subset of r = 0.38 (Fig 2B). Although the r-value of this model is still modest, this is the first time to our knowledge that enhancers have been shown to be predictors of gene expression through their histone modification levels. Critically, when the expression data for learning the model were randomized, the performance was poor (r = –0.15, Fig 2B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Performance and variable importance of enhancer and promoter Hi-C-top predictive models in ESCs.

Predicted expression of the test subset of genes calculated by the models versus their measured expression by RNA-seq. Model performances are represented by the Pearson’s correlation (r) between predicted and measured expression values. (A) Left, the model trained on the promoter regions associated to at least one enhancer using the top significant interactions of Hi-C (Hi-C–top promoter model). Right, the performance of the same model after randomizing the expression of the training subset of genes. The color bar represents the density of dots. (B) Left, the model trained on the enhancer regions associated to at least one promoter using the top significant interactions of Hi-C (Hi-C–top enhancer model). Right, the performance of the same model after randomizing the expression of the training subset of genes. The color bar represents the density of dots. (C) Importance of each histone modification used to train the Hi-C–top promoter predictive model. Importance is defined as the contribution of each variable in the linear regression predictive model and corresponds to the absolute value of the t-statistic for each model parameter. (D) As for C, but for the Hi-C–top enhancer predictive model. (E) As for B, but the model is trained without H3K27me3 as predictive variable. (F) As for D, but the model is trained without H3K27me3 as predictive variable.

https://doi.org/10.1371/journal.pcbi.1009368.g002

We hypothesized that the modest performance of the Hi-C–all enhancer model could be due to some enhancer–promoter associations that were simply regions in close 3D proximity but not functionally linked. To enrich the set of interactions for functional promoter–enhancer loops, we applied a more restrictive threshold on the Hi-C significant interactions (FDR = 0 and ln(p-value) ≤ –100). We obtained a total of 5,555,844 interactions (8% of the total interactions). We then recalculated the enhancer-promoter-gene associations and obtained 1,846 PEs associated to 1,382 BPs and to 1,434 target genes, and 11,777 AEs associated to 7,211 APs and to 8,254 target genes (S8 and S9 Tables). We selected the protein-coding genes included in the new associations (a total of 8,639 genes) to recalculate the predictive models (hereon in termed Hi-C–top models) and random models for promoters (Fig 2A) and enhancers (Fig 2B). While the Hi-C–top promoter model performed similarly to the previous model (r = 0.79 vs. r = 0.81, respectively), the Hi-C–top enhancer model was significantly improved (r = 0.49 vs. r = 0.38). This result further confirmed that enhancers, as well as promoters, possess a quantitative relationship with gene expression.

We now know that an enhancer preferentially interacts with promoters located in the same topologically associating domain (TAD) rather than those located in neighboring domains [32]. Moreover, TADs have an average size of around 1 Mb. Therefore, assigning genes to enhancers located in a distance lower than 1 Mb might seem appropriate. Indeed, when evaluating a new ESC predictive model that associates enhancers to promoters of the same chromatin state that are closer than 1 Mb (1 Mb model; 11,986 protein-coding genes), we achieved a performance of r = 0.34 (S2C Fig). Nonetheless, that performance is lower than the ESC models based on Hi-C data that we have built previously (r = 0.38 and r = 0.49 for Hi-C–all and Hi-C–top enhancer models, respectively). This suggests that matching genes to regulatory elements by 1 Mb distance leads to some false-positive associations, yet maintaining its predictive capacity.

Finally, from the Hi-C–top interactions we selected those involving a distal enhancer (> 5 Kb from a TSS, a total of 5,235 AEs and 696 PEs) to confirm that the predictive capacity was not exclusive of proximal enhancers. We generated a new distal enhancer model (Hi-C–top_distal; 7,925 protein-coding genes) that properly predicted gene expression with a r = 0.41 (S2D Fig).

Enhancers and promoters exhibit similar histone modification patterns. We therefore wondered whether the histone modifications mostly contributing to the prediction of expression were the same ones as well, or whether there were differences in the contributions between the enhancer and the promoter predictive models. To address this, we assessed variable importance in the Hi-C–top model for promoters and enhancers. Notably, H3K27me3—a histone modification associated with transcriptional gene repression—was the prevalent mark in both classes of regulatory elements (Fig 2C and 2D). In contrast to promoters, in which H3K27me3 has a relatively similar importance as other marks (e.g., H2Bub, H3K4me3, and H3K36me3, Fig 2C), H3K27me3 in enhancers represented up to 55% of the total importance (Fig 2D). Therefore, even though promoters and enhancers contribute to predict gene expression mostly through H3K27me3, this contribution seems to be uniquely driven by H3K27me3 in the enhancer predictive model and shared by other histone marks in the promoter predictive model.

Interestingly, H3K27ac—considered the canonical marker of enhancer activation [33]—had little importance in the enhancer model. H3K27me3 and H3K27ac have antagonistic effects and are generally mutually exclusive marks since they occur on the same lysine and are chemically prohibited [34]. Thus, we reasoned that H3K27ac importance was masked by H3K27me3 presence in the predictive model, as both histone modifications are not independent variables. Indeed, when excluding H3K27me3 as a predictive variable, the enhancer model maintains its predictive capacity (r = 0.36; Fig 2E), and the most important variable is in fact H3K27ac (Fig 2F). However, in this case, H3K27ac has a relatively similar importance as H3K27me2. Therefore, H3K27me3 seems more predictive than H3K27ac at enhancers in relationship with the rest of the histone modifications.

LOESS normalization of ChIP-seq and RNA-seq data from heterogeneous sources

We next wanted to determine whether enhancers are predictive of gene expression in other cellular contexts besides ESCs, and whether a predictive model learned in one cell type could predict gene expression in another. As true colocalization of H3K27me3 with H3K4me1 or p300 in the same DNA fragment has been only studied in ESCs, it is still not clear whether PEs exist in other developmental scenarios [20]. To address this issue, we took advantage of the capacity of PEs and BPs to switch into an active state for certain cell types, in a lineage-specific manner during differentiation from ESCs, while remaining inactive in others [20]. We hypothesized that differentiation data could be used to obtain predictive models exclusively from PEs and BPs. We focused on several time points for two cell differentiation mouse models: i) cardiac lineage: mesoderm, cardio precursors, and cardiomyocytes [35]; and ii) neural lineage: neural precursors and cortical neurons [30].

We downloaded RNA-seq and ChIP-seq data of five histone modifications (H3K27me3, H3K4me3, H3K27ac, H3K4me1, and H3K36me3) that were available in the literature for both differentiation models. To remove potential biases (e.g., due to the source of data generation or to a batch effect), we normalized the sequencing samples of the same feature at all the available time points. For this, we applied a normalization based on a local regression (LOESS) that was originally proposed for the pairwise normalization of expression microarrays [36] but generalized for multiple arrays [37]. LOESS normalization is based on a MA methodology, where M is the log₂ ratio of the intensities of the samples, and A is the log₂ of the average intensity. It assumes that the intensities of the two samples should be equal, therefore M = 0. Finally, corrections based on a LOESS are applied to obtain a MA plot in where the regression line approximates M = 0.

We first applied this normalization method over the expression of a set of 20,706 protein-coding genes in the mouse genome, at each differentiation time point, and ESCs (MA plots of each differentiation time point against ESCs before and after LOESS normalization are shown in S3A Fig). As LOESS normalization assumes that expression is equal in all samples, a general balance in global expression distribution of all time points is expected (S3B Fig). To further confirm the normalization efficiency, we tested its performance on two different subsets of genes: housekeeping genes and bivalent genes. We hypothesized that housekeeping genes would show a balanced distribution of expression after normalization, while bivalent genes would increase their expression globally during differentiation (as some of them are activated). For this, we extracted a list of mouse housekeeping genes across 14 mouse tissues from the literature [38] to check their expression. We also evaluated the normalization on our list of bivalent genes (e.g., those associated to BPs). Indeed, after LOESS normalization, the expression of housekeeping genes was correctly balanced (Fig 3A), whereas the expression of bivalent genes maintained the characteristic pattern of increased expression across time (Fig 3B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. RNA-seq and ChIP-seq data before and after LOESS normalization.

(A) Raw and normalized expression of 3,277 housekeeping genes along cardiac and neural differentiation from ESCs. (B) Raw and normalized expression of 3,459 bivalent genes along the same time points as A. (C) Raw and normalized H3K4me3 ChIP-seq signal levels at all 2-Kb bins of the genome. (D) Raw and normalized H3K4me3 ChIP-seq signal levels at 3,344 BPs. (E, F), same as C and D, respectively, but for H3K27me3. CM, cardiomyocytes; CN, cortical neurons; CP, cardio precursors; MES, mesoderm; NPC, neural precursors.

https://doi.org/10.1371/journal.pcbi.1009368.g003

We then ran the same normalization method on the ChIP-seq samples for H3K27me3, H3K4me3, H3K4me1, and H3K36me3 at the same time points (MA plots before and after LOESS normalization over the full set of bins of 2 Kb at Chr19 are shown in S4 Fig). Similar to expression data, we evaluated our normalization method for the H3K4me3 and H3K27me3 ChIP-seq levels across differentiation on two different sets of genomic regions: the whole collection of bins of 2 Kb in which the genome is segmented, and the coordinates of our collection of BPs. We hypothesized that global ChIP-seq levels of the whole genome would become balanced irrespectively of the particular histone modification analyzed, while BPs should present a different pattern for H3K4me3 and H3K27me3 (e.g., increase and decrease of signal along differentiation time points respectively, as a subset of the bivalent genes becomes activated during time). Indeed, after LOESS normalization, the levels of H3K4me3 along the whole genome were balanced (Fig 3C), whereas the same histone mark at BPs presented a clear pattern of increase during differentiation (Fig 3D). Notably, for H3K27me3, we observed the same balance in the levels along the whole genome after LOESS normalization (Fig 3E), while BPs exhibited a pattern of decreased signal across differentiation, in contrast to that observed for H3K4me3 (Fig 3F). In all cases, therefore, there is a substantial improvement after applying this normalization approach, while the analysis solely based on data before normalization would be misleading.

Poised enhancers and bivalent promoters are good predictors of gene expression during differentiation

After normalizing the data on expression and histone modifications across differentiation, we next generated predictive models using PEs and BPs for each differentiation time point. We used the Hi-C–top interactions involving PEs, BPs, and target genes in ESCs (1,846 PEs and 1,382 BPs associated with 1,434 target genes). From this dataset, a total of 1,063 protein-coding genes were used in the analysis. As the number of genes is smaller than in the previous gene sets, we decided to build the models at each cell type from the full set of genes to be evaluated in the rest of the differentiation time points. This approach has the advantage of allowing us to check whether the relationship between gene expression and histone marks in PEs is universal. This would be true if a model trained in a specific cellular context has a good performance in predicting gene expression in another one. We hypothesized that, as shown previously for promoters and gene bodies [9–11,14], there is a universal relationship between gene expression and histone modifications at PEs.

As a control, we randomized expression data and calculated predictive models for each time point. Next, we evaluated the performance of the randomized models on each differentiation dataset. We observed that predictive models for PEs and BPs obtained a significantly higher performance than randomized controls (Figs 4A, S5 and S6A). Surprisingly, all PE models achieved the best performance in cardiomyocytes (Table 1), suggesting that cardiomyocyte gene expression is easier to predict than the gene expression of the other time points. Moreover, all PE models had similar performances at each time point (Table 1). These observations are also true for the BP models (Table 2). Taken together, our results indicate that there is a universal quantitative relationship between gene expression and histone modifications at PEs and BPs across cardiac and neural differentiation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. PE models trained in differentiation time points.

(A) Performance of each differentiation enhancer model on the rest of the differentiation time points as compared to performance over random models. Performance is represented as Pearson’s correlation (r) between predicted expression and measured expression. Significance was assessed using a paired Student’s t-test between the performance of the models and the performance of the random models paired by the differentiation test set (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05). CM, cardiomyocytes; CN, cortical neurons; CP, cardio precursors; MES, mesoderm; NPC, neural precursors. (B) Importance of the histone modifications for each differentiation enhancer model. Importance is defined as the contribution of each variable in the linear regression predictive model and corresponds to the absolute value of the t-statistic for each model parameter.

https://doi.org/10.1371/journal.pcbi.1009368.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Performance of each PE differentiation model at every differentiation time point.

https://doi.org/10.1371/journal.pcbi.1009368.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Performance of each BP differentiation model at every differentiation time point.

https://doi.org/10.1371/journal.pcbi.1009368.t002

In order to confirm the predictive capacity of distal PEs (>5 Kb from a TSS, a total of 696 PEs), we generated new distal PE models (S7 Fig). Indeed, distal PEs maintained the predictive capacity. In this case, 486 protein-coding genes were included in the modelling.

Finally, we assessed the variable importance of the PE models for identifying differences in the contribution of each histone modification to the predictive models in different cellular contexts (Fig 4B). Strikingly, we observed that, in general, the two most important variables are H3K27me3 and H3K36me3. H3K36me3 was the most important histone modification for cardiac differentiation, whereas H3K27me3 was the most important for neural differentiation. In general, H3K27ac followed the above-mentioned histone modifications. H3K4me1 had a relatively low relevance to the predictive models, which suggests that it is involved in delimitating the enhancer regions rather than in contributing to its function. H3K4me3, which was vastly associated with promoter activity, is accordingly the least informative mark for the prediction of gene expression using PEs, suggesting that H3K4me3 is not associated to enhancer activity. The variable importance in the BPs showed that H3K27me3, H3K27ac and, importantly, H3K4me3 were the most informative variables (S6B Fig). Our results suggest that the quantitative relationship between histone modifications varies according to their location in PEs or BPs. Critically, even though there is a universal quantitative relationship between histone modifications and gene expression, this relationship can vary depending on the cellular context.

As for ESCs, we also assessed the effect of H3K27me3 absence over H3K27ac importance in the current scenario. Therefore, we generated new PE predictive models without H3K27me3 (S8A Fig). Interestingly, H3K27ac was generally the most predictive variable for these predictive models, followed by H3K36me3 (S8B Fig). Therefore, we wondered whether in absence of H3K27ac, H3K27me3 would be also more important than H3K36me3. Indeed, when generating new PE models in absence of H3K27ac (S9A Fig), H3K27me3 had the highest importance in all the cases (S9B Fig). Therefore, H3K27me3 seems the most informative variable for predicting gene expression from PEs.

The importance of H3K36me3 seen in Fig 4B agrees with the fact that almost 60% of the PEs that entered the modelling are located within gene bodies. As H3K36me3 is located in the gene body of active genes [39,40], intragenic enhancers also become marked when the genes start to be expressed during differentiation. Therefore, we divided PEs into two groups, intragenic or intergenic, and built new PE predictive models. Both, intragenic and intergenic models were capable of predicting gene expression (S10A and S11A Figs). When assessing for variable importance, we observed that, as expected, H3K36me3 maintained its high contribution in the intragenic predictive models (S10B Fig). However, H3K36me3 importance was reduced in the intergenic predictive models (S11B Fig). On the contrary, H3K27me3 maintained its importance in both, intergenic and intragenic models. Thus, H3K27me3, and not H3K36me3, behaves as a truly universal predictor of PE activity.

Poised enhancers and bivalent promoters are good predictors of gene expression in mouse embryonic tissues

In order to extend our findings from in vitro differentiation to in vivo, we learnt predictive PE and BP models from mouse developmental stages at different tissues. We downloaded ChIP-seq data of H3K27me3, H3K27ac, H3K36me3, H3K4me3 and H3K4me1 and RNA-seq on mouse embryos (heart tissue from 10.5 embryonic day, liver tissue from 11.5 embryonic day, neural tube tissue from 12.5 embryonic day, kidney tissue from 14.5 embryonic day, and lung tissue from 15.5 embryonic day) from ENCODE [41]. We first normalized the ChIP-seq and expression data following the LOESS approach. A total of 1,087 protein-coding genes entered the analysis. We observed that PEs were also predicting gene expression during mouse embryo development (Fig 5A). Again, when assessing for variable importance, we found that H3K27me3 was contributing the most, followed by H3K27ac and H3K36me3 (Fig 5B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. PE models trained using developmental stages.

(A) Performance of each differentiation PE model on the rest of the developmental stages as compared to the performance over the random models. Performance is represented as Pearson’s correlation (r) between predicted expression and measured expression. Significance was assessed using a paired Student’s t-test of the performance of the models or of the random models paired by a differentiation test set (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05). (B) Importance of histone modifications for each development PE model. Importance is defined as the contribution of each variable in the linear regression predictive model and corresponds to the absolute value of the t-statistics for each model parameter. Heart10.5, heart tissue from 10.5 embryonic day; Kidney14.5, kidney tissue from 14.5 embryonic day; Liver11.5, liver tissue from 11.5 embryonic day; Lung15.5, lung tissue from 15.5 embryonic day; NeuralTube12.5, neural tube tissue from 12.5 embryonic day.

https://doi.org/10.1371/journal.pcbi.1009368.g005

Next, we confirmed that BPs were also predicting gene expression during mouse embryo development (S12A Fig). In this case, the most predictive variable was H3K27ac followed by H3K27me3 and H3K4me3 (S12B Fig). Therefore, we further confirmed that differences in variable importance between PE models and BP models exist, which suggests that different histone modifications relate better to PE and BP function respectively.

Cell culture

E14Tg2A ESCs were cultured feeder-free on 15-cm plates coated with 0.1% gelatin. Plates were coated with gelatin for 15 min at 37°C, and then non-bound gelatin was removed. ESCs were cultured with Glasgow minimum essential medium (Sigma) supplemented with β-mercaptoethanol, sodium pyruvate, penicillin–streptomycin, non-essential amino acids, GlutaMAX, 20% fetal bovine serum (Hyclone), and leukemia inhibitory factor (LIF).

Chromatin immunoprecipitation

Cells were grown in 15-cm plates until 70% confluency and crosslinked in 1% formaldehyde in growth medium for 10 min at room temperature in a shaker. To stop fixation, glycine was added to a final concentration of 0.125 M and incubated for 5 min at room temperature. Cells were then washed twice with ice-cold PBS and harvested by gently scrapping plates (on ice) in PBS plus protease inhibitors. Cells from two 15-cm plates were pooled together and centrifuged at 3,400 × g at 4°C for 5 min. Cell pellets were frozen at –80°C until use.

Chromatin was prepared by resuspending the crosslinked pellet in 1.3 ml ice cold ChIP buffer [1× volume SDS buffer (100 mM NaCl, 50 mM Tris-HCl pH 8.1, 5 mM EDTA pH 8.0, and 0.5% SDS) and 0.5 × volume Triton dilution buffer (100 mM NaCl, 100 mM Tris-HCl pH 8.6, 5 mM EDTA pH 8.0, and 5% Triton X-100)] plus proteinase inhibitors. Samples were sonicated 40 cycles (30 seconds on/30 seconds off) in a Bioruptor Pico (Diagenode) and centrifuged at 16,000 × g at 4°C for 20 min to remove the cell debris. To check chromatin size, a 25-μl aliquot was mixed with 175 μl of PBS plus 5 μl of 20 mg/ml proteinase K, and de-crosslinked for 5 h at 65°C. DNA was purified using the QIAquick PCR purification kit (Qiagen), quantified in Nanodrop, and checked by electrophoresis on a 1.2% agarose gel.

ChIP experiments were performed using 30 μg of chromatin (DNA) and 5 μg of antibody in a final volume of 500 μl ChIP buffer. Aliquots of 5 μl were removed as input material (1%). ChIP samples were incubated overnight at 4°C on rotation, and then Protein A agarose beads (Diagenode) (42 μl per ChIP) were blocked 30 min with 0.05% BSA, washed, and added to the ChIP reaction. Samples were incubated for 2 h at 4°C with rotation. After incubation, beads were washed three times with 1 ml of low-salt buffer (140 mM NaCl, 50 mM HEPES pH 7.5, and 1% Triton X-100) and once with 1 ml high-salt buffer (500 mM NaCl, 50 mM HEPES pH 7.5, and 1% Triton X-100). ChIPed material was eluted from the beads in 200 μl freshly prepared elution buffer (1% SDS, 100 mM NaHCO3) at 65°C in a shaker (1000 rpm) for 1 h. Input samples were also brought to 200 μl with elution buffer. After addition of 8 μl of 5 M NaCl to the eluted chromatin and input samples, samples were de-crosslinked overnight at 65°C. The next day, samples were treated with proteinase K [1 μl of 20 mg/ml Proteinase K, plus 4 μl 0.5 M EDTA, and 8 μl Tris-HCl pH 6.5] for 1 h at 45°C. ChIPed DNA and inputs were purified using the QIAquick PCR purification kit (Qiagen) and eluted in 60 μl. The following antibodies were used in the ChIP experiments: H3K27me3 (Millipore, #07–449); H3K4me3 (Diagenode, C15410003); H3K4me1 (Abcam, ab8895); H3K27Ac (Millipore, #07–360); H3 (Abcam, Ab1791); H3K36me3 (Abcam, ab9050); H3K27me1 (Active Motif, #61015); H3K27me2 (Cell Signaling, #9728); H3K79me2 (Abcam, ab3594); H2Bub (Cell Signaling, #5546); and H4K20me3 (Abcam, ab9053). Library preparation for ChIP-seq experiments was performed at the UPF/CRG Genomics Unit. Libraries were sequenced using Illumina HiSeq2000 sequencer.

Input datasets

Raw files and processed data from experiments performed in this study are available at the Gene Expression Omnibus (GEO) under the accession number GSE150633. Raw data of multiple samples from the literature was downloaded and reanalyzed to be included in the study. RNA-seq data of mouse ESCs was extracted from a previous publication from our lab (GEO accession number: GSE79606) [23]. ChIP-seq data of p300 in ESCs was obtained via GEO (GEO accession number: GSE89211) [21]. ChIP-seq data of H3K27me3, H3K4me3, H3K27ac, H3K4me1, and H3K36me3, and RNA-seq data of cardiac differentiation (mesoderm, cardioprecursors and cardiomyocytes), were obtained from https://b2b.hci.utah.edu/gnomex/ (accession numbers: 44R and 7R2) [35]. ChIP-seq data of H3K27me3, H3K4me3, H3K27ac, H3K4me1, and H3K36me3, RNA-seq data of neural differentiation (neural precursors and cortical neurons), and Hi-C data of ESCs were retrieved from GEO (GEO accession number: GSE96107) [30]. PRO-seq data of mouse ESC was obtained from a previous publication from our lab (GEO accession number: GSE99530) [51]. ChIP-seq data of H3K27me3, H3K4me3, H3K27ac, H3K4me1, and H3K36me3, and RNA-seq of mouse developmental stages (heart tissue from 10.5 embryonic day, liver tissue from 11.5 embryonic day, neural tube tissue from 12.5 embryonic day, kidney tissue from 14.5 embryonic day, lung tissue from 15.5 embryonic day) were obtained from ENCODE project [41]. The list of ENCODE accession numbers can be found in S10 Table. When replicates were available, pooling was done except for the ChIP-seq samples of H3K4me3 of neural precursors (replicate 1 was used) and H3K27ac of neural precursors (replicate 2 was used).

ChIP-seq analysis

The sequence reads of ChIP-seq data were mapped to the mm10 version of the mouse genome with the BOWTIE software [52], setting the option–m 1, which eliminates reads that align in more than one region. The ChIP-seq profiles were obtained using the function buildChIPprofile from SeqCode (https://github.com/eblancoga/seqcode). For the p300 ChIP-seq, peak calling against input was performed using MACS [53] with the option—shiftsize 100, which shifts tags to their midpoint. Information about the total number of reads and read mapping of each ChIP-seq experiment produced in this study can be found in Table 4.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Information on ChIP-seq experiments produced in this study.

https://doi.org/10.1371/journal.pcbi.1009368.t004

Chromatin segmentation

ChromHMM [54] was used to obtain a chromatin segmentation model for ESCs using the default parameters. The input data were ChIP-seq experiments of H3K4me3, H3K27me3, H3K27ac, and H3K4me1, using ChIP-seq of H3 as control. First, the function BinarizeBam was used to binarize the input mapped data. Next, the LearnModel function was ran to learn different chromatin segmentation models of ESCs, using from 4 to 16 states; the 9-state model was selected because it showed the higher number of states with no redundancy.

RNA-seq analysis

The pair-end sequence reads of RNA-seq data were mapped to the mm10 version of the mouse genome with TopHat [55], setting the options—mate-inner-dist 100, which is the expected mean distance between mate pairs, and -g 1, which eliminates those reads which align in more than one region. The RNA-seq profiles were obtained using the function buildChIPprofile from SeqCode. The FPKMs (fragments per kilobase of transcript per million mapped reads) of each gene in the RefSeq catalogue [24] of the mouse genome were calculated using Cufflinks [56], setting the option—max-bundle-frags 5,000,000, which specifies the maximum genomic length for the bundles.

PRO-seq analysis

The single-end sequence reads of PRO-seq data was mapped to the mm10 version of the mouse genome with TopHat [55], setting the options—mate-inner-dist 100, which is the expected mean distance between mate pairs, and -g 1, which eliminates those reads which align in more than one region. The normalized count of reads of PRO-seq at intergenic enhancers averaged by the length of the region was calculated by recoverChIPlevels from SeqCode.

Hi-C analysis

Hi-C data were processed with TADbit [57]. Briefly, sequencing reads were mapped to the reference genome (mm10) by applying a fragment-based strategy, which is dependent on the GEM mapper [58]. Mapped reads were filtered to remove those resulting from unspecified ligations, errors, or experimental artefacts. Specifically, seven different filters were applied using the default parameters in TADbit: self-circles, dangling ends, errors, extra dangling-ends, over-represented, duplicated, and random breaks [57]. After pooling replicates, Hi-C data were normalized with OneD correction [59] at 5 kb of resolution to remove known biases. Significant Hi-C interactions were called with the analyzeHiC function of HOMER software suit [60], binned at 5 kb of resolution, and with the default p-value threshold of 0.001.

Gene expression predictive model

The regression linear models were built to predict gene expression by adjusting the following formula: where y_i is the log₂ of the FPKMs of gene i, with a pseudo count of 0.1. x_i1 to x_in are the log₂-normalized count of reads of each ChIP-seq signal at the defined promoters or enhancers averaged by the length of the region calculated by recoverChIPlevels from SeqCode, plus a pseudo count of 0.1. β₀ to β_n are the coefficients that we would like to calculate and ε is the error. The predictive models were trained on protein-coding genes. The set of data was randomly divided into two subsets, a training subset with the 80% of entries, and a test subset with the remaining 20% of entries. In the case of differentiation, each of the time points was used as training subsets and then the predictive models were evaluated in the rest. A 10-fold cross-validation was repeat three times to verify that the quantitative relationship between expression and histone modifications was not specific for a subset of the data. The following functions were used: trainControl to perform the 10-fold cross-validation, train to train the models, and varImp to calculate the variable importance, from the R package caret [45]. For models trained on enhancers, genes were introduced into the dataset as many times as the number of associated enhancers they had. To evaluate the specificity of our predictive models, we randomly shuffled the expression values of all the genes in the mouse genome. Thus, for each initial predictive model a random model was also obtained, where all the values y_i were shuffled, maintaining the values of x_i1 to x_in intact. This operation generates a new table of gene expression assignments in which the putative relationship between histone marking and expression of genes (if any) is completely lost. The random models were next generated following the same procedure as the models.

LOESS normalization

The FPKMs of all protein-coding genes and ChIP-seq levels of PEs and BPs were normalized for ESCs, cardiomyocytes, cortical neurons, cardio precursors, mesoderm, and neural precursors; and also, for heart tissue from 10.5 embryonic day, liver tissue from 11.5 embryonic day, neural tube tissue from 12.5 embryonic day, kidney tissue from 14.5 embryonic day, and lung tissue from 15.5 embryonic day. To normalize the ChIP-seq levels of H3K27me3, H3K4me3, H3K27ac, H3K4me1, and H3K36me3 on the PEs and BPs, the genome was first divided into 2 Kb bins (note that bin size reflects the average size of PEs and BPs). Next, the count of reads, normalized by total number of reads and averaged by the length, was calculated with recoverChIPlevels function from SeqCode. Finally, the normalization parameters were calculated in those bins and applied to the count of reads normalized by the total number of reads and averaged by the length of PEs and BPs. The normalize.loess function of the R package affy [61] was used to normalize ChIP-seq data and expression data. Genes and bins with a 0 in any columns were discarded, as it was not possible to determine whether it was due to a sequencing error or a real absence of signal.