Functional annotation & region-level analysis

Annotation and Interpretation

At this point the EWAS is done: we have a per-CpG table of effect sizes and p-values, BACON-adjusted for residual inflation (notebook 06), and — if we chose the stratified route — a METAL meta-analysis of the strata (notebook 07). Next we introduce functional annotation and region-level analysis. Previously it was mentioned that the snakemake pipeline can do differentially methylated region analysis in addition to the EWAS. But, it wasn’t really explained what that is or why you would want to do it. Here we more fully detail functional annotation and secondary analysis (including regional analyses) that are a necessary part in interpretation of results.

1. CpG-level functional annotation — attach genomic coordinates, gene context, CpG-island relationship, and expression links to every probe, so the top hits can be read as “this CpG sits in the promoter of gene X” rather than “cg16340178”.
2. Region-level analysis (DMRs) — ask whether spatially clustered CpGs move together. A differentially methylated region (DMR) can be highly significant even when no single CpG in it survives genome-wide correction, because neighboring probes are correlated and that correlation carries signal (Pedersen et al. 2012; Peters et al. 2015).

For functional annotation of the EWAS results we use Wanding Zhou’s hg38 EPIC annotation resources (Zhou et al. 2017). Zhou recently (July of 2026) made a major update of his website and released new versions of his resources (v8.1). This tutorial implements some of the newest version of Zhou’s files in addition to some ‘legacy’ files. The Snakemake pipeline uses entirely legacy files at the moment, but will be updated to use the latest versions in the next release. We only use hg38 annotations in this tutorial, but hg19 resources are available if necessary. The three tables used for annotation, all downloadable from the zhou-lab GitHub are:

• EPIC.hg38.manifest.gencode.v41.tsv.gz — gene model: gene names, transcript biotypes, and signed distance to the nearest TSS.
• EPIC.hg38.manifest.gencode.v36.tsv.gz — the v36 table additionally carries the CpG-island (CGI) relation column we use for island context.
• EPIC.hg38.coord.tsv.gz — remapped hg38 coordinates. Note this table is keyed by genomic position and carries no probeID column, so it cannot be merged onto results by probe; the v41 manifest above already supplies CpG_chrm/CpG_beg/CpG_end and is what we actually join.

Methylation is often a mediator of change. The environment can cause changes to methylation, which subsequently causes change in gene expression. Methylation sites that cause changes in gene expression are called expression quantitative trait methylation sites (eQTMs). So, one additional external layer, the BIOS cis-eQTM catalog (whole-blood methylation→expression links) was added to expand the functional annotation. This helps identify which (if any) effector gene may be behind methylation change. I recommend reading more into the BIOS eQTM catalog more before using and reporting the results as there are some distinct limitations (only 450k methylation array probes were tested, European-ancestry only cohort, original analysis uses hg19 mapping). This part of the functional annotation was added to the EWAS Snakemake upon request from a reviewer for the information.

NoteThis notebook and the Snakemake pipeline

The krferrier/EWAS pipeline automates everything here — its annotate.smk and dmr.smk rules run the same annotation join and the same comb-p region detection we walk through manually. We reproduce the steps by hand so the logic is clearly layed out.

library(data.table)      # fread()/fwrite() and the := joins used all through this chapter
library(ggplot2)         # the two context bar charts in Part 1

Part 1 — CpG-level functional annotation

What annotation resources exist, and which to trust

An EPIC probe ID like cg16340178 is the Illumina name for a CpG. To interpret it we need to know where it lands in the genome and what functional element it tags. Historically people reached for the hg19 Illumina/minfi manifest for this, but its coordinates are a genome build behind current references and its gene mapping is generous (it lists any transcript overlapping the probe). Zhou’s resources supersede it:

Zhou hg38 table	What it gives	Column(s) we use
manifest.gencode.v41	gencode v41 gene model	genesUniq, transcriptTypes, distToTSS
manifest.gencode.v36	earlier gene model that also carries CGI relation	CGIposition (Island / N_/S_Shore / N_/S_Shelf / open sea)
coord	remapped hg38 coordinates	CpG_chrm, CpG_beg, CpG_end
BIOS cis-eQTM (Bonder et al. 2017) (external)	genes whose expression tracks the CpG’s methylation in whole blood	likely downstream effector gene

The annotation join

We start from the BACON-adjusted top-table saved in chapter 06 and left-join the Zhou tables on probe ID. The v41 gene table can list several TSS distances per probe (one per transcript); we keep the nearest, which is what defines the promoter call below.

NoteWhat this chapter needs

This chapter starts from data/06_ewas_bacon_toptable.csv.gz, written by chapter 06. If you did not run that chapter, fetch the published checkpoint instead:

Terminal

./get_data.sh F_ewas_results

# The two Zhou gencode manifests joined here are tens of MB each and are not
# tracked in the repository: they come from the zhou-lab release linked above,
# and `./get_data.sh H_annotation` puts them under ../masks/. Everything in
# Part 1 rendered below reads a table this chunk, or the eQTM chunk after it,
# writes.

# 1. EWAS results: the BACON-adjusted toptable from chapter 06
tt <- fread("data/06_ewas_bacon_toptable.csv.gz")

# 2. Zhou gencode v41 gene model: gene names, biotypes, signed distance to TSS
v41 <- fread(cmd = "zcat ../masks/EPIC.hg38.manifest.gencode.v41.tsv.gz")
# a probe with several transcripts has several distToTSS values (';'-joined);
# keep the one nearest the TSS — this is what the promoter definition uses
nearest_tss <- function(x) {
  v <- suppressWarnings(as.numeric(strsplit(x, ";")[[1]]))
  if (all(is.na(v))) return(NA_real_)
  v[which.min(abs(v))]
}
# one function call per probe over the whole manifest: the slow line here
v41[, distTSS := vapply(distToTSS, nearest_tss, numeric(1))]

# 3. CpG-island relation — only the v36 gene table carries CGIposition
v36 <- fread(cmd = "zcat ../masks/EPIC.hg38.manifest.gencode.v36.tsv.gz",
             select = c("probeID","CGIposition"))
v36[, island := fifelse(is.na(CGIposition) | CGIposition == "", "OpenSea", CGIposition)]

# 4. Join everything on probe ID. The v41 manifest already carries the hg38
#    coordinates (CpG_chrm/CpG_beg/CpG_end), so it serves as both the coordinate
#    and the gene-model source. Note that the separate `EPIC.hg38.coord.tsv.gz`
#    table is keyed by position only — it has no probeID column — so it cannot be
#    joined to results by probe and is not used here.
#    all.x = TRUE: every tested CpG stays in the table, annotated or not.
annotated <- tt |>
  merge(v41[, .(probe = probeID, CpG_chrm, CpG_beg, CpG_end,
                genesUniq, transcriptTypes, distTSS)], by = "probe", all.x = TRUE) |>
  merge(v36[, .(probe = probeID, island)], by = "probe", all.x = TRUE)

# 5. Zhou's gene-feature convention: the gene model runs from 1.5 kb upstream of
#    the TSS through the termination site. A probe within +/-1.5 kb of a TSS is
#    promoter-associated; genic but beyond that is gene body; no gene is intergenic.
annotated[, feature := fifelse(is.na(genesUniq) | genesUniq %in% c("", "NA"), "Intergenic",
                        fifelse(!is.na(distTSS) & abs(distTSS) <= 1500, "Promoter", "Gene body"))]
# rank by the BACON p-value once, here: every display table below takes its rows
# off the top of this object, and the merges above do not preserve order
setorder(annotated, bacon.p)

# The deposited table is written a few sections down, after the eQTM columns are
# joined on, so that the file on disk carries every annotation this chapter builds.

On our Grady PTSD subset this annotates all 756,251 CpGs with hg38 coordinates and maps 647,034 (85.6%) to at least one gencode v41 gene. The remaining ~14% are intergenic. This is expected since a large fraction of the EPIC array deliberately targets enhancers and open sea away from genes.

The annotated top hits

# The ten rows rendered below, derived from `annotated` rather than typed out so
# that the rounding, the column names and the ranking all come from the join.
t10 <- head(annotated, 10)[, .(
  probe, chr = CpG_chrm, pos = CpG_beg,
  # an intergenic probe has no symbol; print a dash rather than NA
  gene = fifelse(is.na(genesUniq) | genesUniq == "", "-", genesUniq),
  distTSS = round(distTSS), feature, island,
  # delta_beta is chapter 06's case−control β difference; ×100 puts it in
  # percentage points, which is how the prose discusses effect sizes
  delta_beta = round(100 * delta_beta, 2),
  bacon_p = signif(bacon.p, 3), bacon_FDR = signif(bacon.adj.P, 3))]
setnames(t10, c("probe","chr","pos","gene","distToTSS","feature","island",
                "\u0394\u03b2 (pp)","BACON p","BACON FDR"))
fwrite(t10, "data/08_annotation/top10_display.csv")

# written by the `top10-derive` chunk above
t10 <- fread("data/08_annotation/top10_display.csv")
knitr::kable(t10, caption = "Top 10 CpGs by BACON-adjusted p-value, with hg38 coordinates and gene/island context.")

Top 10 CpGs by BACON-adjusted p-value, with hg38 coordinates and gene/island context.

probe	chr	pos	gene	distToTSS	feature	island	Δβ (pp)	BACON p	BACON FDR
cg15434749	chr8	6391843	-	NA	Intergenic	OpenSea	-2.54	0.0e+00	0.0224
cg20600436	chr2	180911725	SCHLAP1	219622	Gene body	OpenSea	1.38	1.0e-07	0.0255
cg25717994	chr13	77919495	EDNRB;OBI1-AS1	-12	Promoter	S_Shore	0.57	3.0e-07	0.0627
cg02045051	chr8	134784106	ENSG00000279518;ENSG00000289405	1773	Gene body	OpenSea	-0.76	3.0e-07	0.0627
cg05825555	chr4	25926237	SMIM20	11929	Gene body	OpenSea	2.37	1.7e-06	0.2360
cg24228058	chr3	37620225	ITGA9	168085	Gene body	OpenSea	-1.26	2.0e-06	0.2360
cg18321729	chr15	72231987	PKM	-168	Promoter	S_Shore	-0.56	2.4e-06	0.2360
cg16340178	chr17	7315701	ENSG00000261915;GPS2;NEURL4	-137	Promoter	S_Shore	-0.89	2.5e-06	0.2360
cg26643724	chr5	81343724	ACOT12	479	Promoter	OpenSea	0.68	3.6e-06	0.2750
cg01824108	chr10	109457377	-	NA	Intergenic	S_Shore	-0.84	3.6e-06	0.2750

