17  Case study: DNA methylation

This chapter follows one RRBS experiment through EpiNexus: the library type, quality control, differentially methylated regions (DMRs) and where they lie. It also explains how dmrseq decides what is significant, and why a careful analysis can find no DMRs at all.

17.1 The question

DNA methylation is added to CpGs by enzymes called DNA methyltransferases. Two of them, Dnmt3a and Dnmt3b, add methyl groups to CpGs that were not methylated before. This is called de novo methylation (Zocher et al. 2021). When the stem cells of the adult brain turn into neurons, which CpGs need these two enzymes? We compare neurons made from stem cells that lack both enzymes with neurons made from normal stem cells.

17.2 The dataset

The data come from a study by Zocher and colleagues (GEO GSE167954) (Zocher et al. 2021). In their mice, the Dnmt3a and Dnmt3b genes can be cut out by an enzyme called Cre. Some mice also carry the Cre gene (Nestin::CreERT2), which becomes active in neural stem cells after a drug, tamoxifen. The authors grew neural stem cells from the hippocampus of both kinds of mice, turned them into neurons in culture, and sorted the neurons. They then made RRBS libraries (Chapter 1) from three cell lines per genotype (Table 17.1).

Table 17.1: The six RRBS libraries used in this chapter. Single-end reads of 76 bases, from the European Nucleotide Archive.

Sample	GEO	Genotype	Reads
Neuron_WT1	GSM5116896	control	13.86 million
Neuron_WT2	GSM5116897	control	10.41 million
Neuron_WT3	GSM5116898	control	13.93 million
Neuron_KO1	GSM5116899	without Dnmt3a/b	13.54 million
Neuron_KO2	GSM5116900	without Dnmt3a/b	15.66 million
Neuron_KO3	GSM5116901	without Dnmt3a/b	15.43 million

Together the files are about 3.0 GB. The study used mm10, so we use it too.

NoteIn EpiNexus 1.2

1. Create a project with the genome Mouse (mm10), and import the six GSM numbers. The import uses GEO’s genotype labels, such as nestin::CreERT2(+/-)-Dnmt3a-2lox-Dnmt3b-2lox for the knockout.
2. Start a run: DNA Methylation, Full Pipeline, genome mm10. On the Samples & Genome step, make one group per genotype and check the comparison: knockout first, control second.
3. On the Parameters step, keep Library type at Detect from the reads (recommended). Under Advanced Parameters, the defaults are Minimum reads per CpG, in every sample = 4 and Region FDR (q-value) = 0.05. This run used 10 reads.

The run used in this chapter is c6a0187f, made with EpiNexus 1.2.4. It is a Run again of an earlier run made with EpiNexus 1.2.2; both found the same numbers of candidates, swap regions and DMRs (Chapter 18).

Figure 17.1 and Figure 17.2 show the two steps of the run wizard.

Figure 17.1: Groups & Contrasts on the Samples & Genome step. Each sample’s group is its genotype label from GEO; the labels are long, so the boxes cut them off. The reference is the control group (-/-), so the one contrast is the knockout (+/-) against the control. EpiNexus 1.2 interface; not covered by the book’s licence.

Figure 17.2: The Parameters step for a DNA-methylation run. Library type is set to detect RRBS or WGBS from the reads. Under Advanced Parameters are the two settings of the DMR test, here at their defaults: 4 reads per CpG in every sample and a region FDR of 0.05. The run in this chapter used 10 reads. EpiNexus 1.2 interface; not covered by the book’s licence.

17.3 RRBS or WGBS?

The two kinds of bisulfite library need different processing (Chapter 5), so EpiNexus first checks which kind it has. RRBS cuts DNA with the enzyme MspI at CCGG sites, so “the first three bases of (almost) all RRBS reads are either CGG or TGG” (Krueger 2026b). EpiNexus reads the start of each sample’s first reads. Here 90–96 % of reads start with CGG or TGG, so EpiNexus treated the samples as RRBS. (EpiNexus calls a sample RRBS above 50 % and WGBS below 20 %; these lines are its own.)

RRBS changes two steps:

