14 Case study: H3K27ac ChIP-seq in the heart
Chapters 4 to 8 used one ChIP-seq dataset to show each step: the first run, read processing, quality control, peaks and signal tracks. This chapter takes the same dataset further. It asks the biological question, and follows EpiNexus’s answer through differential analysis, Gene Ontology, motifs, super-enhancers and enhancer–gene links.
14.1 The question
Arid1a encodes a subunit of the SWI/SNF chromatin-remodelling complex (Boogerd et al. 2023). This complex moves nucleosomes, so that other proteins can reach the DNA. Boogerd and colleagues removed Arid1a from heart muscle cells (cardiomyocytes) in mice. They used a Cre enzyme that is only made in cardiomyocytes (αMHC-Cre), and studied the hearts at postnatal day 7 (Boogerd et al. 2023).
H3K27ac marks active enhancers and promoters (Chapter 1). So H3K27ac ChIP-seq can show which regulatory regions become more or less active without Arid1a. We ask:
- Which regions lose or gain H3K27ac in the knockout?
- Which transcription factors bind there?
- Which genes might these regions control?
14.2 The dataset and the run
GEO series GSE155100 has four H3K27ac libraries (control and knockout, two replicates each) and one input per library (Boogerd et al. 2023). The reads are single-end, 75 bp. Chapter 2 describes the design and Chapter 4 the run.
The run used the default ChIP-seq settings: MACS3 with q 0.05 against each sample’s own input, TMM normalisation, and differential regions at FDR < 0.1 and |log2 fold change| > 0.5. The comparison was aMHC-Cre;Arid1a f/f cKO vs Arid1a f/f Control, so a positive fold change means “higher in the knockout”.
The run in this chapter is cafc042f, made with EpiNexus 1.2.0 by Run again on the first run of this project (93a23cb5, EpiNexus 1.1.0, used in Chapters 5 to 8). Quality control, peaks, super-enhancers and enhancer–gene links are the same in both runs. The differential results differ, because version 1.2 changed how H3K27ac regions are defined and tested.
14.3 Quality control and peaks
The quality checks are described in detail in Chapter 6 (Table 6.2). In short:
- All four H3K27ac libraries have a FRiP of 36–45 %, so most reads fall in peaks.
- EpiNexus marks each library with a warning, not a failure: 15–17.5 million reads are below ENCODE’s target of 20 million.
- The inputs show no enrichment, as they should.
MACS3 found 54,372 to 66,954 peaks per library (Chapter 7). Where are they? EpiNexus assigned every peak to a genomic feature (Table 14.1).
| Feature | Peaks | Share |
|---|---|---|
| Intron | 104,211 | 43 % |
| Promoter (≤ 1 kb) | 52,503 | 22 % |
| Distal intergenic | 44,737 | 18 % |
| Promoter (1–3 kb) | 26,319 | 11 % |
| Exon, UTRs, downstream | 14,745 | 6 % |
About a fifth of the peaks lie within 1 kb of a gene start, and a third within 3 kb, the promoter window EpiNexus uses. Most of the rest lie in introns and between genes, where enhancers are.
EpiNexus also tested each library’s peaks for Gene Ontology terms (Chapter 11). The top terms were “focal adhesion”, “actin filament organization”, “contractile muscle fiber” and “myofibril”. These describe heart muscle. They tell us that the peaks are where we expect them in heart tissue. They say nothing about Arid1a: a list of all peaks reflects the tissue, not the experiment.
GO terms of a sample’s full peak list describe the cell type. To learn about the experiment, test the regions that changed between conditions.
14.4 Differential H3K27ac
H3K27ac often covers a few kilobases around an active enhancer. So, for H3K27ac and other sharp histone marks, EpiNexus 1.2 tests whole peak regions, not short windows. It calls peaks once on all four H3K27ac libraries pooled together, against the four pooled inputs. Calling peaks on pooled libraries keeps the error rate of the later test under control (Lun and Smyth 2014). EpiNexus then merges the peaks and splits regions longer than 5 kb. This gave 92,052 regions, with a median width of 710 bp. It counted the reads in each region, and tested the difference with DESeq2 (Love et al. 2014) (Chapter 9).
2,461 regions changed: 2,267 lost H3K27ac in the knockout and 194 gained it (Figure 14.1). Losses outnumber gains more than ten to one. Without Arid1a, many regulatory regions of the heart become less active.
The ten strongest changes are all losses. Two lie near genes for ion channels: the potassium channel gene Kcnj12 (log2 fold change −3.9) and the calcium channel gene Cacna1g (−1.7). Others lie near genes with little known function, such as 4631405J19Rik. The strongest region of all, near 4631405J19Rik, has an FDR of 6 × 10−43.
With two replicates in each group, EpiNexus warns that outlier samples cannot be detected, and that the p-values lean on the trend across all regions (Chapter 2). The strongest changes are the most likely to be real. Changes near the threshold are less certain. More replicates would find more of them.
14.5 Gene Ontology of the changed regions
EpiNexus tested the lost and the gained regions for GO terms, against all tested regions (Chapter 11).
- Lost regions (2,267): 284 significant terms. The top one is “response to denervation involved in regulation of muscle adaptation”.
- Gained regions (194): 5 significant terms, such as “type 2 fibroblast growth factor receptor binding”.
Open the GO table before you read these terms. For each term it shows how many regions and how many genes stand behind it, and the result of a second test that counts genes instead of regions (Chapter 11).
- The top term for the lost regions comes from 26 regions near only four genes: Dag1, Fbxo32, Hdac4 and Scn5a. Its gene test is not significant (FDR 0.12).
- The top term for the gained regions comes from four regions near one gene, Fgf7.
A term that rests on one or a few genes says more about those genes than about a pathway. Stronger terms pass both tests. For the lost regions, 38 terms do, and the top ones describe the heartbeat (Table 14.2).
| GO term (lost regions) | Regions | Genes | FDR, gene test |
|---|---|---|---|
| cell communication involved in cardiac conduction | 37 | 15 | 0.0002 |
| cardiac muscle cell contraction | 47 | 16 | 0.0006 |
| cardiac muscle cell action potential involved in contraction | 41 | 13 | 0.002 |
| ventricular cardiac muscle cell action potential | 28 | 10 | 0.002 |
So, without Arid1a, H3K27ac falls near genes that conduct the electrical signal of the heart and make it contract. The gained regions are few. Three of their five terms rest on a single gene, and the other two on three or four genes. They do not form a clear theme.
14.6 Motifs in the changed regions
EpiNexus looks for known motifs in the lost and the gained regions. It compares them with regions of the same GC content that did not change, using SEA (Bailey and Grant 2021) and the motifs of JASPAR 2024 (Rauluseviciute et al. 2024) (Chapter 12). For histone regions it tests 500 bp around the middle of each region, because the factors bind between the nucleosomes that carry the mark.
| Regions | Significant motifs | Top motif families (best motif) |
|---|---|---|
| Lost (2,267) | 48 | MEF2 (MEF2D, q 10−25); steroid hormone receptors (ESRRB, q 2 × 10−8); Tal-related (FIGLA, q 6 × 10−8) |
| Gained (194) | 1 | Hairy-related (HES2, q 0.03) |
| Lost, against gained | 2 | Nuclear factor 1 (NFIB, q 0.002) |
The lost regions carry MEF2 motifs (Table 14.3). MEF2 is a family of transcription factors best known for controlling muscle genes. The ESRR motifs belong to nuclear receptors that control energy production in cells. So the regions that lose activity without Arid1a may be bound by muscle and metabolic regulators.
A second method, monaLisa (Machlab et al. 2022), sorts all tested regions into 17 groups by their fold change, and tests each group against the others. It agrees. The MEF2 motifs are enriched among regions with moderate losses (log2 fold change between −0.8 and −1.1). The ESRR motifs are enriched among the regions that lost the most.
The gained regions show almost nothing. With 194 regions, one motif barely passes the threshold. A small set has little power, and a weak result here is not evidence that no factor binds.
A motif shows which factors could bind. It does not show that they do, or that they cause the change. To show binding, you need ChIP-seq or CUT&RUN for the factor itself.
14.7 Super-enhancers
Some enhancers are much stronger than all others, and they often lie near genes that define the cell type. Whyte and colleagues called them super-enhancers (Whyte et al. 2013). EpiNexus finds them in each library in three steps (Chapter 13):
- It joins peaks that lie within 12.5 kb of each other, leaving out peaks within 2.5 kb of a gene start.
- It ranks the joined regions by their H3K27ac signal.
- It draws the ranked signal as a curve (Figure 14.2). The curve is flat for most regions and then rises steeply. Regions above the bend are super-enhancers.
| Library | Joined regions | Super-enhancers |
|---|---|---|
| Control, rep 1 | 14,292 | 767 |
| Control, rep 2 | 16,254 | 793 |
| Knockout, rep 1 | 13,171 | 715 |
| Knockout, rep 2 | 14,816 | 977 |
The strongest super-enhancers in the first control library lie near Sorbs1, Tns1, Rbm20, Gata4 and Mical2. Gata4 encodes a key transcription factor of the heart, and Rbm20 controls the splicing of heart muscle genes. So the strongest enhancers point to the identity of the tissue, as expected.
Two cautions:
- Counts do not compare conditions. The two knockout libraries have 715 and 977 super-enhancers (Table 14.4). The difference between replicates is larger than the difference between conditions. To compare conditions, use the differential analysis.
- The nearest gene is a guess. Several top regions are named after genes such as 4930449E18Rik or Gm26978. These are genes with little or no known function that happen to lie closest. The gene that an enhancer controls can lie further away.
14.8 Enhancer–gene links
To guess which genes an enhancer controls, EpiNexus uses the Activity-by-Contact (ABC) model (Fulco et al. 2019) (Chapter 13). It scores each pair of enhancer and gene by the enhancer’s activity (here its H3K27ac reads) times its expected contact with the gene. The contact is estimated from the distance. Each library got 111,126 to 129,705 links to 32,416–35,889 genes.
These links are an approximation. The full ABC model also uses chromatin accessibility (ATAC-seq or DNase-seq). In a large test against CRISPR experiments in human cells, ABC with H3K27ac alone found true enhancer–gene pairs only a little better than distance alone, and less well than ABC with DNase-seq data (Gschwind et al. 2026). Use the links to choose candidates, not as proof.
Here is one example. Chapter 8 showed the gene Atp1a2 in the genome browser (Figure 8.2). Its strongest link in both control libraries is a region about 2 kb from its start site, with an ABC score of 0.17. In the knockout, this region has much less H3K27ac: its activity falls from 7.9 and 10.6 in the controls to 1.2 in the first knockout library, and in the second knockout library it is no longer linked at all.
14.9 What this case study shows
- Check the basics first. Good QC and peaks in the expected places (Chapter 6, Chapter 7) give confidence in everything that follows.
- Without Arid1a, regulatory regions mostly lose activity. 2,267 regions lost H3K27ac and 194 gained it.
- The lost regions point to the heartbeat. They carry MEF2 and ESRR motifs, and lie near genes for cardiac conduction and contraction.
- Read GO terms through their genes. Terms that rest on one or a few genes are weak evidence; terms that also pass the gene test are stronger.
- Ask questions of the changes, not of all peaks. GO terms of all peaks describe the tissue.
- Super-enhancers point to cell identity. But compare conditions with the differential analysis, not with super-enhancer counts.
- Enhancer–gene links are candidates. With H3K27ac alone, they are an approximation of the ABC model.