Now the list has biological context. The top hit cg15434749 is intergenic open sea on chr8 with no nearby gencode v41 gene. cg25717994 sits essentially on the EDNRB TSS (distToTSS = −12 bp, so promoter by Zhou’s ±1.5 kb rule); cg18321729 tags the PKM promoter (−168 bp); cg16340178 is in a CpG-island shore (S_Shore) just upstream of the GPS2/NEURL4 TSS (−137 bp). Four of the top ten fall within 1.5 kb of a TSS. The trend of significant CpGs being in a promotor region is mild in this analysis, but would likely be stronger in a full scale, well-powered EWAS. Only the top CpG clears the Bonferroni significance threshold (BACON p ≈ 3.0 × 10⁻⁸ against a threshold of ≈ 6.6 × 10⁻⁸); two pass FDR < 0.05. Read the rest of this table as ‘nominally significant’ rather than ‘statistically significant’.

In a real study this table is the hand-off point to the enrichment tools in Part 3 and to comparison against published EWAS catalogs.

Where the top hits sit relative to genes and islands

A useful sanity check is to compare the distribution of the top-ranked CpGs across gene features and island context against the array background. A well-powered signal often enriches in promoters and island shores.

# The two figures below, and the percentages their captions quote. Needs
# `annotated` from the join above; the PNGs it writes are the ones included
# underneath.
TOPN <- 1000
top  <- head(annotated, TOPN)          # `annotated` is already ordered by bacon.p

# One helper for both panels: tabulate a column in the top set and in the whole
# tested set, then convert to a percentage *within* each set, because the two
# sets differ in size by three orders of magnitude and raw counts would not be
# comparable.
mk <- function(col, lev) {
  a <- data.table(set = "Top 1,000",  grp = factor(top[[col]], levels = lev))
  b <- data.table(set = "All tested", grp = factor(annotated[[col]], levels = lev))
  d <- rbind(a, b)[!is.na(grp), .N, by = .(set, grp)]
  d[, pct := 100 * N / sum(N), by = set]
  # top set first, so the series the reader is asking about leads the legend
  d[, set := factor(set, levels = c("Top 1,000", "All tested"))]
  d
}
fd <- mk("feature", c("Promoter", "Gene body", "Intergenic"))
# island levels in genomic order: island, its two shores, its two shelves, then
# everything else — alphabetical order would scramble that
id <- mk("island", c("Island","N_Shore","S_Shore","N_Shelf","S_Shelf","OpenSea"))

ctx_plot <- function(d, xlab) {
  ggplot(d, aes(grp, pct, fill = set)) +
    geom_col(position = position_dodge(width = 0.8), width = 0.7) +
    scale_fill_manual(values = c(`Top 1,000`  = "#1A6B75",   # teal
                                 `All tested` = "#B8873F"),  # sand
                      name = NULL) +
    labs(x = xlab, y = "% of CpGs")
}
ggsave("data/08_annotation/08_feature_distribution.png",
       plot = ctx_plot(fd, "gene feature (Zhou gencode v41)"),
       width = 7, height = 3.2, dpi = 300)
ggsave("data/08_annotation/08_island_distribution.png",
       plot = ctx_plot(id, "CpG-island relation (Zhou CGI)"),
       width = 7, height = 3.2, dpi = 300)

Figure 1: Gene-feature context of the top 1,000 CpGs versus all 756,251 tested, using Zhou’s promoter (±1.5 kb of TSS) / gene-body / intergenic scheme. The top set leans slightly toward promoters (54.2% vs 52.2%) and away from gene bodies (31.7% vs 33.4%), with intergenic essentially unchanged (14.1% vs 14.4%) — the direction a real signal moves, but a small effect at this sample size.

Figure 2: CpG-island context, same comparison, using the Zhou CGI relation. A visible lean toward islands in the top set (22.3% vs 19.2%) and away from open sea (51.4% vs 55.6%). Hold this impression loosely: the formal test in Part 3 (Know Your CpG, another Zhou resource) does not find it significant.

The distributions move only modestly between the top 1,000 and the background (Figure 1, Figure 2). The top-set leans toward promoters and islands, but it’s not clear if that pattern is more than chance. The bar chart above is a useful visual indicator, but not a formal enrichment test, which is covered in Part 3 of this chapter.

The BIOS eQTM layer: methylation that tracks expression

Annotated methylation data often includes a column with the ‘nearest gene’ to the CpG. This is a useful annotation for identifying where in the genome the CpG is physically located and can hint at potential functional impact. For example, evidence of methylation at CpG in a promotor region of a gene could indicate a cis-effect between the CpG and expression of that gene. However, in reality, DNA sequence is not linear, it has 3-dimensional structure. In 3D space, the ‘nearest gene’ in terms of base-pairs may not be the gene that the methylation site is having an effect on. The answer we really want is is whether a CpG actually affects gene expression. The BIOS consortium answered this directly by pairing whole-blood methylation with RNA-seq in ~3,000 samples and testing every cis CpG–transcript pair, producing a catalog of cis-eQTMs (Bonder et al. 2017). This information provides the likely effector gene behind a methylation change — which is not necessarily the nearest gene.

NoteBIOS consortium meQTL catalogs

The BIOS consortium also provides catalogs for cis- and trans-meQTLs (methylation Quantitative Trait Loci), which provide evidence of genetic loci that cause changes to methylation status. This annotation information was not incorporated into the snakemake pipeline and how to merge it in will not be covered in this tutorial as it would require lifting over that data from hg19 to hg38 and only represents meQTLs for CpGs from the 450k array. If you have both methylation and genetic data for your sample, which is often the case, you can perform cis- and trans-meQTL testing yourself. The eQTM relationships require having gene expression and methylation data, which is less common and why the effort was made to liftover and update the meQTL data from the BIOS consortium.

The catalogs are keyed on Ensembl gene IDs; the pipeline’s prepare_bios_eqtm_annotation.R helper re-maps those to current HGNC symbols for eQTM data.

NoteWhat this chapter needs

This section reads data/08_annotation/bios_eqtm_hgnc_annotated.tsv, written by the pipeline helper shown in the chunk below. If you do not want to run the helper, fetch the published copy instead:

Terminal

./get_data.sh H_annotation

# Extends `annotated` from the join above. The two files it writes at the end
# are the ones every later section of this chapter reads.

# One-time prep: Ensembl -> current HGNC symbol (pipeline helper)
#   Rscript prepare_bios_eqtm_annotation.R \
#     --bios-eqtm 2015_09_02_cis_eQTMsFDR0.05-CpGLevel.txt \
#     --hgnc-complete-set hgnc_complete_set.txt \
#     --output bios_eqtm_hgnc_annotated.tsv
eqtm <- fread("data/08_annotation/bios_eqtm_hgnc_annotated.tsv")

# collapse to one row per CpG: effector gene(s) + strongest eQTM FDR
eq <- eqtm[HGNCName_GRCh38 != "" & !is.na(HGNCName_GRCh38),
           .(BIOS_eQTM_genes  = paste(sort(unique(HGNCName_GRCh38)), collapse = ";"),
             BIOS_eQTM_minFDR = min(as.numeric(FDR), na.rm = TRUE)),
           by = .(probe = SNPName)]

annotated <- merge(annotated, eq, by = "probe", all.x = TRUE)
setorder(annotated, bacon.p)      # merge() re-sorts by the join key; rank it again

# --- what Part 1 hands on -------------------------------------------------
# The CSV is the deposited annotation table (Zenodo tier H_annotation); the RDS
# is the same object, read back by the DMR per-CpG table in Part 2. Both are
# written here, after the eQTM merge, so neither can be missing a column the
# other has.
fwrite(annotated, "data/08_annotation/PTSD_ewas_annotated_zhou.csv.gz")
saveRDS(annotated, "data/08_annotation/annotated.rds")

On our subset the eQTM layer annotates 10,803 of the 756,251 tested CpGs (1.4%) with at least one expression-linked effector gene. I will reiterate, there are so few annotations because the BIOS eQTM catalog was curated using the 450k methylation array, and not all probes used in the 450k array were included in the 850k array. In addition, the BIOS catalog only keeps pairs at FDR < 0.05. That said, if we filter our results to look our top EWAS hits that do carry an eQTM, we see that the nearest-gene and effector-gene can differ:

# The twelve rows rendered below: the best-ranked hits that carry an eQTM.
eqd <- annotated[!is.na(BIOS_eQTM_genes)][1:12]
# Is the effector gene also the nearest gene? TRUE when any eQTM gene appears in
# the probe's own (semicolon-separated) gene list. mapply because both sides are
# strings that have to be split row by row.
eqd[, nearest := fifelse(mapply(function(g, e) {
      if (is.na(g) || g == "") return(FALSE)
      any(strsplit(e, ";")[[1]] %in% strsplit(g, ";")[[1]])
    }, genesUniq, BIOS_eQTM_genes), "yes", "no")]
eqdisp <- eqd[, .(probe,
                  nearest_gene = fifelse(is.na(genesUniq) | genesUniq == "", "-", genesUniq),
                  eQTM_gene = BIOS_eQTM_genes, `nearest?` = nearest,
                  eQTM_FDR = signif(BIOS_eQTM_minFDR, 3),
                  bacon_p = signif(bacon.p, 3))]
fwrite(eqdisp, "data/08_annotation/eqtm_display.csv")

# the same hits, untrimmed and with every annotation column, for anyone who
# wants to look further down the ranking than the table below shows
fwrite(annotated[!is.na(BIOS_eQTM_genes)][1:200],
       "data/08_annotation/eqtm_top_hits.csv")

# written by the `eqtm-top-derive` chunk above
knitr::kable(fread("data/08_annotation/eqtm_display.csv"),
  caption = "Top EWAS hits (by BACON p) that carry a BIOS cis-eQTM. The 'nearest?' column flags where the expression-linked effector gene is not the nearest gene by distance.")

Top EWAS hits (by BACON p) that carry a BIOS cis-eQTM. The ‘nearest?’ column flags where the expression-linked effector gene is not the nearest gene by distance.