• Trimming. Trim Galore runs with --rrbs, which removes the bases that library preparation fills in at the end of each fragment (Krueger 2026b).
• No duplicate removal. Every RRBS fragment starts at an MspI site, so reads from different DNA molecules start at the same places. Removing “duplicates” would throw real data away. Bismark’s authors write that “deduplication is not recommended for RRBS, amplicon or other target enrichment-type libraries” (Krueger 2026a).

17.4 Quality control

Table 17.2: Bisulfite QC from EpiNexus’s report (run c6a0187f). MspI starts: reads starting with CGG or TGG. Mean depth: reads per covered CpG. mCpG, mCHH: share of methylated calls in each context.

Sample	MspI starts	Reads	Mapped	Mean depth	mCpG	mCHH	Conversion
control 1	92 %	13.8 M	73.0 %	15.9×	22.1 %	0.7 %	99.3 %
control 2	96 %	10.4 M	71.1 %	12.3×	27.5 %	0.7 %	99.3 %
control 3	91 %	13.9 M	73.5 %	17.5×	18.8 %	0.7 %	99.3 %
knockout 1	90 %	13.5 M	71.2 %	12.3×	25.2 %	0.7 %	99.3 %
knockout 2	92 %	15.6 M	70.9 %	14.1×	26.2 %	0.7 %	99.3 %
knockout 3	91 %	15.4 M	70.4 %	13.5×	26.7 %	0.7 %	99.3 %

The libraries are good (Table 17.2):

• Conversion is 99.3 % in every sample. Bisulfite must convert every unmethylated cytosine. In mammals, cytosines outside CpGs (CHH) are almost never methylated, so the share of “methylated” CHH (0.7 %) shows the cytosines that bisulfite missed. ENCODE asks for a conversion rate of at least 98 % (ENCODE Project Consortium 2026).
• Depth. Each covered CpG was read 12 to 18 times on average.
• Methylation values have the expected shape: in every sample, most CpGs are near 0 % methylated and a second group is near 100 %.
• Replicates. In the sample PCA, the first component (76.2 % of the variation) separates the two genotypes.

The Quality Control page shows the same table (Figure 17.3).

Figure 17.3: The bisulfite QC table on the Quality Control page (run c6a0187f). The table is wider than the page and is scrolled, so the sample names (Neuron_KO1 … Neuron_WT3) are cut at the left. The line above it says that the library is RRBS and was not deduplicated; Conversion is graded against ENCODE’s 98 %. EpiNexus 1.2 interface; not covered by the book’s licence.

WarningPitfall: comparing samples by their overall mCpG

The mCpG column ranges from 18.8 % to 27.5 %. The three lowest values and the three highest values each mix both genotypes. This number averages all CpG calls, weighted by reads. RRBS reads gather in CpG islands, which are mostly unmethylated. So the value moves with where each library’s reads fell, not only with methylation. Compare samples CpG by CpG, as the DMR test does.

The M-bias plot shows methylation by position along the read (Figure 17.4). After the first few bases, it is flat, as it should be. The first three positions are higher in every sample. These positions read the CpG of the MspI site where the DNA was cut: “every read provides information on at least one CpG right in the start” (Krueger 2026b). They measure a different set of CpGs from the rest of the read, and the pattern is the same in all six samples.

Figure 17.4: M-bias: CpG methylation by position in the read, for each sample. After the first three positions the lines are flat at about 15–25 %. Figure made by EpiNexus 1.2.4 (run c6a0187f).

17.5 Calling DMRs with dmrseq

EpiNexus finds differentially methylated regions with dmrseq (Korthauer et al. 2019). It works in four steps:

1. Choose the CpGs to test. EpiNexus sums the two strands of each CpG, then keeps CpGs with at least 10 reads in every sample. That left 520,396 CpGs. With 1 read, 1.2 million CpGs would be covered in every sample; each extra read required removes some.
2. Find candidate regions. dmrseq smooths the methylation difference along the genome and looks for runs of at least 5 CpGs where the difference is above 10 points. It found 4,132 candidates.
3. Build a null. dmrseq swaps samples between the two groups, so that each new “group” mixes both genotypes, and looks for candidates again. Any region it finds then cannot come from the genotype. With three samples per group, there are 9 such swaps. They gave 1,656 to 2,104 candidates each, 17,298 in all.
4. Compare. Each real candidate gets a p-value: the share of the 17,298 swap regions that are at least as strong. The p-values are corrected for testing many regions, and regions with q ≤ 0.05 are DMRs: 3,188.