probe	nearest_gene	eQTM_gene	nearest?	eQTM_FDR	bacon_p
cg26713179	ZNF766	ZNF528;ZNF528-AS1;ZNF880	no	0.000000	4.00e-06
cg12218747	CBR1-AS1;SETD4	CBR1	no	0.000000	1.80e-05
cg09138671	TMCO4	TMCO4	yes	0.000232	6.34e-05
cg25535873	HLA-DOA	HLA-DOA	yes	0.028900	1.23e-04
cg08632701	CBR1-AS1;SETD4	CBR1	no	0.000000	1.43e-04
cg27530239	CALCOCO2	CALCOCO2	yes	0.013700	3.19e-04
cg19351838	MMP17	MMP17	yes	0.000000	3.55e-04
cg26246807	ZIK1	ZNF154;ZNF772	no	0.001740	3.59e-04
cg11947245	-	SLC46A3	no	0.005990	4.59e-04
cg05281708	ZKSCAN7-AS1;ZNF35	ZNF35	yes	0.035000	6.13e-04
cg05226462	ZNF205;ZNF213-AS1	ZNF205	yes	0.000000	6.58e-04
cg01947224	ZIK1;ZNF416	ZIK1;ZNF134;ZNF154;ZNF211;ZNF304;ZNF416;ZNF419;ZNF530;ZNF547;ZNF549;ZNF551;ZNF671;ZNF772;ZNF773	yes	0.000000	7.54e-04

For cg06083200 the nearest gene is PPP1R11, but the CpG’s methylation tracks the expression of TRIM26. For cg26713179 the nearest gene is ZNF766 while the eQTM points at the neighboring ZNF528/ZNF880 cluster. Four of these twelve hits show an effector that differs from the nearest gene. When reporting results, the effector gene is usually the one people want to know about if the information is available.

Part 2 — Region-level analysis (DMRs)

Why look at regions at all

Single-CpG EWAS tests each probe independently, but methylation is not evenly spaced throughout the genome, and methylation arrays tend to target CpG islands and gene promotor regions where there are likely to be multiple methylation sites close together. Methylation sites in close proximity to each other are likely to have the same methylation status. A region with several unmethylated CpGs is a stronger biological indicator that is an actively expressed region than a single unmethylated CpG, and vice versa for methylated regions. This phenomenon in statistics is called spatial autocorrelation. A cluster of five or six adjacent CpGs each at p ≈ 10⁻⁴ when tested individually (not reaching statistical significance in an EWAS) is, collectively, a statistically significant region.

DMR tools

The choice of DMR tool is constrained by what upstream modeling you did.

Importantcomb-p works on summary statistics; DMRcate needs line-level data

comb-p (Pedersen et al. 2012) takes only a BED-style formatted file of per-CpG p-values and coordinates. It never sees the samples or line-level methylation or phenotype data. That means it works downstream of any EWAS model, however complex. If your analysis plan involves stratification and meta-analysis, comb-p is the tool to use. This is why the Snakemake pipeline uses comb-p.

DMRcate (Peters et al. 2015) takes line-level β/M values plus a single design formula and re-fits the model internally to compute a smoothed, region-level statistic. It is powerful and self-contained, but it cannot consume meta-analyzed results — there is no single design matrix that represents “stratify by sex, then inverse-variance meta-analyze.” If you used the stratified-meta workflow, DMRcate is not an option; if you fit one pooled model with sex as a covariate, it is.

Running comb-p on our EWAS results

comb-p’s pipeline command chains four steps: estimate the spatial autocorrelation function (ACF) of the p-values, use it to compute autocorrelation-adjusted regional p-values (the Stouffer–Liptak–Kechris correction), FDR-correct those, and finally stitch adjacent low-p CpGs into regions. First we write the BED (this is the pipeline’s make_bed.R). The p-values come straight from the chapter 06 BACON fit; the coordinates come from the same Zhou v8.1 release the chapter 03 mask came from, so the BED and the mask cannot drift apart:

NoteWhat this chapter needs

This section starts from data/06_ewas.rds, written by chapter 06. If you did not run that chapter, fetch the published checkpoint instead:

Terminal

./get_data.sh F_ewas_results

It also reads ../masks/EPIC.ordering.tsv.gz and ../masks/EPIC.hg38.coord.tsv.gz, the probe-order and hg38 coordinate tables from the same Zhou v8.1 release chapter 03 used. If you do not have those two files, the same command fetches them:

Terminal

./get_data.sh H_annotation

# This chunk executes: it writes the file comb-p consumes, so rendering the
# chapter regenerates the BED from exactly the p-values discussed above.
tt  <- as.data.table(readRDS("data/06_ewas.rds")$tt)
# The YAME probe order and the coordinate table are two halves of one file pair:
# row i of the coordinate table is the probe named on row i of the ordering
# file, which is why they are assembled side by side and not merged.
ord <- fread("../masks/EPIC.ordering.tsv.gz", select = "Probe_ID")  # YAME probe order
crd <- fread("../masks/EPIC.hg38.coord.tsv.gz")                     # same order, one row per probe
map <- data.table(probe = ord$Probe_ID,
                  CpG_chrm = crd$CpG_chrm, CpG_beg = crd$CpG_beg)

bed <- merge(tt, map, by = "probe")[
  # comb-p walks ordered positions along a chromosome, so keep the assembled
  # chromosomes only: a probe on a scaffold has no neighbors to borrow from
  grepl("^chr([0-9]+|X|Y)$", CpG_chrm),
  # comb-p's column contract: chrom/start/end, the p-value it smooths, and an
  # id it carries through untouched. The '#' on the first name is what makes
  # the header a comment line to the tool.
  .(`#chrom` = CpG_chrm, start = CpG_beg, end = CpG_beg + 2L,
    pvals = bacon.p, cpgid = probe)]
setorder(bed, `#chrom`, start)
fwrite(bed, "data/08_annotation/PTSD_ewas_annotated_results.bed", sep = "\t")

That writes 756,234 of the 756,251 tested CpGs — the 17 that fall on unplaced or alt contigs are dropped, because comb-p works on ordered positions along a chromosome and a scaffold has no neighbors to borrow from.

Then the region detection, with the pipeline’s default parameters (config.yml: min_pvalue = 1e-4, window_size = 200, region_filter = 0.05):

comb-p pipeline \
    --seed 1e-4 \         # a CpG must reach this p to seed a region
    --dist 200 \          # CpGs within 200 bp can extend a region
    --region-filter-p 0.05 \
    -p data/08_annotation/dmr/PTSD_dmr \
    data/08_annotation/PTSD_ewas_annotated_results.bed

That one command writes seven files, all sharing the -p prefix. Nothing prints the answer to screen, so it is worth knowing which file holds what:

file	contents
PTSD_dmr.args.txt	the exact invocation, the comb-p version, and the run date — keep this for your methods section
PTSD_dmr.acf.txt	the fitted autocorrelation function: lag_min, lag_max, correlation, N, p
PTSD_dmr.slk.bed.gz	every input CpG with its smoothed regional p: #chrom, start, end, p (your input p), region-p (SLK-smoothed)
PTSD_dmr.fdr.bed.gz	the same rows plus region-q, the Benjamini–Hochberg correction of region-p
PTSD_dmr.regions.bed.gz	the stitched regions, minimal and headerless: chrom, start, end, min_p, n_probes
PTSD_dmr.regions-p.bed.gz	the regions with statistics: #chrom, start, end, min_p, n_probes, z_p, z_sidak_p
PTSD_dmr.regions-t.bed	the regions that pass --region-filter-p — here the single region with z_sidak_p < 0.05

regions-p.bed.gz is the file to read: it lists every candidate region with both the uncorrected and the Šidák-corrected region p. regions-t.bed is that same table filtered down to what survived, which makes it the convenient endpoint for a pipeline. In this run the two files carry identical columns because our input BED has no effect-size column for comb-p to carry through.

The ACF step reports short-range autocorrelation in our p-values (correlation ≈ 0.054 out to ~71 bp, p ≈ 1.3 × 10⁻¹²⁹ over 202,078 CpG pairs). Methylation is spatially structured here, as expected. The SLK/FDR steps then propose three candidate regions, and one of them survives the Šidák correction for the number of regions tested:

region (hg38)	n CpGs	min_p	z_sidak_p
chr6:28,633,493–28,633,701	11	2.3 × 10⁻¹³	2.6 × 10⁻¹⁶
chr3:138,608,544–138,608,573	2	7.8 × 10⁻⁵	0.45
chr20:32,018,064–32,018,066	1	2.1 × 10⁻⁵	1.00

Those are five of the seven columns in regions-p.bed.gz — the three coordinate columns are collapsed into one here, and z_p is omitted (it is the Šidák-uncorrected region p, 3.6 × 10⁻²⁰ for the chr6 region). The surviving region sits at chr6:28.63 Mb, inside the MHC/HLA region.

Warningmin_p in the comb-p output is not a raw probe p-value

The region above is reported with min_p = 2.3 × 10⁻¹³ and z_sidak_p = 2.6 × 10⁻¹⁶. Despite the column name, that min_p is not the smallest input p-value in the region: it is the smallest Stouffer–Liptak–Kechris–smoothed, FDR-corrected regional q among the region’s CpGs (the raw SLK p at those same CpGs goes down to 6.0 × 10⁻¹⁹). Smoothing pools each CpG with its neighbors, so these values are orders of magnitude below anything in the input: the smallest raw BACON p among the eleven CpGs is 1.7 × 10⁻⁵, and the smallest anywhere in the BED is 3.0 × 10⁻⁸. Smoothed regional statistics and probe-level statistics are not on the same scale, and comparing them directly will make a region look far stronger than the data supporting it.

ImportantSet the seed threshold on purpose

--seed is the p-value a CpG must reach before comb-p will even start a region there; --dist then decides how far the region may extend. The seed is not a significance threshold for the result — the Šidák-corrected region p is — so setting it at genome-wide-significance levels defeats the point of doing region analysis at all. The default of 1 × 10⁻⁴ finds the MHC region above. If your DMR run comes back empty, you can adjust the seed p-value to something less stringent, but make sure you have a reason to justify the p-value you do use!

TipThe MHC is a blood-EWAS landmark

That the strongest regional signal falls in the MHC on chromosome 6. The MHC is the most gene-dense, most polymorphic, and most methylation-variable region in the genome, and it surfaces in a large fraction of blood EWAS. It is simultaneously the region where you should be most cautious about SNP-driven artifacts. We already removed common-SNP-affected sites, but it is possible that rare variants, structural variants, or undocumented genetic variation is driving this signal. Further causal inference testing with genetic or gene expression data can help clarify why we are seeing this signal in the MHC.

Figure 3: The detected DMR at chr6:28.63 Mb. Each stem is one CpG; the eleven inside the shaded window (blue) are the region comb-p reports. Every one sits far below the single-CpG genome-wide threshold (−log₁₀ p = 7.2), yet together — because they are adjacent and move in the same direction — they reach Šidák-corrected p = 2.6 × 10⁻¹⁶. This power gain is the entire reason to do region-level analysis.

ENmix::combp: the same method, but with R

comb-p is a Python command-line tool, which means an extra environment to install for one step of one chapter. If you would rather stay in R, ENmix::combp (Xu et al. 2021) implements the same Stouffer–Liptak–Kechris approach, and it takes the BED we just built with only a column rename. (The same package also offers ipdmr (Xu et al. 2020), the authors’ interval-p-value variant.)

NoteWhat this chapter needs

This section starts from data/08_annotation/PTSD_ewas_annotated_results.bed, written by the make-bed chunk above. If you did not run that chunk, fetch the published copy instead:

Terminal

./get_data.sh H_annotation

# ENmix is not in `envs/methyl.yml` and comb-p itself lives in `envs/combp.yml`
# (Python 2.7), so this cannot run in the render environment. Shown here; the
# DMR outputs it produces come from the H_annotation tier.
library(ENmix)

bed <- data.table::fread("data/08_annotation/PTSD_ewas_annotated_results.bed")
# ENmix expects its own column names, in this order
data.table::setnames(bed, c("chr", "start", "end", "p", "probe"))
bed <- as.data.frame(bed)          # ENmix wants a plain data.frame
bed$chr <- as.character(bed$chr)

combp(data = bed,
      dist.cutoff = 750,           # must exceed bin.size (see the warning below)
      bin.size    = 310,           # ACF bin width
      seed        = 1e-4,          # FDR threshold for seeding a region
      nCores      = 4,
      region_plot = FALSE, mht_plot = FALSE)
# → resu_combp.csv  (chr, start, end, p, fdr, sidak, nprobe, probe)

It runs in about a minute on our 756,234-row BED and finds the same lead region: a 12-CpG DMR at chr6:28,633,494–28,633,743, the MHC-proximal locus command-line comb-p reports as 11 CpGs at chr6:28,633,493–28,633,701. The probe lists overlap almost completely; ENmix extends the region by one CpG (cg12763978) at the 3′ end.

WarningTwo ways this differs from the command-line tool:

bin.size must be smaller than dist.cutoff. ENmix builds its autocorrelation bins by seq(bin.size, dist.cutoff, bin.size), so calling it with the shipped command-line settings (--dist 200 with the default bin.size = 310) fails with wrong sign in 'by' argument rather than a message naming the problem. The dist.cutoff = 750 above is chosen to exceed the default bin width.

The two implementations do not return the same region count. With the settings above ENmix reports 8 candidate regions of which 7 clear Šidák < 0.05, where command-line comb-p at --dist 200 --seed 1e-4 reports 3 candidates of which 1 clears Šidák. Both single out the same chr6 locus as by far the strongest — ENmix calls 12 CpGs there against comb-p’s 11 — but ENmix is much more permissive about everything else. Some of that is the wider dist.cutoff, and some is that ENmix describes itself as a modified comb-p. Treat it as an R-native implementation of the same idea. Either method is okay, just make sure to report which one you used and its parameters.

The DMR numbers in this chapter come from the command-line tool at the parameters shown above, which is what chapter 07’s pipeline invokes.

DMRcate: the line-level alternative

If you are using DMRcate, it operates on the M-value matrix directly:

RStudio Console

# Shown, not run: `M` and `design` are yours to supply. This chapter never
# builds them, because our EWAS was run stratified (see the callout above),
# so there is no single design matrix to hand DMRcate.
library(DMRcate)
# M = CpG x sample matrix; design = model.matrix(~ PTSD + sex + age + cellcounts + SVs)
myannotation <- cpg.annotate("array", M, arraytype = "EPICv1",
                             analysis.type = "differential",
                             design = design, coef = "PTSD")
dmrs <- dmrcate(myannotation, lambda = 1000, C = 2)   # Gaussian-kernel smoothing
results.ranges <- extractRanges(dmrs, genome = "hg38")

Note what it requires that comb-p does not: the full M matrix in memory and a single design/coef. If your EWAS was residualized or meta-analyzed, you cannot express it as one coef here.

Annotating the DMRs

Once regions exist, the pipeline’s dmr_annotation.R maps each DMR back to genes and islands by overlapping its coordinates against UCSC refGene, cpgIslandExt, and the HGNC BigBed (fetched by the annotate.smk cache rule from hgdownload.soe.ucsc.edu, all hg38). The DMR’s constituent CpGs are looked up in the EWAS BED so each region carries its probe list and direction of effect. The overlap is a coordinate join — bedtools intersect in the pipeline, or GenomicRanges::findOverlaps() if you are staying in R — with no modeling involved. We run it here for our one region, the eleven CpGs of the comb-p MHC region at chr6:28,633,493–28,633,701:

# eval: false because the two UCSC BEDs it reads are the `annotate.smk` cache
# downloads, which are not tracked in this repository.
library(GenomicRanges)   # GRanges(), IRanges(), overlapsAny()
library(rtracklayer)     # import(), for the cached UCSC BED files
# DMR CpGs (from the comb-p region) as a GRanges
dmr <- GRanges("chr6", IRanges(28633493, 28633701))

# refGene / HGNC / cpgIslandExt hg38 BEDs (annotate.smk cache; here read locally)
refGene <- import("refGene.bed.gz")          # rtracklayer::import
hgnc    <- import("hgnc.bed.gz")
# genes whose body/±5 kb overlaps the DMR
refGene[overlapsAny(refGene, dmr + 5000)]    # -> character(0)
hgnc[overlapsAny(hgnc, dmr)]                  # -> COX8CP1 (pseudogene)

This is a prime example region because the three annotation sources disagree, and the pipeline’s default gene table is the one that comes up empty. While often the sources do agree, it is not uncommon for annotation disagreements to occur:

# The comparison rendered below. The third row is computed from our own
# probe-level annotation; the first two are one-line readings of the two overlap
# queries in the chunk above, which is why they are written out as text — an
# empty overlap has no gene symbol to print.
# It needs annotated.rds, written in Part 1. The two UCSC overlap queries in
# the chunk above are not re-run here -- their one-line outcomes are the text in
# rows one and two.
ann <- as.data.table(readRDS("data/08_annotation/annotated.rds"))
# the region's CpGs by coordinate, not by a pasted probe list, so the region
# definition lives in one place
reg <- ann[CpG_chrm == "chr6" & CpG_beg >= 28633493 & CpG_beg <= 28633701]
setorder(reg, CpG_beg)
stopifnot(nrow(reg) == 11L)     # the eleven comb-p reported; fail loudly if the map changed

# every gencode v41 gene the eleven probes are assigned to, de-duplicated
gencode_genes <- paste(sort(unique(unlist(strsplit(reg$genesUniq, ";")))),
                       collapse = ", ")

dmr_src <- data.table(
  `Annotation source` = c("UCSC refGene (pipeline default gene table)",
                          "UCSC HGNC BigBed (annotate.smk cache)",
                          "Zhou GENCODE v41 (our probe-level annotation)"),
  `What it names for chr6:28,633,493–28,633,701` = c(
    # refGene[overlapsAny(refGene, dmr + 5000)] returned nothing at all
    "nothing within \u00b15 kb (nearest ZBED9 ~17 kb up, LINC00533 ~15 kb down)",
    # hgnc[overlapsAny(hgnc, dmr)] returned the pseudogene
    "COX8CP1 (COX8C pseudogene, directly overlapping)",
    paste0(gencode_genes, " (promoter-associated, non-coding)")))
fwrite(dmr_src, "data/08_annotation/dmr_source_comparison.csv")

# written by the `dmr-sources-derive` chunk above
knitr::kable(fread("data/08_annotation/dmr_source_comparison.csv"),
  caption = "Three hg38 gene annotations, three answers for the same DMR. refGene — the pipeline's default `gene_table` — names nothing nearby; the region only acquires a label from the pseudogene track or from GENCODE.")

Three hg38 gene annotations, three answers for the same DMR. refGene — the pipeline’s default gene_table — names nothing nearby; the region only acquires a label from the pseudogene track or from GENCODE.

Annotation source	What it names for chr6:28,633,493–28,633,701
UCSC refGene (pipeline default gene table)	nothing within ±5 kb (nearest ZBED9 ~17 kb up, LINC00533 ~15 kb down)
UCSC HGNC BigBed (annotate.smk cache)	COX8CP1 (COX8C pseudogene, directly overlapping)
Zhou GENCODE v41 (our probe-level annotation)	ENSG00000271440, ENSG00000287279 (promoter-associated, non-coding)

refGene is a curated set of well-supported transcripts, so a DMR sitting in a gene-poor stretch of the MHC returns nothing within 5 kb — the nearest curated genes (ZBED9, LINC00533) are >15 kb away. The HGNC BigBed, which includes pseudogenes and genes from difficult to map regions, overlaps COX8CP1 directly, and our probe-level GENCODE v41 annotation assigns the CpGs to two non-coding transcripts as promoter hits. None of these is “wrong”; they reflect different inclusion criteria. The lesson for a DMR walkthrough is to read the coordinates, not just the gene symbol: an empty gene_table overlap does not mean the region is intergenic junk — here it is a CpG-island shore over a pseudogene in an immune-gene-dense locus, exactly the kind of place methylation differences cluster. Our eleven CpGs, all hypomethylated in cases, look like this:

# The eleven-row table rendered below, and the fuller per-CpG file that backs
# it. Both come from `reg` — the region's CpGs in the Part 1 annotation layer,
# built in the chunk above.
disp <- reg[, .(CpG = probe,
                `hg38 position (chr6)` = formatC(CpG_beg, big.mark = ",", format = "d"),
                # the annotation layer's label is "Promoter (±1.5 kb TSS)";
                # trim the parenthetical so the column stays narrow
                `GENCODE feature` = sub(" \\(.*", "", feature),
                `CGI context` = island,
                # β difference on the 0-1 scale: -0.022 reads as 2.2 points hypo
                `Δβ (case − control)` = round(delta_beta, 3),
                `BACON p` = signif(bacon.p, 3))]
fwrite(disp, "data/08_annotation/dmr_cpg_display.csv")

# the same eleven CpGs with everything the annotation layer knows about them,
# including the transcript list the display table has no room for
full <- reg[, .(probe, CpG_beg, CpG_end, genesUniq, transcriptTypes, distTSS,
                feature, island, BIOS_eQTM_genes,
                delta_beta = round(delta_beta, 4),
                bacon.es   = round(bacon.es, 3),
                bacon.p    = signif(bacon.p, 3))]
fwrite(full, "data/08_annotation/dmr_cpg_annotation.csv")

# written by the `dmr-cpgs-derive` chunk above
knitr::kable(fread("data/08_annotation/dmr_cpg_display.csv"),
  caption = "The eleven CpGs of the DMR. All fall in the same CpG-island north shore and are hypomethylated in PTSD cases — Δβ between −2.2 and −5.6 percentage points, every one in the same direction. That coherence is what makes a region call meaningful; none of the eleven individual BACON p-values comes close to genome-wide significance.")