Step 4 has a floor. With 17,298 swap regions, the smallest possible p-value is 1 / 17,299, about 0.00006. The strongest candidates all reach it. That is why the top of the volcano plot is flat (Figure 17.5): 289 DMRs share the smallest q-value, 0.0008.

Figure 17.5: EpiNexus’s DMR result: each point is a candidate region. Blue: lower methylation without Dnmt3a/b (3,148); orange: higher (40); grey: not significant (944). The flat rows at the top are regions that reached the smallest possible p-value. Figure made by EpiNexus 1.2.4 (run c6a0187f).

TipBest practice

Test regions with a method that accounts for the differences between replicates and controls the false discovery rate for regions. dmrseq builds its null from sample swaps, “a pooled null distribution, that can be implemented even when as few as two samples per population are available”, and “accurately controls the false discovery rate” (Korthauer et al. 2019). No publication sets the minimum reads per CpG for such a test. EpiNexus’s default is 4 reads in every sample.

17.6 Results: focal loss of methylation

3,188 DMRs: 3,148 lose methylation without Dnmt3a/b, and 40 gain it. The losses are large: the median DMR loses 36 points, and 2,814 DMRs change by 20 points or more. The median DMR is small: 7 CpGs over 163 bp.

Yet the average over all covered CpGs barely moves: 46.3 % in control neurons, 44.3 % without the enzymes (Figure 17.6). The loss is focal. Dnmt3a/b methylate specific places, not the whole genome.

Figure 17.6: The Overview tab of the Methylation page (run c6a0187f): 3,188 DMRs (40 gained, 3,148 lost), 520,396 CpGs tested with at least 10 reads in every sample, the mean methylation of each group, and the CpG-island context of the DMR CpGs against all tested CpGs, with odds ratios. EpiNexus 1.2 interface; not covered by the book’s licence.

Where. EpiNexus places every DMR on the CpG-island map of the genome: an island, a shore (up to 2 kb from an island), a shelf (2–4 kb) or the open sea beyond (Irizarry et al. 2009; Cavalcante and Sartor 2017). It compares the CpGs in DMRs with all tested CpGs (Figure 17.7). In RRBS, most tested CpGs lie in islands (64 %). The CpGs that lose methylation are less often in islands (49 %) and twice as often in shores (24 % against 12 %).

Figure 17.7: CpG-island context of the tested CpGs (grey), the CpGs in DMRs that lose methylation (blue) and those that gain it (orange). Figure made by EpiNexus 1.2.4 (run c6a0187f).

Genes. Some of the strongest losses lie in CpG islands of genes that neurons use, such as Egr3 (21 CpGs, −78 points) and Neurod1 (−63 points), and of Hic1 (24 CpGs, −81 points). These regions lose most of their methylation without Dnmt3a/b.

The gains. 22 of the 40 gains lie on the X chromosome, 21 of them in CpG islands of X-linked genes such as Zic3, Efnb1 and Rps4x. Of the 3,148 losses, only 28 (0.9 %) are on the X. All six samples come from female mice (GEO). In female cells, one X chromosome is silenced, and the CpG islands on the silent X are methylated. So gains on the X more likely reflect the X-chromosome state of the cell lines than the missing enzymes. This is an interpretation, not a test.

WarningPitfall: changes on the sex chromosomes

When changes cluster on the X or Y chromosome, check the sex and the X-chromosome state of each sample before reading them as an effect of your treatment. Cultured female cells can lose an X chromosome or partly re-activate the silent one. EpiNexus 1.2 does not yet check this; look at the DMRs per chromosome on the Methylation page.

NoteIn EpiNexus 1.2

The Methylation page shows the DMR table, the volcano plot, the CpG-island context and a Samples tab with the sample PCA and correlation. The bisulfite QC table and the M-bias plots are on the Quality Control page. The export has the DMRs (dmrs.csv), every candidate region with its p- and q-value (candidates.csv), the figures, and each sample’s CpG calls.