The eleven CpGs of the DMR. All fall in the same CpG-island north shore and are hypomethylated in PTSD cases — Δβ between −2.2 and −5.6 percentage points, every one in the same direction. That coherence is what makes a region call meaningful; none of the eleven individual BACON p-values comes close to genome-wide significance.

CpG	hg38 position (chr6)	GENCODE feature	CGI context	Δβ (case − control)	BACON p
cg07017437	28,633,493	Promoter	N_Shore	-0.022	3.00e-04
cg25653641	28,633,534	Promoter	N_Shore	-0.056	1.69e-05
cg22572476	28,633,546	Promoter	N_Shore	-0.056	2.15e-04
cg19488431	28,633,551	Promoter	N_Shore	-0.053	3.23e-05
cg27535677	28,633,587	Promoter	N_Shore	-0.035	1.45e-03
cg22497095	28,633,597	Promoter	N_Shore	-0.041	7.45e-04
cg26865747	28,633,599	Promoter	N_Shore	-0.042	2.08e-04
cg03759229	28,633,639	Promoter	N_Shore	-0.038	2.32e-03
cg13565129	28,633,641	Promoter	N_Shore	-0.041	3.89e-04
cg14654363	28,633,665	Promoter	N_Shore	-0.026	3.09e-03
cg00990380	28,633,699	Promoter	N_Shore	-0.024	1.23e-03

From regions to pathways: goregion

The gene-set tools in Part 3 take a list of CpGs. When your unit of discovery is a region, missMethyl::goregion (Phipson et al. 2016) is the region-level analog of gometh: give it a GRanges of your DMRs and it maps each region to the CpGs it contains, collapses those to genes, and runs the same probe-number-bias-corrected GO/KEGG test. It is the natural way to ask “are my comb-p regions collectively enriched for a biological process?”

# error: true, not eval: false -- this call stops with an error on purpose, and
# the point of the section is to show that error, so we let it run and let
# Quarto print what R says. Needs the missMethyl stack and the hg19 EPIC
# annotation package.
library(missMethyl)
# goregion works in the array's native annotation build (hg19 for EPIC),
# so lift the hg38 DMR coordinates back with the EPIC manifest's own map.
anno <- minfi::getAnnotation("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
# the eleven probes the region table above derived, spelled out so this chunk
# stands on its own
dmr_cpgs <- c("cg07017437","cg25653641","cg22572476","cg19488431",
              "cg27535677","cg22497095","cg26865747","cg03759229",
              "cg13565129","cg14654363","cg00990380")
reg <- GenomicRanges::reduce(GenomicRanges::GRanges(
  anno[dmr_cpgs, "chr"],
  IRanges::IRanges(anno[dmr_cpgs, "pos"], width = 2)))

# all.cpg is the tested universe: `tt` here is the chapter 06 top-table read in
# by the make-bed chunk above
gr_go <- goregion(reg, all.cpg = tt$probe, collection = "GO",
                  array.type = "EPIC")

Error in `getMappedEntrezIDs()`:
! There are no genes annotated to the significant CpGs

Run on our single MHC DMR, goregion returns an error: “There are no genes annotated to the significant CpGs.”. The region sits over a pseudogene (COX8CP1) with no protein-coding gene in the EPIC annotation, so there is nothing to map to a GO term. Two lessons fall out of this:

1. Enrichment inherits the annotation’s gene assignments. If the annotation maps a region to no gene, a gene-centric lens for pathway enrichment will not have anything to query or report.
2. goregion needs many regions to be meaningful. One region is never an “enrichment”; a real DMR analysis feeds it the dozens-to-hundreds of regions comb-p or DMRcate returns genome-wide. We show the mechanism here even though it doesn’t produce a meaningful result so the call is in your toolkit.

Part 3 — Pathway, trait, and functional enrichment

Annotation names the genes near your hits; enrichment asks whether that gene list is collectively over-represented in any biological process, pathway, or previously reported trait. This is where an EWAS stops being a list of CpGs and starts pointing at mechanism. Two things make methylation enrichment different from the gene-expression version you may know:

1. Probe-number bias. Genes are not covered evenly on the array — a large gene like PRKG1 carries hundreds of probes, a small one only a handful. A naive test that treats “gene has ≥1 significant CpG” as a hit will therefore flag long, densely-tiled genes purely because they had more chances to be picked. Enrichment for methylation must correct for the number of probes per gene.
2. Multiple CpGs map to one gene, and CpGs in intergenic space map to none. The tool has to collapse probes to genes sensibly before testing.

missMethyl gometh: the reproducible, bias-corrected default

missMethyl::gometh() (and its generic gsameth() for custom gene sets) is the field-standard R implementation. It maps your significant CpGs to genes, models the per-gene probe count as the bias covariate, and runs a corrected hypergeometric test against GO or KEGG (Phipson et al. 2016). Enrichment has not yet been added to the Snakemake pipeline, but will in a future release.

# The Wallenius test runs over every GO term, so this is slow, and it needs the
# missMethyl stack plus the hg19 EPIC annotation package. The RDS written at the
# end is what the comparison table further down reads.
library(missMethyl)
tt     <- data.table::fread("data/06_ewas_bacon_toptable.csv.gz")
all.cpg <- tt$probe                          # tested universe = every CpG you analyzed
sig.cpg <- tt[order(bacon.p)][1:1000]$probe  # your significant set

# GO enrichment, corrected for probes-per-gene. array.type matches your array.
go <- gometh(sig.cpg = sig.cpg, all.cpg = all.cpg,
             collection = "GO", array.type = "EPIC", plot.bias = TRUE)

go <- go[order(go$P.DE), ]

# KEGG is a one-word change
kegg <- gometh(sig.cpg = sig.cpg, all.cpg = all.cpg,
               collection = "KEGG", array.type = "EPIC")

# Both collections in one object. The method comparison below reads `kegg` back
# out of it, and keeping `go` alongside means the GO panel of the figure can be
# redrawn without re-running the test.
saveRDS(list(go = go, kegg = kegg), "data/08_annotation/08_gometh.rds")

Two arguments carry all the weight. all.cpg is the full tested set — it must be every CpG that survived your QC and filtering, not the whole array manifest, or the test is biased toward whatever your pipeline happened to remove. plot.bias = TRUE draws the probe-number bias curve; glance at it once to confirm the correction is actually engaging. Defining “significant” by a fixed top-N (here 1,000) is a pragmatic choice when nothing clears FDR; with real signal you would instead pass the CpGs at your chosen significance threshold.

Run on our teaching subset, the corrected results are:

Figure 4: GO and KEGG enrichment of the top 1,000 CpGs by BACON p, computed with gometh (probe-number–corrected hypergeometric test). Bars show nominal \(-\log_{10} P_{DE}\); labels are DE-genes / gene-set size. Nothing survives FDR at n = 87 — the expected outcome when no CpG reaches epigenome-wide significance. On a real study these panels are where mechanism appears.

The top nominal terms here are a grab-bag — basement membrane assembly (GO:0070831, P = 1.8 × 10⁻⁵), extracellular matrix assembly, viral translational termination-reinitiation, regulatory T cell differentiation — with no coherent theme and no term at FDR < 0.05 in either collection. The smallest GO FDR is 0.41 and the smallest KEGG FDR is 0.38 (Basal cell carcinoma, nominal P = 0.0010). This is exactly what enrichment of a null gene list should look like. Reading this result as “PTSD methylation remodels the extracellular matrix” would be false: with 1,000 arbitrarily-chosen CpGs, a handful of small gene sets always reach nominal p < 0.01 by chance — note that the leading GO term has only 15 genes in it, 7 of which the top set touched. The FDR column, not the nominal p, is the result.

methylGSA: the same question, a different bias model

gometh is not the only way to correct for the probe-number bias, and the choice of correction can change the outcome. methylGSA (Ren and Kuan 2019) offers two alternatives to gometh’s Wallenius hypergeometric test:

• methylglm models each gene set with a logistic regression in which the per-gene number of probes enters as a covariate — a continuous adjustment rather than a weighted urn, and it takes the whole ranked p-value vector, not a hard top-N cutoff.
• methylRRA aggregates the CpG p-values to a per-gene score by robust rank aggregation, then runs either an over-representation test or a GSEA-style rank test — no significance cutoff at all.

NoteThis section runs in its own environment

On the R-only route, BiocManager installs methylGSA alongside everything else and there is nothing to do. Conda users need envs/methylgsa.yml (conda env create -f envs/methylgsa.yml) because bioconda’s newest methylGSA build is for an older R than the rest of the stack. Nothing is shared at run time: the only input below is the BACON top-table CSV written by chapter 06.

# This chunk needs `envs/methylgsa.yml`, a separate R 4.2 environment -- bioconda
# has no r45 build of methylGSA, and pinning the whole tutorial back to R 4.2 for
# one section is the trade this project deliberately refuses. So it cannot run in
# the render environment and is shown instead; its two output tables are
# committed and read below.
library(methylGSA)
library(data.table)   # fwrite(), for the two top-tables written at the end
# methylGSA needs the array annotation LOADED, not merely installed --
# methylglm/methylRRA stop with an error if it is missing from the search path.
library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)

tt <- data.table::fread("data/06_ewas_bacon_toptable.csv.gz")
pv <- setNames(tt$bacon.p, tt$probe)          # NAMED vector: names = CpG IDs

# logistic-regression enrichment (bias covariate = probes per gene)
mglm <- methylglm(cpg.pval = pv, array.type = "EPIC",
                  GS.type = "KEGG", minsize = 10, maxsize = 500)

# GSEA on rank-aggregated per-gene scores
mrra <- methylRRA(cpg.pval = pv, array.type = "EPIC", method = "GSEA",
                  GS.type = "KEGG", minsize = 10, maxsize = 500)

# --- what leaves this environment ----------------------------------------
# Both fits in one object, plus the two flat top-tables the comparison below
# reads. `Method` is carried in each file so the three tools' tables can be
# row-bound, and the KEGG id keeps its leading zeros by staying a string.
saveRDS(list(methylglm = mglm, methylRRA = mrra),
        "data/08_annotation/08_methylgsa.rds")
glm_top <- head(as.data.table(mglm)[order(pvalue)], 10)[
  , .(Method = "methylglm", KEGG = ID, Pathway = Description, Size, pvalue, padj)]
fwrite(glm_top, "data/08_annotation/08_methylgsa_glm_top.csv")
rra_top <- head(as.data.table(mrra)[order(pvalue)], 5)[
  , .(Method = "methylRRA-GSEA", KEGG = ID, Pathway = Description, Size,
      NES, pvalue, padj)]
fwrite(rra_top, "data/08_annotation/08_methylgsa_rra_top.csv")

We ran all three tools on the same BACON top-table, against the same KEGG collection, and recorded both how many pathways each calls significant at FDR < 0.05 and what each one puts at the top:

# The comparison table rendered below, rebuilt from the three result files the
# two chunks above write rather than typed out by hand.
g  <- readRDS("data/08_annotation/08_gometh.rds")
# gometh returns the KEGG id in the row names, which fread() cannot see later
kg <- as.data.table(g$kegg, keep.rownames = "ID")[order(P.DE)]
# KEGG ids are zero-padded ("04140"), so they must not be read as integers. Note
# the two files disagree on what they call things: methylglm writes ID and
# Description, methylRRA-GSEA writes KEGG and Pathway. Reading the wrong name
# gives NULL, which c() silently drops -- and a three-row table built from a
# two-element vector recycles, pairing each method with another method's pathway
# while the p-values stay correctly aligned. Name the columns each file uses.
gl <- fread("data/08_annotation/08_methylgsa_glm_top.csv", colClasses = c(ID = "character"))
rr <- fread("data/08_annotation/08_methylgsa_rra_top.csv", colClasses = c(KEGG = "character"))

cmp <- data.table(
  Method = c("gometh (Wallenius hypergeometric)",
             "methylglm (logistic regression)",
             "methylRRA-GSEA (rank aggregation)"),
  Input  = c("top 1,000 CpGs vs universe",
             "full named p-value vector",
             "full named p-value vector"),
  `Bias model` = c("probes per gene (Wallenius)",
                   "probe count as covariate",
                   "gene-level RRA score"),
  # the two methylGSA files hold their own top rows ordered by p, so a zero here
  # is a zero over the whole collection: the smallest FDR is in the first row
  `KEGG FDR<0.05` = c(sum(kg$FDR < 0.05), sum(gl$padj < 0.05), sum(rr$padj < 0.05)),
  `Top KEGG pathway` = c(kg$Description[1], gl$Description[1], rr$Pathway[1]),
  `Its nominal p` = signif(c(kg$P.DE[1], gl$pvalue[1], rr$pvalue[1]), 3),
  `Its FDR`       = signif(c(kg$FDR[1],  gl$padj[1],   rr$padj[1]),   3)
)
fwrite(cmp, "data/08_annotation/08_gsa_method_comparison.csv")

# written by the `gsa-compare-derive` chunk above
knitr::kable(fread("data/08_annotation/08_gsa_method_comparison.csv"),
  caption = "Three probe-bias corrections, one dataset (top: gometh from the section above; methylGSA's two models below). All three return zero KEGG pathways at FDR < 0.05, but each nominates a different pathway as its leader — the ranking, not just the count, depends on the bias model and the input convention.")

Three probe-bias corrections, one dataset (top: gometh from the section above; methylGSA’s two models below). All three return zero KEGG pathways at FDR < 0.05, but each nominates a different pathway as its leader — the ranking, not just the count, depends on the bias model and the input convention.

Method	Input	Bias model	KEGG FDR<0.05	Top KEGG pathway	Its nominal p	Its FDR
gometh (Wallenius hypergeometric)	top 1,000 CpGs vs universe	probes per gene (Wallenius)	0	Basal cell carcinoma	0.001010	0.3770
methylglm (logistic regression)	full named p-value vector	probe count as covariate	0	Phosphatidylinositol signaling system	0.004710	0.4480
methylRRA-GSEA (rank aggregation)	full named p-value vector	gene-level RRA score	0	Hedgehog signaling pathway	0.000463	0.0976

On this dataset the three agree on the count — zero KEGG pathways at FDR < 0.05 from all three — and disagree on everything else. Each nominates a different leader (Basal cell carcinoma, Phosphatidylinositol signaling, Hedgehog signaling), and the smallest FDR any of them reaches spans an order of magnitude: 0.38 for gometh, 0.45 for methylglm, 0.098 for methylRRA-GSEA. That last number is the informative one. methylRRA’s GSEA mode uses the full ranking and no cutoff, so it is the most sensitive and the most prone to calling pathways in a null dataset — it comes within a factor of two of significance here, on data where nothing is there to find. Had our nominal p been a little smaller, or the gene set a little larger, it alone would have produced a “significant pathway” to write up.

# written by the `methylgsa` chunk above
rra <- fread("data/08_annotation/08_methylgsa_rra_top.csv")
# NES is only present when the file came from a GSEA run, so take the columns
# that are there instead of failing the render on a missing one
rra_cols <- intersect(c("KEGG", "Pathway", "Size", "NES", "pvalue", "padj"), names(rra))
knitr::kable(rra[, ..rra_cols],
  digits = 4,
  caption = "The five top-ranked methylRRA-GSEA KEGG pathways. None survives FDR < 0.05, but Hedgehog signaling reaches FDR = 0.098 on a dataset with no epigenome-wide signal — the sensitivity/specificity trade-off of a cutoff-free rank test, made concrete.")

The five top-ranked methylRRA-GSEA KEGG pathways. None survives FDR < 0.05, but Hedgehog signaling reaches FDR = 0.098 on a dataset with no epigenome-wide signal — the sensitivity/specificity trade-off of a cutoff-free rank test, made concrete.

KEGG	Pathway	Size	NES	pvalue	padj
4340	Hedgehog signaling pathway	55	1.2015	0.0005	0.0976
910	Nitrogen metabolism	21	1.2071	0.0105	0.5180
620	Pyruvate metabolism	39	1.1511	0.0146	0.5180
5217	Basal cell carcinoma	55	1.1316	0.0078	0.5180
4974	Protein digestion and absorption	74	1.1198	0.0052	0.5180

ImportantThe tool is a modeling choice, not a lookup

Report the method and its settings, not just “KEGG enrichment.” A reviewer who reruns your CpGs through a different tool will get a different list, and that is expected — gometh, methylglm, and methylRRA encode different assumptions about how probe density biases the test. Pick one a priori (we default to gometh for its conservative, well-validated behavior), and if you report a second for sensitivity, say so and show both. Never run all three and present the one with the prettiest pathway.

KnowYourCG (KYCG): CpG-centric enrichment on regulatory elements

gometh asks a gene-set question — it collapses your CpGs to genes and tests GO/KEGG. But methylation acts through regulatory elements — CpG islands, transcription-factor binding sites, chromatin states, histone marks — and collapsing to genes throws that resolution away. KnowYourCG (KYCG), part of Zhou’s sesame/knowYourCG ecosystem (Goldberg et al. 2025), keeps the CpG as the unit of analysis and tests your hit set directly against curated databases of these elements. It is the CpG-centric complement to gometh, not a replacement, and it comes from the same lab as the annotation we used in Part 1, so the coordinates and databases are build-consistent by construction.

# knowYourCG is not in `envs/methyl.yml`: it needs Bioconductor 3.18+ and the
# environment is pinned at r-base=4.5/Bioc 3.22 for the rest of the stack, so
# this call is shown rather than run. Its output is the committed
# data/08_annotation/08_kycg_results.csv that the next three chunks read, so
# every number below is derived from a real run of exactly this code.
library(knowYourCG)
library(sesameData)
tt      <- data.table::fread("data/06_ewas_bacon_toptable.csv.gz")
all.cpg <- tt$probe                              # universe = every CpG you tested
sig.cpg <- tt[order(bacon.p)][1:1000]$probe   # your query set

# What CpG-feature databases are available for this platform?
# ALWAYS list first: the catalog names carry datestamps and they change
# between releases. A wrong name fails at load time, not at test time.
listDBGroups("EPIC")           # EPIC knowledgebases: CGI, TFBS, chromHMM, HM, ...

# Cache the databases you want from ExperimentHub, then test.
dbs <- c("KYCG.EPIC.CGI.20210713",               # CpG-island relation
         "KYCG.EPIC.chromHMM.20211020",          # chromatin states
         "KYCG.EPIC.TFBSconsensus.20211013",     # TF binding sites
         "KYCG.EPIC.HMconsensus.20211013")       # histone marks
res <- testEnrichment(query = sig.cpg, databases = dbs,
                      universe = all.cpg, platform = "EPIC")
res <- res[order(res$FDR), ]

# All 863 tested features, one row each: the deposited result table, and the
# input to every KYCG number quoted below.
data.table::fwrite(res, "data/08_annotation/08_kycg_results.csv")

testEnrichment runs a Fisher test of your query CpGs against each database, using your tested set as the universe (the same universe discipline that matters for gometh). The databases are versioned knowledgebases pulled from ExperimentHub, so the run is reproducible.

The run confirms the universe it is testing against and how many features it tested: universe: 756251 | query: 1000, 863 features across the four knowledgebases. That number matters for the interpretation below — the four databases are wildly different sizes.

On our teaching subset KYCG returns 41 features at FDR < 0.05 where gometh returned nothing. That looks like a decisive win for the CpG-centric lens. It is partly real and partly an artifact of the test’s shape, and separating the two is the lesson of this section.

Figure 5: The twelve most significant KYCG features for the top 1,000 CpGs. Four EPIC knowledgebases were tested; only three are represented here, and the CGI knowledgebase contributes no significant feature at all. Bars are \(\log_2\) odds ratio; every bar shown clears FDR < 0.05. Note what is absent: no CpG-island-relation feature appears, and every significant enrichment is positive and modest (log₂OR 0.36–1.28).

# This chunk executes: it reads the result table the chunk above wrote and
# derives the three summaries the rest of this section discusses, so none of
# those numbers can drift away from the run that produced them.
k <- fread("data/08_annotation/08_kycg_results.csv")
# testEnrichment names the knowledgebase in `group` and the feature in `dbname`,
# which reads backwards; rename once here
setnames(k, c("group", "dbname"), c("db", "feature"))
# drop the "KYCG.EPIC." prefix and the release datestamp: "CGI", "chromHMM", ...
k[, kb := sub("^KYCG\\.EPIC\\.", "", sub("\\.[0-9]+$", "", db))]

# per-knowledgebase: how much was tested, how much survived, and the leader
ksum <- k[, .(
  `Features tested` = .N,
  `Nominal p<0.05`  = sum(p.value < 0.05, na.rm = TRUE),
  `FDR<0.05`        = sum(FDR < 0.05, na.rm = TRUE),
  `Top feature`     = feature[which.min(FDR)],
  `Top log2(OR)`    = round(estimate[which.min(FDR)], 2),
  `Top FDR`         = signif(min(FDR, na.rm = TRUE), 3)
), by = .(Knowledgebase = kb)][order(-`FDR<0.05`)]
fwrite(ksum, "data/08_annotation/08_kycg_db_summary.csv")