17.7 When there are no DMRs

A careful analysis can find no DMRs. The most common reason is the design: few samples, or replicates that differ from each other.

dmrseq can only use swaps that mix the genotypes. With three samples per group there are 9, as here. With two samples per group there are only 2. The null then holds the regions of only two swaps, and a real region must stand out against both. If the two replicates of a group differ, the swaps produce strong regions of their own, and few real candidates beat them.

EpiNexus warns about this on the Review step of the run wizard: “… has 2 samples per group: dmrseq can compare it with only 2 sample-swapped versions of the data”. If a run finds no DMRs:

• Look at the candidates. EpiNexus lists the strongest candidate regions, marked as not significant, with the counts of candidates and swap regions. A real change that is too weak to pass shows up there.
• Look at the replicates. If the PCA does not separate the groups, the swaps find as much as the real comparison does.
• Add samples. A third sample per group raises the number of swaps from 2 to 9.

A test that pools the reads of each group, and ignores how much replicates differ, reports more differences. But it cannot tell them from differences between replicates of the same group (Chapter 2).

17.8 What this case study shows

• The library type matters. RRBS needs --rrbs trimming and no duplicate removal; EpiNexus detects it from the reads.
• Check conversion first. 99.3 % here, above ENCODE’s 98 %.
• Losses of de novo methylation are focal. 3,148 regions lose a median of 36 percentage points of methylation, while the average over all CpGs moves by 2 points.
• Read a DMR test through its null. The swaps set the smallest p-value, and with two samples per group there are only two swaps.
• Check the sex chromosomes before interpreting changes there.

17.9 Further reading

• Zocher et al. (2021): the study that produced these data, with the role of de novo methylation in adult neurogenesis.
• Korthauer et al. (2019): dmrseq and its permutation test.
• Meissner et al. (2008): RRBS.
• Irizarry et al. (2009): CpG-island shores.

Cavalcante, Raymond G., and Maureen A. Sartor. 2017. ‘Annotatr: Genomic Regions in Context’. Bioinformatics 33 (15): 2381–83. https://doi.org/10.1093/bioinformatics/btx183.

ENCODE Project Consortium. 2026. Whole-Genome Bisulfite Sequencing (WGBS) Data Standards. Https://www.encodeproject.org/data-standards/wgbs/.

Irizarry, Rafael A., Christine Ladd-Acosta, Bo Wen, Zhijin Wu, Carolina Montano, et al. 2009. ‘The Human Colon Cancer Methylome Shows Similar Hypo- and Hypermethylation at Conserved Tissue-Specific CpG Island Shores’. Nature Genetics 41 (2): 178–86. https://doi.org/10.1038/ng.298.

Korthauer, Keegan, Sutirtha Chakraborty, Yuval Benjamini, and Rafael A. Irizarry. 2019. ‘Detection and Accurate False Discovery Rate Control of Differentially Methylated Regions from Whole Genome Bisulfite Sequencing’. Biostatistics 20 (3): 367–83. https://doi.org/10.1093/biostatistics/kxy007.

Krueger, Felix. 2026a. Bismark Documentation: Deduplication. Https://felixkrueger.github.io/Bismark/usage/deduplication/.

Krueger, Felix. 2026b. Trim Galore: RRBS Guide, Directional Libraries. Https://trimgalore.com/rrbs/directional/.

Meissner, Alexander, Tarjei S. Mikkelsen, Hongcang Gu, Marius Wernig, Jacob Hanna, et al. 2008. ‘Genome-Scale DNA Methylation Maps of Pluripotent and Differentiated Cells’. Nature 454 (7205): 766–70. https://doi.org/10.1038/nature07107.

Zocher, Sara, Rupert W. Overall, Gabriel Berdugo-Vega, Nicole Rund, Anne Karasinsky, et al. 2021. ‘De Novo DNA Methylation Controls Neuronal Maturation During Adult Hippocampal Neurogenesis’. The EMBO Journal 40 (18): e107100. https://doi.org/10.15252/embj.2020107100.