# every feature that clears FDR < 0.05, best first
ktop <- k[FDR < 0.05][order(FDR, p.value),
  .(Knowledgebase = kb, Feature = feature,
    `log2(OR)` = round(estimate, 2), `Query CpGs` = overlap,
    `DB size` = nD, p = signif(p.value, 3), FDR = signif(FDR, 3))]
fwrite(ktop, "data/08_annotation/08_kycg_top.csv")

# the CGI knowledgebase in full — all five rows, significant or not, because the
# section below reads it as a negative result and needs the near-misses visible
kcgi <- k[kb == "CGI"][order(FDR),
  .(Feature = feature, `log2(OR)` = round(estimate, 2),
    `Query CpGs` = overlap, `DB size` = nD,
    p = signif(p.value, 3), FDR = signif(FDR, 3))]
fwrite(kcgi, "data/08_annotation/08_kycg_cgi.csv")

# written by the `kycg-derive` chunk above
knitr::kable(fread("data/08_annotation/08_kycg_db_summary.csv"),
  caption = "KYCG results per knowledgebase. The four databases differ in size by two orders of magnitude, and almost all of the significant hits come from the largest one.")

KYCG results per knowledgebase. The four databases differ in size by two orders of magnitude, and almost all of the significant hits come from the largest one.

Knowledgebase	Features tested	Nominal p<0.05	FDR<0.05	Top feature	Top log2(OR)	Top FDR
TFBSconsensus	783	211	35	ING2	0.55	0.00328
HMconsensus	60	13	5	H3K36me2	0.95	0.01440
chromHMM	15	1	1	1_TssA	0.49	0.00328
CGI	5	1	0	Island	0.27	0.08350

Read that table before reading the hits. Three things fall out of it:

• The 41 hits are not spread evenly — 35 of 41 come from TFBSconsensus alone, which contributes 783 of the 863 features tested. Five come from histone marks, one from chromatin state.
• 1_TssA (active TSS) is the single chromatin-state hit (log₂OR = 0.49, FDR = 0.0033) — consistent with the mild promoter lean Figure 1 showed descriptively, and the one result here that lines up with an independent observation elsewhere in the chapter.
• CpG-island relation returns nothing at FDR < 0.05. Islands are enriched in the expected direction (log₂OR = 0.27, nominal p = 0.0078) but land at FDR = 0.084. The visual island lean in Figure 2 (22.3% vs 19.2%) is therefore not corroborated by the formal test — the descriptive impression survives, the inference does not.

# written by the `kycg-derive` chunk above
knitr::kable(head(fread("data/08_annotation/08_kycg_top.csv"), 15),
  caption = "The 15 most significant of KYCG's 41 features at FDR < 0.05. `Query CpGs` is how many of the 1,000 fell in that feature; `DB size` is how many of the 756,251 tested CpGs the feature covers.")

The 15 most significant of KYCG’s 41 features at FDR < 0.05. Query CpGs is how many of the 1,000 fell in that feature; DB size is how many of the 756,251 tested CpGs the feature covers.

Knowledgebase	Feature	log2(OR)	Query CpGs	DB size	p	FDR
chromHMM	1_TssA	0.49	229	131945	6.70e-06	0.00328
TFBSconsensus	ING2	0.55	170	92693	7.60e-06	0.00328
TFBSconsensus	RBL1	0.49	192	109245	2.28e-05	0.00553
TFBSconsensus	ING5	0.50	180	101596	2.82e-05	0.00553
TFBSconsensus	ZNF547	1.10	37	13355	3.20e-05	0.00553
TFBSconsensus	UHRF2	0.64	93	46736	7.46e-05	0.01070
TFBSconsensus	HMGN3	0.55	121	65192	1.16e-04	0.01310
TFBSconsensus	TRAF7	0.65	85	42409	1.21e-04	0.01310
TFBSconsensus	RAG2	0.43	196	115584	1.40e-04	0.01340
HMconsensus	H3K36me2	0.95	39	15602	1.67e-04	0.01440
HMconsensus	H3K4me3B	0.54	114	61523	1.97e-04	0.01540
TFBSconsensus	ZNF600	0.46	151	86776	3.13e-04	0.02250
TFBSconsensus	ZNF146	1.10	25	8960	5.22e-04	0.03230
TFBSconsensus	SMAD2/3	0.87	38	16020	5.37e-04	0.03230
TFBSconsensus	FAM208A	0.51	107	58729	5.61e-04	0.03230

Warning41 significant features is not 41 findings

Before you write these up, check three things that the FDR column alone will not tell you.

1. Count the nominal hits. 226 of 863 features reach nominal p < 0.05, against ~43 expected under the null. That excess is real — the top-1,000 set is genuinely non-random with respect to regulatory annotation — but it is distributed across hundreds of overlapping features, not concentrated in a few.

2. The features are not independent. The TFBSconsensus databases are consensus binding sites derived from ChIP-seq, and they overlap each other heavily: the median significant feature covers 7.4% of the tested universe (56,061 of 756,251 CpGs), and the largest covers 17%. ING2, RBL1, ING5, RAG2, and ZNF335 each overlap 90,000–120,000 CpGs. Benjamini–Hochberg assumes independence (or positive dependence); with sets this large and this correlated, “41 significant” is closer to a handful of distinct signals counted repeatedly. The more honest description is one broad enrichment in open, promoter-proximal, TF-bound chromatin, not 41 transcription factors.

3. Big databases win by construction. Every significant feature here has a positive log₂OR between 0.36 and 1.28 — small effects that only clear FDR because the overlap counts are large enough to make the Fisher test precise. A feature covering 5,000 CpGs needs a much bigger odds ratio to reach the same p. That is why TFBSconsensus, with the most and largest sets, supplies 85% of the hits, and why a database of five island categories supplies none.

Both tests see the same underlying fact — the top CpGs lean toward open, promoter-proximal chromatin — but they differ in what they can resolve. gometh has to collapse CpGs to genes and then test gene sets of a few dozen to a few hundred members, which at this sample size leaves it no power at all. KYCG keeps the CpG as the unit and tests sets of tens of thousands, so the same weak lean becomes formally detectable. The price is resolution: KYCG tells you which kind of chromatin, not which pathway, and its overlapping sets make the hit count a poor guide to how many independent things you have found.

Run both. Report gometh’s null as a null, and report KYCG’s 41 features as what they are — a single, modest, broadly-distributed enrichment in regulatory chromatin, at a sample size where no individual CpG is convincing.

EWAS Atlas / EWAS Open Platform Toolkit: the copy-paste web route

When you want a fast look without writing R — or want to cross-reference against what other studies have already reported — the EWAS Toolkit, part of the EWAS Open Platform (formerly EWAS Atlas) at ngdc.cncb.ac.cn/ewas/toolkit, takes a pasted list of CpG probe IDs and runs seven enrichment modalities in one pass (Xiong et al. 2022; Li et al. 2019). The value is that several of these draw on curated resources you can’t easily reproduce locally:

• Trait enrichment — tests your CpG list against the platform’s manually curated catalog of associations from thousands of published EWAS, using a weighted Fisher’s exact test where each probe is weighted by the number of studies reporting that probe–trait association. This is the analysis you can’t reproduce locally without their database, and it is the most useful: it tells you whether your hits overlap CpGs already linked to smoking, aging, BMI, or a specific disease.
• Genomic location enrichment — assigns each probe to 13 location categories defined relative to gene and CpG island (each probe gets one gene-relative and one island-relative category) and tests over-representation with a hypergeometric test (p-value + odds ratio). This overlaps what knowYourCG and our own feature/island tables already give you — but note the Toolkit’s location annotation is hg19-based (450K: IlluminaHumanMethylation450kanno.ilmn12.hg19; EPIC: IlluminaHumanMethylationEPICanno.ilm10b4.hg19), so it is not build-consistent with the hg38 annotation we built in Part 1.
• GO & KEGG enrichment — implemented on top of missMethyl::gometh, so it is literally the same probe-bias-corrected test you ran locally; expect concordance with your gometh result, and read a disagreement as a gene-universe difference (the web tool uses the full array as background).
• Chromatin state & histone modification enrichment — overlaps your CpGs against Roadmap Epigenomics data (chromatin state via the core 15-state model; 11 histone/accessibility marks: DNase, H2A, H3K27ac, H3K27me3, H3K4me1, H3K4me2, H3K4me3, H3K79me2, H3K9ac, H3K9me3, H4K20me1), again by hypergeometric test. This is the same question knowYourCG’s chromHMM/histone databases and eFORGE answer — a regulatory-context check — but built on hg19 Roadmap tracks.
• Tissue methylation — reports whether your CpGs are tissue-specific across 31 core tissues from the EWAS Data Hub (a probe counts as tissue-specific when abs(tau_hyper − tau_hypo) > 0.7, a cutoff set by permutation). Useful as a sanity check that whole-blood hits aren’t dominated by a signal that is really tissue-of-origin driven.
• Expression regulation — reports Pearson correlations between methylation and expression across six tissues (brain, colon, kidney, liver, stomach, testis). This is the same idea as the BIOS eQTM layer we joined in Part 1, but computed in solid tissues rather than whole blood — so for a blood EWAS, prefer the BIOS blood eQTMs and read this panel as supporting, cross-tissue context.
• Motif enrichment — runs HOMER on a 500 bp window centered on each CpG to find over-represented transcription-factor binding motifs, a lead on which TFs might read or be blocked by the methylation change.

The workflow is literally copy-paste: take the probe column of your top hits (or the DMR constituent CpGs), paste into the query box, submit. Practical notes:

• It is EPIC/450K-aware but check the array setting, and paste cg######### IDs, not coordinates.
• Mind the genome build. The location, chromatin-state, and histone panels are annotated on hg19; the rest of this tutorial is deliberately hg38. That is fine for the categorical questions these panels answer (is my set enriched for enhancers? for island shores?), but don’t copy an hg19 coordinate from the Toolkit into an hg38 result table.
• Treat it as a curated cross-check, not a replacement for the reproducible R path. Its trait and tissue panels are the parts you genuinely can’t run locally; its GO/KEGG and location panels duplicate gometh/knowYourCG, so use those in the pipeline for the record and the Toolkit interactively to mine the catalogs.

eFORGE: is a hit list driven by cell composition?

For a whole-blood EWAS the single most important adversarial check is whether your top CpGs are simply marking a shift in cell-type proportions rather than a phenotype-specific signal. eFORGE (Breeze et al. 2016) answers this directly: it takes a set of CpGs (your top hits, or a DMR’s constituent probes) and tests — by matched-background resampling — whether they are over-represented in cell-type-specific regulatory elements: DNase I hypersensitive sites and the Roadmap/ENCODE core-15 chromatin states across ~40 blood and other primary cell types.

The result you want to see for a well-adjusted blood EWAS is no strong, single-cell-type enrichment — that is evidence your cell-composition adjustment (chapter 04) did its job. The result that should stop you is a list of neutrophil- or T-cell-specific enhancers: a signature that the “phenotype” effect is really a composition effect that leaked past the model. Because it works on cg######### IDs against pre-built reference panels, eFORGE runs from a web form or the eForge R package with no local reference data to assemble.

NoteWhy this check matters here specifically

Our subset carries the confounds chapter 05 mapped — cell composition loads heavily on the leading PCs. An eFORGE cell-type enrichment on the top hits is the natural closing sanity check: it asks the composition question of the final CpG list, after all modeling, in the vocabulary (chromatin states, DHS) that a reviewer of a blood EWAS will expect to see addressed.

Other tools worth knowing, and how to choose

No single tool answers every interpretation question. A short map of the landscape, roughly in the order you’d reach for them:

Question	Tool	Notes
GO / KEGG, bias-corrected, reproducible	missMethyl::gometh / gsameth	Default gene-set test. gsameth takes any custom gene-set list (e.g. MSigDB Hallmark).
CpG-centric enrichment on regulatory elements	knowYourCG::testEnrichment	CGI, TFBS, chromatin state, histone marks — tests CpGs directly, no gene collapse. Build-consistent with the Zhou annotation.
Region-level GO enrichment from DMRs	missMethyl::goregion	Feeds comb-p/DMRcate regions in directly, correcting for probes-per-region.
“Has this CpG been linked to a trait before?”	EWAS Toolkit / EWAS Catalog / MRC-IEU EWAS Catalog	Curated look-ups; the fastest external validation.
Broad gene-set / pathway diagrams, no bias model	ShinyGO, g:Profiler, Enrichr, GREAT	Built for gene lists generally; do not correct for methylation probe bias, so feed them a gene list you already trust, not raw CpG hits.
Cell-type / tissue enrichment (is a hit driven by composition?)	eFORGE, LOLA	Tests whether your CpGs cluster in cell-type-specific regulatory elements — a direct check on residual confounding.
Chromatin-state / enhancer context	Roadmap/ENCODE overlaps via LOLA or annotatr	Places hits in enhancers, promoters, and TF-binding regions.
Causal direction (does methylation cause the phenotype?)	Mendelian randomization with mQTLs (e.g. GoDMC)	The rigorous follow-up when a hit is real; well beyond a single EWAS.
“Is my phenotype associated with epigenetic age?”	Epigenetic clocks: methylclock, Horvath (Horvath 2013) / Hannum (Hannum et al. 2013) predictors	A composite outcome computed from a fixed CpG set, not an enrichment of your hits. Fit it as a separate phenotype-level analysis; don’t read it off your EWAS table.

Interpreting the result: a short checklist

Before writing “our EWAS implicates pathway X,” walk through:

1. Did anything survive FDR at the CpG level? If not (as here), say so plainly and treat enrichment as hypothesis-generating, not confirmatory.
2. Is the enrichment robust? Re-run with the full-array background vs your filtered set; a term that only appears with a permissive background is fragile.
3. Does trait enrichment corroborate the biology? Overlap with previously reported CpGs is stronger evidence than a lone GO term.
4. Could composition explain it? Check with eFORGE or by confirming your cell-count adjustment (chapter 04) held.
5. Is the direction consistent? Hyper- vs hypo-methylation across a region or pathway should tell a coherent story; scattered directions suggest noise.
6. How big are the feature sets, and do they overlap? A CpG-centric test on very large, correlated annotation sets will produce many significant features from one underlying signal. Count distinct signals, not rows.
7. Are the hits you are annotating stratum-specific? If they came from a stratified meta-analysis, look up the per-stratum statistics first (notebook 07) — annotating a CpG that only one arm supports is annotating one stratum’s noise.

ImportantTwo caveats this walkthrough inherits, and you should carry forward

Everything in this chapter is downstream of choices made in notebooks 05–07, and two of them limit how far the annotation can be pushed.

The surrogate variables were estimated on the whole cohort. When those same SVs are used inside a single-sex stratum they still carry sex-defined structure, because they were built to describe variation across both sexes. A cleanly stratified analysis would re-estimate SVs within each stratum — which the Snakemake pipeline is not currently set up to do (notebook 07). Treat stratum-specific hits as provisional for this reason alone.

Per-stratum power is very low. The stratified arms spend the same 24 design columns on 45 and 42 samples, leaving roughly 22 and 19 residual degrees of freedom. That is enough to fit the model and not enough to estimate small effects precisely.

Summary

• Annotation makes probe IDs biologically interpretable. Take coordinates, gene model, and CpG-island relation all from one genome build annotation set; don’t mix annotations built from hg19 with those from hg38. The BIOS eQTM catalog provides some annotation for likely effector genes, which can differ from the nearest gene.
• Region analysis recovers signal that single-CpG testing misses by leveraging the spatial correlation of methylation. On our subset an 11-CpG MHC cluster reaches Šidák p = 2.6 × 10⁻¹⁶ while none of its CpGs comes close to the single-CpG EWAS significance threshold.
• Match the DMR tool to your upstream model. comb-p runs on summary statistics and therefore works after stratified meta-analysis. DMRcate needs line-level data and one design formula, so it cannot consume meta-analyzed results.
• Run enrichment multiple ways. missMethyl::gometh for reproducible, probe-bias-corrected gene-set (GO/KEGG) tests; knowYourCG::testEnrichment for CpG-centric tests against regulatory-element databases (CGI, TFBS, chromatin state); and the EWAS Toolkit for the trait catalog you can’t reproduce locally. Read the FDR column, not the nominal p, and always feed the tool your real tested CpG-set. On this subset all three gene-set tools (gometh, methylglm, methylRRA) return nothing at FDR < 0.05, while KYCG returns 41 of 863 regulatory features — a difference driven as much by set size as by biology, and one whose honest reading is a single broad enrichment in promoter-proximal, TF-bound chromatin rather than 41 discoveries.
• Real outputs from this notebook — PTSD_ewas_annotated_zhou.csv.gz (full Zhou-annotated table), the KYCG results table, and the comb-p DMR files — are the biological deliverables of the whole tutorial. The single-CpG hits don’t survive correction on this teaching subset, which is the expected result at n = 87.

NoteOther arrays

Zhou publishes the same hg38 annotation set for every array — swap EPIC for HM450 or EPICv2 in the manifest/gene/SNP filenames and in the platform/array.type arguments to testEnrichment/gometh. The comb-p and DMRcate calls are unchanged apart from the arraytype argument.

References

Bonder, Marc Jan, René Luijk, Daria V Zhernakova, et al. 2017. “Disease Variants Alter Transcription Factor Levels and Methylation of Their Binding Sites.” Nature Genetics 49 (1): 131–38. https://doi.org/10.1038/ng.3721.

Breeze, Charles E, Dirk S Paul, Jenny van Dongen, et al. 2016. “eFORGE: A Tool for Identifying Cell Type-Specific Signal in Epigenomic Data.” Cell Reports 17 (8): 2137–50. https://doi.org/10.1016/j.celrep.2016.10.059.

Goldberg, David C., Lucas Kaminsky, Timothy J. Peters, and Wanding Zhou. 2025. “KnowYourCG: Supervised Learning of CpG Methylation States with Automated Feature-Set Enrichment.” Science Advances 11. https://zhou-lab.github.io/kycg/.

Hannum, Gregory, Justin Guinney, Ling Zhao, et al. 2013. “Genome-Wide Methylation Profiles Reveal Quantitative Views of Human Aging Rates.” Molecular Cell 49 (2): 359–67. https://doi.org/10.1016/j.molcel.2012.10.016.

Horvath, Steve. 2013. “DNA Methylation Age of Human Tissues and Cell Types.” Genome Biology 14 (10): R115. https://doi.org/10.1186/gb-2013-14-10-r115.

Li, Mengwei, Dong Zou, Zhenglin Li, et al. 2019. “EWAS Atlas: A Curated Knowledgebase of Epigenome-Wide Association Studies.” Nucleic Acids Research 47 (D1): D983–88. https://doi.org/10.1093/nar/gky1027.

Pedersen, Brent S, David A Schwartz, Ivana V Yang, and Katerina J Kechris. 2012. “Comb-p: Software for Combining, Analyzing, Grouping and Correcting Spatially Correlated p-Values.” Bioinformatics 28 (22): 2986–88. https://doi.org/10.1093/bioinformatics/bts545.

Peters, Timothy J, Michael J Buckley, Aaron L Statham, et al. 2015. “De Novo Identification of Differentially Methylated Regions in the Human Genome.” Epigenetics & Chromatin 8: 6. https://doi.org/10.1186/1756-8935-8-6.

Phipson, Belinda, Jovana Maksimovic, and Alicia Oshlack. 2016. “missMethyl: An r Package for Analyzing Data from Illumina’s HumanMethylation450 Platform.” Bioinformatics 32 (2): 286–88. https://doi.org/10.1093/bioinformatics/btv560.

Ren, Xu, and Pei Fen Kuan. 2019. “methylGSA: A Bioconductor Package and Shiny App for DNA Methylation Data Length Bias Adjustment in Gene Set Testing.” Bioinformatics 35 (11): 1958–59. https://doi.org/10.1093/bioinformatics/bty892.

Xiong, Zhuang, Fengan Yang, Mengwei Li, et al. 2022. “EWAS Open Platform: Integrated Data, Knowledge and Toolkit for Epigenome-Wide Association Study.” Nucleic Acids Research 50 (D1): D1004–9. https://doi.org/10.1093/nar/gkab972.

Xu, Zongli, Liang Niu, and Jack A. Taylor. 2021. “The ENmix DNA Methylation Analysis Pipeline for Illumina BeadChip and Comparisons to 7 Other Preprocessing Pipelines.” Clinical Epigenetics 13 (1). https://doi.org/10.1186/s13148-021-01207-1.

Xu, Zongli, Changchun Xie, and Jack A. Taylor. 2020. “ipDMR: Identification of Differentially Methylated Regions with Interval P-Values.” Bioinformatics 37 (5): 711–13. https://doi.org/10.1093/bioinformatics/btaa732.

Zhou, Wanding, Peter W Laird, and Hui Shen. 2017. “Comprehensive Characterization, Annotation and Innovative Use of Infinium DNA Methylation BeadChip Probes.” Nucleic Acids Research 45 (4): e22. https://doi.org/10.1093/nar/gkw967.