EWAS

We now have everything an epigenome-wide association study needs: a normalized, filtered methylation matrix (chapter 02–chapter 03), cell-composition estimates (chapter 04), and surrogate variables that absorb technical structure (chapter 05). This notebook fits the association model, checks that it is well-calibrated, and provides guidance for interpreting results. It also lays out two modeling conventions you will meet in the literature — methylation-as-outcome vs. phenotype-as-outcome — and a residualization workflow that is common in practice.

The worked example tests PTSD case/control status against methylation at each CpG, adjusting for sex, age, cell composition, and surrogate variables.

What this chapter needs

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

Terminal

./get_data.sh E_model_inputs

Code

library(limma)   # lmFit(), eBayes(), topTable()
library(bacon)   # bacon(), and the four diagnostic plots in section 4

# The matrix the model is fit to: the stratified-ComBat M-values chapter 05
# wrote, 756,251 probes x 87 samples. It is ~510 MB on disk and about the same
# again in memory, so this read takes several seconds.
Mcb <- readRDS("data/05_mvals_combat.rds")

# Its small companion, which is committed to the repository: the 87-sample
# phenotype table, the six cell proportions from chapter 04, the smoking proxy,
# and the surrogate variables with the k that chapter 05's sweep settled on.
sv_pieces <- readRDS("data/05_sva.rds")
mdk   <- sv_pieces$mdk            # phenotype table, rows in Mcb's column order
propk <- sv_pieces$propk          # CD8T, CD4T, NK, Bcell, Mono, Neu
smoke <- sv_pieces$smoke          # continuous smoking proxy score
K     <- sv_pieces$k_selected     # the k chapter 05 settled on
SV    <- sv_pieces$SV[, seq_len(K), drop = FALSE]

# lmFit() pairs matrix columns with design rows by position and does not check
# that they describe the same samples, so confirm the counts agree here.
stopifnot(ncol(Mcb) == nrow(mdk))

# Control is the reference level, which is what makes the tested coefficient
# `ptsdCase` and makes its sign read as case minus control rather than the
# reverse. Left to alphabetical ordering, "Case" would become the reference.
ptsd    <- relevel(factor(mdk$ptsd), ref = "Control")
is_case <- ptsd == "Case"

# `ew` is the handoff object this chapter saves as data/06_ewas.rds and that
# chapter 08 reads back. The sample counts are known before any model is fit;
# the design size, the results table and lambda are added below as they are
# computed, and the whole list is written out in one place at the end.
ew <- list(n = ncol(Mcb), ncase = sum(is_case), nctrl = sum(!is_case),
           k_sv = K, n_tested = nrow(Mcb),
           panel_excluded = sv_pieces$excluded_n)

cat("EWAS matrix:", nrow(Mcb), "probes x", ncol(Mcb), "samples | k =", K, "\n")

EWAS matrix: 756251 probes x 87 samples | k = 6

1. Two ways to write the model — and why the outcome matters

There are two ways to set up the per-CpG regression, and they answer subtly different questions.

(a) Methylation as the outcome:

\[\text{methylation}_{\text{CpG}} \sim \text{PTSD} + \text{sex} + \text{age} + \text{cell composition} (+ \text{SVs})\]

(b) Phenotype as the outcome:

\[\text{PTSD} \sim \text{methylation}_{\text{CpG}} + \text{sex} + \text{age} + \text{cell composition} (+ \text{SVs})\]

Both are used, but methylation-as-outcome is the field standard for EWAS for several reasons:

Methylation is the measured, continuous response. Treating the β/M-value as the dependent variable matches the data-generating process: you measured methylation and ask what it varies with.
It generalizes to any exposure type. With methylation as outcome, the predictor of interest can be binary (PTSD), continuous (BMI, age, pack-years), or ordinal — the model form doesn’t change. Phenotype-as-outcome forces a different link function for each exposure type (logistic for binary, linear for continuous), complicating a genome-wide pipeline.
Effect sizes are biologically legible. A coefficient in the methylation-outcome model is “change in methylation per unit exposure” — a Δβ you can read as a percent methylation difference.
It is what meta-analysis tools expect. METAL and comparable tools combine per-CpG effect sizes and standard errors on a common scale; the methylation-outcome β/SE is that common currency across cohorts.

When would you flip it to phenotype-as-outcome? Mainly when methylation is conceived as a predictor/biomarker of the phenotype — e.g. building a methylation risk score, or asking “does methylation at this locus predict disease?” rather than “does disease associate with methylation here?”.

Test on M-values, report on β-values

Regardless of the outcome convention, the statistical test is run on M-values (logit-transformed β), which are approximately homoscedastic and better satisfy the model’s variance assumptions (Du et al. 2010). Effect sizes are then reported as Δβ because a β difference is biologically interpretable — a Δβ of 0.05 is a “5% methylation change.”

2. The covariate set

Our design includes covariates for everything we established as a driver of methylation:

sex and age — strong, ubiquitous methylation covariates (chapter 05 confirmed both drive leading PCs).
cell composition (6 estimated proportions) the single largest source of variance in whole blood (chapter 04).
a smoking proxy — this cohort has no recorded smoking variable, so chapter 05 builds one from a published panel of replicated smoking-associated CpGs. The panel probes themselves are excluded from testing.
array position — the chip row (R01C01…R08C01), a recorded technical factor that is estimable here and so enters as a fixed effect (7 dummy columns).
surrogate variables — the 6 SVs retained in chapter 05, which absorb structure that none of the measured covariates name.

This is the step where the sample count changes. Chapters 01–05 worked with all 96 QC-passing arrays, because none of those steps needs the exposure. The model above does, so the nine samples with no recorded PTSD status drop out and the EWAS fits n = 87 (32 cases, 55 controls). No sample was excluded for quality — see phenotype completeness.

Smoking

Smoking is well established to cause distinct, systemic methylation changes, so it is important to include as a covariate in modeling if you have that information available. If you don’t, there are methods to predict smoking status from a subset of methylation sites (these sites would then be excluded from the EWAS). Alternatively, you can exclude CpGs identified as strongly associated with smoking from previous literature prior to or post EWAS (prior is recommended if using this option). Lastly, surrogate variables can capture latent variation due to smoking if you don’t have a smoking variable included in the null model, but use this option with caution as the surrogate variables may attenuate more biological signal than intended. As a last resort, you can run the EWAS without accounting for smoking, just keep this in mind while interpreting results.

Code

# Array position is an 8-level factor on an EPIC chip, so it costs 7 dummy
# columns. Rebuilding the factor from as.character() drops any level the sample
# sheet defines but no modeled sample occupies; an empty level would add an
# all-zero dummy column and leave the design rank-deficient.
pos <- factor(as.character(mdk$array_pos))

# Every model term goes into one data frame so that model.matrix() can name the
# columns it builds: `propk` expands to the six cell proportions by their own
# names, and `SV` to SV1...SV6. check.names = FALSE keeps those names intact.
dd <- data.frame(ptsd = ptsd, sex = factor(mdk$sex), age = mdk$age,
                 smoke = smoke, pos = pos, propk, SV = SV, check.names = FALSE)
design <- model.matrix(~ ptsd + sex + age + smoke + pos +
                         CD8T + CD4T + NK + Bcell + Mono + Neu + SV, data = dd)
stopifnot("ptsdCase" %in% colnames(design))
ew$design_ncol <- ncol(design)
ew$resid_df    <- ncol(Mcb) - ncol(design)

# 756,251 least-squares fits in a single call -- which is why limma is the tool
# here rather than a loop over lm(). This is the expensive step in the chapter:
# roughly a minute, and several GB while the per-probe coefficients and residual
# variances are held in memory.
fit <- lmFit(Mcb, design)        # per-probe weighted least squares
fit <- eBayes(fit)               # empirical-Bayes moderation of the variance

# sort.by = "P" fixes the row order of the saved table to increasing p-value;
# topTable()'s default sorts by the log-odds column B instead.
#
# limma:: is required, not decoration. bacon exports its own topTable method,
# and because bacon is attached after limma above, a bare topTable() call
# resolves to bacon's -- whose signature is (object, number, adjust.method,
# sort.by) and which fails with "unused argument (coef = ...)". Qualifying the
# call is the fix that does not depend on library() order.
tt <- limma::topTable(fit, coef = "ptsdCase", number = Inf, sort.by = "P")
tt$probe <- rownames(tt)

# Delta-beta: the case-minus-control difference on the beta scale, taken from
# the same ComBat-adjusted M-values the model saw, so the effect reported and
# the effect tested describe one matrix. The back-transform is a second copy of
# 756,251 x 87 doubles, so it is dropped again as soon as the means are taken.
bcb <- 2^Mcb / (2^Mcb + 1)
db  <- rowMeans(bcb[, is_case, drop = FALSE],  na.rm = TRUE) -
       rowMeans(bcb[, !is_case, drop = FALSE], na.rm = TRUE)
tt$delta_beta <- db[tt$probe]
rm(bcb); invisible(gc())

# Genomic inflation factor from the p-values: the conventional definition, and
# the same one chapter 05's k-selection sweep used, so the two are comparable.
lam    <- function(p) median(qchisq(1 - p, 1), na.rm = TRUE) / qchisq(0.5, 1)
lambda <- lam(tt$P.Value)

# Fix the column order of the table that gets saved and deposited.
tt <- tt[, c("probe", "logFC", "delta_beta", "AveExpr", "t", "P.Value",
             "adj.P.Val", "B")]

limma’s eBayes step is what makes this appropriate for an EWAS on a modest sample size: it moderates each probe’s variance toward a pooled estimate, stabilizing the t-statistics for the ~700k probes tested (Ritchie et al. 2015).

Code

cat(sprintf("EWAS sample size: n = %d  (%d cases, %d controls)\n", ew$n, ew$ncase, ew$nctrl))

EWAS sample size: n = 87  (32 cases, 55 controls)

Code

cat(sprintf("Design columns (covariates + intercept): %d\n", ew$design_ncol))

Design columns (covariates + intercept): 24

Code

cat(sprintf("Residual degrees of freedom per probe:   %d\n", ew$resid_df))

Residual degrees of freedom per probe:   63

Code

cat(sprintf("Surrogate variables retained:            %d\n", ew$k_sv))

Surrogate variables retained:            6

Code

cat(sprintf("CpGs tested: %s  (%d smoking-panel probes excluded)\n",
            format(ew$n_tested, big.mark = ","), ew$panel_excluded))

CpGs tested: 756,251  (20 smoking-panel probes excluded)

3. Is the model well-calibrated? The QQ plot and genomic inflation

Before looking at which CpGs are top-ranked, check whether the p-values behave. The genomic inflation factor λ compares the median observed test statistic to its null expectation. λ ≈ 1 means the model is well-calibrated; λ ≫ 1 signals residual confounding (unmodeled batch, cell composition, or population structure) inflating the statistics; λ ≪ 1 can signal over-correction, but may also indicate that there truly is very little differential methylation or the sample size is underpowered to detect true signals. Due to the biological differences between genetic variants and methylation and the technical design differences in how they are measured in arrays, EWAS results tend to have inflated λ compared to GWAS results. Additional bias and inflation correction is performed post-EWAS to account for this. The lambda post-bias-and-inflation correction is what should be used to judge the model fit.

4. Correcting residual bias and inflation with BACON

BACON (Iterson et al. 2017) estimates the empirical null distribution of the test statistics with a Bayesian three-component mixture (one null component plus two alternative components for positive and negative effects), then rescales every statistic to that estimated null. This corrects both bias (a non-zero center of the null — a systematic shift) and inflation (a null wider than N(0,1) — the same phenomenon λ measures), and it does so without assuming that all deviation from the diagonal is confounding.

Always run BACON on your EWAS results:

A λ near 1 does not mean there is no residual bias — bias (a shift of the null’s center) barely moves λ, which is driven by the spread of the statistics. BACON estimates and removes both.
Applying the empirical-null correction and inspecting its diagnostics is what lets you distinguish “well-calibrated” from “looks calibrated by one number.”

BACON works from either t-statistics or effect sizes and their standard errors. We chose effect sizes and standard errors as we can derive these from the output of limma and these estimates are needed down stream if a meta-analysis is performed.

Code

# The Gibbs sampler is stochastic, so the seed is what makes the inflation,
# the bias and every adjusted p-value below reproducible run to run.
set.seed(42)

es <- tt$logFC              # effect size per CpG
se <- tt$logFC / tt$t       # standard error = beta / t-statistic

# 5,000 iterations with 2,000 burn-in (bacon's defaults) over all 756,251
# statistics: about a minute. Passing NULL for the test statistics tells bacon
# to work from the effect size / standard error pair instead.
bc <- bacon(NULL, effectsizes = es, standarderrors = se)

# es() and se() below are bacon's accessor functions, not the two numeric
# vectors of the same name defined above: R resolves a name used in function
# position to a function and skips the vector.
inflation(bc)   # estimated inflation (sigma.0); 1 = none

 sigma.0 
1.004763

Code

bias(bc)        # estimated bias (mu); 0 = none

       mu.0 
0.003272855

Code

tt$bacon.p  <- pval(bc)     # empirical-null-adjusted p-values
tt$bacon.es <- es(bc)       # adjusted effect sizes
tt$bacon.se <- se(bc)       # adjusted standard errors

# BH across the adjusted p-values. The results section thresholds on this
# column, so it has to be computed from bacon.p rather than reused from limma.
tt$bacon.adj.P <- p.adjust(tt$bacon.p, method = "BH")

Code

# The four calibration numbers, collected into the object saved at the end of
# the chapter as data/06_bacon_summary.rds. lambda_raw is limma's, lambda_bacon
# is the same statistic recomputed on the empirical-null-adjusted p-values.
bs <- list(inflation    = inflation(bc),
           bias         = bias(bc),
           lambda_raw   = lambda,
           lambda_bacon = lam(tt$bacon.p),
           n            = ew$n)
cat(sprintf("Estimated inflation (sigma.0): %.3f\n", bs$inflation))

Estimated inflation (sigma.0): 1.005

Code

cat(sprintf("Estimated bias (mu):           %.3f\n", bs$bias))

Estimated bias (mu):           0.003

Code

cat(sprintf("lambda before BACON:           %.3f\n", bs$lambda_raw))

lambda before BACON:           1.022

Code

cat(sprintf("lambda after  BACON:           %.3f\n", bs$lambda_bacon))

lambda after  BACON:           1.023

Here the model was already well-calibrated (λ ≈ 1.02), so BACON changes little — the estimated inflation is ≈ 1 and the adjusted λ stays close to 1. That is the reassuring outcome: the correction is available and applied, and it confirms rather than rescues the analysis. On a poorly-specified model you would instead see inflation ≫ 1 pulled back toward 1, and the diagnostics below would show it.

Always inspect the four BACON diagnostic plots

The single number can hide a mis-specified or non-converged fit. BACON’s Gibbs sampler produces four diagnostic plots, and you should look at all four every time. (The interpretation below follows the guidance in the pipeline repository README, where these are documented in full.)

Code

# The four plots below are drawn from `bc` and written to disk here, so the
# figures that follow display this run rather than a stored copy of an older
# one. They are base and lattice graphics rather than ggplots, so each goes to
# its own jpeg() device instead of through ggsave().
dir.create("data/06_bacon", showWarnings = FALSE, recursive = TRUE)

# fit() is bacon's mixture-fit plot. The limma object named `fit` above does not
# get in the way, for the same reason es() and se() still worked: a name in
# function position resolves to a function.
jpeg("data/06_bacon/fit.jpg", width = 1600, height = 1200, res = 200)
print(fit(bc, n = 100))                  # n = 100 histogram bins

NULL

Code

dev.off()

png 
  2

Code

jpeg("data/06_bacon/qq.jpg", width = 1600, height = 1200, res = 200)
print(plot(bc, type = "qq"))             # a ggplot, hence the print()
dev.off()

png 
  2

Code

# traces() and posteriors() draw straight onto the device, so no print() here.
# burnin = FALSE shows the chain after burn-in only, which is the part you
# judge convergence on.
jpeg("data/06_bacon/traces.jpg", width = 1600, height = 1400, res = 180)
traces(bc, burnin = FALSE)
dev.off()

png 
  2

Code

jpeg("data/06_bacon/posteriors.jpg", width = 1600, height = 1200, res = 200)
posteriors(bc)
dev.off()

png 
  2

**Fit plot.** Histogram of the observed z-scores with the fitted three-component mixture overlaid: the black curve is the overall fit, red is the estimated null component, and blue/green are the two alternative components. The black curve tracks the histogram closely — the model is well-specified. A large gap between the fitted (black) and observed (histogram) densities would signal mis-specification or non-convergence.

**QQ plots, before (uncorrected) and after (corrected) BACON.** Points should hug the diagonal, breaking upward into a tail spike only for genuine associations. Systematic departure from the line *before* the spike is bias/inflation; after correction the bulk of points should sit closer to the diagonal. Here the two panels are nearly identical because the input was already calibrated.

**Traces plot.** The Gibbs sampler’s value for each parameter (p, mu, sigma for each of the three components) across iterations. A well-mixed, converged chain looks like a fuzzy horizontal band with a narrow range in estimates — a ‘hairy caterpillar’ — with no drift or trend. A sharp early sweep followed by settling is normal (the prior was far from the estimate); a chain that never settles means non-convergence.

**Posteriors plot.** Scatter of the sampled posterior draws for the null inflation (sigma.0) and the null proportion (p.0), with 75/90/95% probability ellipses. A dense, roughly elliptical cloud sitting inside the contours indicates the sampler converged on a stable estimate; a sparse or lopsided cloud spilling outside the ellipses indicates it did not.

In this run the fit tracks the histogram, the QQ panels are near-diagonal, the traces are hairy caterpillars, and the posteriors form a dense cloud inside the ellipses — all four say the sampler converged and the empirical null is trustworthy.

Make BACON a default, and read the diagnostics

Report BACON-adjusted effect sizes and p-values as your primary EWAS results. It is not usually necessary to report the diagnostic plots, but I strongly recommend reporting the adjusted lambda.

5. Results

Code

bonf   <- 0.05 / nrow(tt)
n_bonf <- sum(tt$bacon.p < bonf)
n_fdr  <- sum(tt$bacon.adj.P < 0.05)
cat(sprintf("Probes tested: %s\n", format(nrow(tt), big.mark = ",")))

Probes tested: 756,251

Code

cat(sprintf("Bonferroni-significant (p < %.2g): %d\n", bonf, n_bonf))

Bonferroni-significant (p < 6.6e-08): 1

Code

cat(sprintf("FDR < 0.05: %d\n", n_fdr))

FDR < 0.05: 2

Code

knitr::kable(head(tt[order(tt$bacon.p),
                     c("probe","logFC","delta_beta","P.Value","bacon.p","bacon.adj.P")], 6),
             digits = c(0, 3, 4, 8, 8, 4), row.names = FALSE,
             col.names = c("probe", "logFC", "delta_beta", "P (limma)",
                           "P (BACON)", "FDR (BACON)"),
             caption = "Top 6 CpGs, ranked by BACON-adjusted p-value. logFC is the M-value effect; delta_beta is the case−control β difference (interpretable as percent methylation).")

Top 6 CpGs, ranked by BACON-adjusted p-value. logFC is the M-value effect; delta_beta is the case−control β difference (interpretable as percent methylation).
probe	logFC	delta_beta	P (limma)	P (BACON)	FDR (BACON)
cg15434749	-0.394	-0.0254	5.00e-07	2.967175e-08	0.0224
cg20600436	0.172	0.0138	8.60e-07	6.731565e-08	0.0255
cg25717994	0.308	0.0057	2.38e-06	2.825682e-07	0.0627
cg02045051	-0.176	-0.0076	2.73e-06	3.313880e-07	0.0627
cg05825555	0.213	0.0237	8.71e-06	1.653689e-06	0.2360
cg24228058	-0.231	-0.0126	1.03e-05	2.002865e-06	0.2360

Only 1 CpG passes Bonferroni and 2 pass FDR < 0.05 — a borderline result that a study of this size cannot settle either way. With n = 87, a single marginal hit is not unexpected.

The next two figures are the same Manhattan plot drawn twice, once on limma’s uncorrected p-values and once on the BACON-adjusted ones. Plotting both is worth the space: it shows exactly what the correction did to your results, and it makes clear that the one CpG which clears Bonferroni does so because of the rescaling, not independently of it. Report the adjusted version.

Manhattan plot of the EWAS, drawn on limma’s uncorrected p-values. Nothing crosses the Bonferroni line (red) on this scale — the single CpG that clears Bonferroni does so only after BACON rescaling, which shifts it just past the threshold. A near-empty Manhattan is the expected result for a subtle psychiatric phenotype in a small teaching subset (n=87); the full-scale analysis (n≈800) and meta-analysis across cohorts is what powers genome-wide discovery for PTSD (Katrinli et al. 2020).

The same Manhattan plot drawn on the BACON-adjusted p-values — the ones the chapter recommends you report. The axes, the coordinates, and the Bonferroni line are identical to the previous figure; only the y-values change. Rescaling by the empirical null lifts `cg15434749` (filled point) just across Bonferroni and leaves `cg20600436` (open ring) passing FDR < 0.05 alone. Comparing the two panels is the point: BACON does not manufacture signal, it re-centres and re-scales the whole test-statistic distribution, and a CpG that was already at the top of the uncorrected plot is the one that moves across.

Volcano plot. The x-axis is the EWAS effect size itself — the `ptsdCase` coefficient from the limma fit, on the M-value scale, adjusted for every covariate in the design. Colour marks the sign of that coefficient: teal for CpGs less methylated in cases, plum for those more methylated. Note that this is the adjusted effect, so its sign can differ from a raw case-minus-control β difference for probes where the covariates carry part of the signal. The three most significant probes are ringed and labelled. Even they move methylation only slightly, illustrating that EWAS effect sizes for complex phenotypes are typically small and require large samples to detect reliably.

6. An alternative in common use: the M-value residual workflow

The design above puts all covariates — biological and technical — into one large regression. A common alternative, and the convention many groups (including this tutorial’s pipeline author) prefer, is a two-stage residualization:

Regress the technical variation out first. Fit each CpG’s M-value on only the technical/nuisance terms (surrogate variables, or recorded batch such as chip and position) and keep the residuals — the methylation signal with technical structure removed.
Run the association on the residuals with only the biological model: \[\text{resid}(\text{M})_{\text{CpG}} \sim \text{PTSD} + \text{sex} + \text{age} + \text{cell composition}\]

# Stage 1: residualize M-values on the technical terms (SVs here; or chip + position)
tech <- model.matrix(~ SV)                       # nuisance-only design
M_resid <- residuals(lmFit(M, tech), M)          # limma residuals, probes x samples

# Stage 2: association model on the residuals, biological covariates only
bio <- model.matrix(~ ptsd + sex + age + Neu + NK + CD4T + CD8T + Bcell + Mono, data = md)
fit <- eBayes(lmFit(M_resid, bio))

Why do this? The usual motivation is degrees of freedom. Our single-stage model spends 24 design columns per probe; at n = 87 that leaves 63 residual df. Moving the technical terms into a first stage — where they are estimated once — leaves the association model lean (PTSD + sex + age + 6 cell props) and returns df to the test that matters.

It is also common to regress out the blood cell proportions in addition to the technical variation prior to running the main model. There is no wrong choice in which or how many covariates you adjust for prior to the EWAS, just make sure to report what exactly was done in your methods.

The caveat: residualization understates uncertainty

Two-stage residualization treats the stage-1 fitted values as known when it runs stage 2, so it doesn’t propagate the uncertainty of estimating the technical terms. In practice, with many samples relative to nuisance terms, the difference from the single-stage model is small and the df savings are worth it — which is why it is widely used. But be aware it is an approximation to the joint model, not identical to it. Report which approach you used, and ideally confirm the top hits are stable under both.

7. When and why to stratify, then meta-analyze

The model so far adjusts for sex as a covariate: it lets men and women have different mean methylation, but forces the PTSD effect to be a single number shared across both. That assumption is often wrong for methylation.

Why sex is a special case in methylation

Sex is one of the largest sources of variation in the methylome, and it is not just an X/Y-chromosome effect — autosomal CpGs show widespread sex differences in mean level, in variance, and in how they respond to exposures and phenotypes. When the effect of your exposure/trait/disease genuinely differs between men and women (a sex interaction), a sex-adjusted model estimates a blend of the two and can miss a real signal that is present in only one sex, or that points in opposite directions. Because sex-specific methylation effects are common, I strongly recommend performing sex-stratified analyses.

More generally, stratify when you have reason to believe the effect itself (not just the baseline) differs across a grouping variable: genetic ancestry, tissue or cell subtype, exposure dose, or age categories can all justify it. The cost is power — each stratum has fewer samples — so stratification trades a cleaner, assumption-light estimate within each group against reduced within-stratum sample size. However, the loss of power from sample size can be recovered with meta-analysis to re-join the strata results.

Stratify → analyze → meta-analyze

The workflow is three steps:

Split the data by the stratifying variable (e.g. sex) and drop that variable from the covariate set within each stratum — it is now constant.
Run the full EWAS separately in each stratum: the same model, the same BACON adjustment, the same diagnostics, per stratum.
Meta-analyze the per-stratum results with an inverse-variance-weighted fixed-effect meta-analysis, which combines the two effect estimates weighting each by its precision (1/SE²). This recovers power from both strata while still allowing the underlying effects to differ — and it gives you a heterogeneity test (Cochran’s Q, I²) that formally asks whether the sexes differ at each CpG.

Inverse-variance meta-analysis of a single CpG is short enough to show by hand:

# per-CpG effect + SE from each stratum's BACON-adjusted results
# (bF, seF) from females; (bM, seM) from males
wF <- 1 / seF^2                       # inverse-variance weights
wM <- 1 / seM^2
beta_meta <- (wF*bF + wM*bM) / (wF + wM)      # combined effect
se_meta   <- sqrt(1 / (wF + wM))              # combined standard error
z_meta     <- beta_meta / se_meta
p_meta     <- 2 * pnorm(-abs(z_meta))

# heterogeneity: does the effect differ between the sexes?
Q  <- wF*(bF - beta_meta)^2 + wM*(bM - beta_meta)^2   # Cochran's Q, 1 df
p_het <- pchisq(Q, df = 1, lower.tail = FALSE)

In practice you run this across all CpGs at once. The standard tool for it is METAL (Willer et al. 2010) — the same command-line meta-analyzer used for GWAS. It reads each stratum’s results file, matches CpGs by ID, does the inverse-variance weighting and the heterogeneity test genome-wide, and writes a combined table. A concordant hit (same direction in both sexes, low heterogeneity) is far more credible than a hit seen in only one stratum. Chapter 07 runs the stratified EWAS and METAL as a managed workflow.

You do not need METAL for this

METAL’s SCHEME STDERR is inverse-variance-weighted fixed-effects meta-analysis. The single-CpG code above vectorizes directly, so the whole genome-wide meta-analysis is a dozen lines of R over the two strata’s BACON-adjusted result tables:

library(data.table)

# each stratum's BACON-adjusted results: one row per CpG
# Written by the stratified pipeline run in chapter 07; paths are relative to
# the tutorial folder. If you have not run the pipeline, `./get_data.sh
# G_pipeline_run` fetches the published run these two files come from.
mF <- fread("../ewas_pipeline/run_grady/F/F_PTSD_ewas_bacon_results.csv.gz")[
        , .(probe = cpgid, bF = bacon.es, seF = bacon.se)]
mM <- fread("../ewas_pipeline/run_grady/M/M_PTSD_ewas_bacon_results.csv.gz")[
        , .(probe = cpgid, bM = bacon.es, seM = bacon.se)]
mt <- merge(mF, mM, by = "probe")          # CpGs tested in both strata

mt[, wF := 1 / seF^2][, wM := 1 / seM^2]                      # precision weights
mt[, beta_meta := (wF * bF + wM * bM) / (wF + wM)]            # combined effect
mt[, se_meta   := sqrt(1 / (wF + wM))]                        # combined SE
mt[, z_meta    := beta_meta / se_meta]
mt[, p_meta    := 2 * pnorm(-abs(z_meta))]

# heterogeneity, per CpG
mt[, Q     := wF * (bF - beta_meta)^2 + wM * (bM - beta_meta)^2]   # Cochran's Q, 1 df
mt[, p_het := pchisq(Q, df = 1, lower.tail = FALSE)]
mt[, I2    := pmax(0, (Q - 1) / Q) * 100]

# direction string, the way METAL reports it
mt[, direction := paste0(ifelse(bF > 0, "+", "-"), ifelse(bM > 0, "+", "-"))]

setorder(mt, p_meta)

This is not an approximation. Run against the same two strata this tutorial meta-analyzes with METAL, across all 756,251 CpGs tested in both, it returns the identical answer: the same top hit (cg22671410, p = 2.027 × 10⁻¹³, direction +-, I² = 0), the same 9 CpGs below the Bonferroni threshold of 6.61 × 10⁻⁸, the same 50 at FDR < 0.05, and a rank correlation of 1 between the two sets of p-values. METAL writes its effect and standard error columns rounded to four decimal places, and once you round the R result the same way, the largest disagreement across all 756,251 CpGs is exactly zero.

So when is METAL worth it?

Not for the arithmetic — for the bookkeeping. METAL handles allele/strand alignment (irrelevant for methylation, central for GWAS), reads many cohort files with inconsistent column names, tracks sample sizes per marker, and applies genomic-control correction on request. For a two-stratum methylation meta-analysis where you control both input files, the R above is easier to audit and keeps everything in one language. For a multi-cohort consortium meta-analysis, use METAL.

8. When you don’t have IDATs

Everything in this notebook operates on a methylation matrix and a phenotype table, so the non-IDAT path is identical once you have M-values:

Convert a deposited β-matrix to M-values (M = log2(beta/(1-beta)); see normalization).
Obtain covariates from the series metadata: sex and age from pData, cell composition via a reference-based estimator that runs on β-matrices (e.g. EpiDISH/RPC) when you can’t run estimateCellCounts2 from IDATs (chapter 04 covers this), and surrogate variables from sva run on the matrix (chapter 05).
Fit the identical limma model. The association step never touches raw intensities.

Saving what we just computed

The fit and the BACON correction are the whole cost of this chapter. Three files carry them forward: the object chapter 08 reads, the calibration summary section 4 reports from, and a flat copy of the results table.

Code

# 1. The handoff object. `tt` is attached only now, after BACON, so that the
#    saved table already carries bacon.p / bacon.es / bacon.se / bacon.adj.P --
#    chapter 08 reads those columns straight out of this file and never reruns
#    the sampler.
ew$tt     <- tt
ew$lambda <- lambda
saveRDS(ew, "data/06_ewas.rds")

# 2. The four calibration numbers, on their own, for anything that needs them
#    without loading a 756,251-row table.
saveRDS(bs, "data/06_bacon_summary.rds")

# 3. A flat copy of the results table, one row per CpG. Gzipped CSV rather than
#    a tab-delimited file in the repository: 756,251 rows is far too large to
#    commit, so this one is deposited in the Zenodo record instead. The write
#    takes about a minute.
write.csv(tt, gzfile("data/06_ewas_bacon_toptable.csv.gz"), row.names = FALSE)

Picking the chapter up from the saved fit instead

Fitting 756,251 probes needs the 510 MB ComBat matrix in memory with a few GB of working space on top of it. If that is not practical on your machine, fetch the published fit and read it instead of running the two chunks above:

Terminal

./get_data.sh F_ewas_results

RStudio Console

ew <- readRDS("data/06_ewas.rds")
tt <- ew$tt                                   # already carries the BACON columns
bs <- readRDS("data/06_bacon_summary.rds")

Every number and table in this chapter will then run. The four diagnostic plots will not, because they are drawn from the bc object rather than from the saved table.

Next: Running the EWAS pipeline — scale this analysis with a reproducible Snakemake workflow that models methylation as the outcome, runs a sex-stratified EWAS, and meta-analyzes the strata with METAL. Then functional annotation of the results.

References

Du, Pan, Xiao Zhang, Chiang-Ching Huang, et al. 2010. “Comparison of Beta-Value and m-Value Methods for Quantifying Methylation Levels by Microarray Analysis.” BMC Bioinformatics 11: 587. https://doi.org/10.1186/1471-2105-11-587.

Iterson, Maarten van, Erik W van Zwet, Bastiaan T Heijmans, and the BIOS Consortium. 2017. “Controlling Bias and Inflation in Epigenome- and Transcriptome-Wide Association Studies Using the Empirical Null Distribution.” Genome Biology 18 (1): 19. https://doi.org/10.1186/s13059-016-1131-9.

Katrinli, Seyma, Adriana Lori, Varun Kilaru, et al. 2020. “Association of HLA Locus Alleles with Posttraumatic Stress Disorder.” Brain, Behavior, and Immunity 87: 37–45. https://doi.org/10.1016/j.bbi.2020.04.038.

Ritchie, Matthew E, Belinda Phipson, Di Wu, et al. 2015. “Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies.” Nucleic Acids Research 43 (7): e47. https://doi.org/10.1093/nar/gkv007.

Willer, Cristen J, Yun Li, and Gonçalo R Abecasis. 2010. “METAL: Fast and Efficient Meta-Analysis of Genomewide Association Scans.” Bioinformatics 26 (17): 2190–91. https://doi.org/10.1093/bioinformatics/btq340.

--- title: "EWAS" --- ```{r} #| label: setup #| include: false source("_setup.R") ``` We now have everything an epigenome-wide association study needs: a normalized, filtered methylation matrix ([chapter 02](02_normalization.qmd)–[chapter 03](03_probe_filtering.qmd)), cell-composition estimates ([chapter 04](04_cell_composition.qmd)), and surrogate variables that absorb technical structure ([chapter 05](05_batch_effects.qmd)). This notebook fits the association model, checks that it is well-calibrated, and provides guidance for interpreting results. It also lays out **two modeling conventions** you will meet in the literature — methylation-as-outcome vs. phenotype-as-outcome — and a residualization workflow that is common in practice. The worked example tests **PTSD case/control status** against methylation at each CpG, adjusting for sex, age, cell composition, and surrogate variables. ::: {.callout-note} ## What this chapter needs This chapter starts from `data/05_mvals_combat.rds`, written by [chapter 05](05_batch_effects.qmd). If you did not run that chapter, fetch the published checkpoint instead: ```{.bash filename="Terminal"} ./get_data.sh E_model_inputs ``` ::: ```{r} #| label: load library(limma) # lmFit(), eBayes(), topTable() library(bacon) # bacon(), and the four diagnostic plots in section 4 # The matrix the model is fit to: the stratified-ComBat M-values chapter 05 # wrote, 756,251 probes x 87 samples. It is ~510 MB on disk and about the same # again in memory, so this read takes several seconds. Mcb <- readRDS("data/05_mvals_combat.rds") # Its small companion, which is committed to the repository: the 87-sample # phenotype table, the six cell proportions from chapter 04, the smoking proxy, # and the surrogate variables with the k that chapter 05's sweep settled on. sv_pieces <- readRDS("data/05_sva.rds") mdk <- sv_pieces$mdk # phenotype table, rows in Mcb's column order propk <- sv_pieces$propk # CD8T, CD4T, NK, Bcell, Mono, Neu smoke <- sv_pieces$smoke # continuous smoking proxy score K <- sv_pieces$k_selected # the k chapter 05 settled on SV <- sv_pieces$SV[, seq_len(K), drop = FALSE] # lmFit() pairs matrix columns with design rows by position and does not check # that they describe the same samples, so confirm the counts agree here. stopifnot(ncol(Mcb) == nrow(mdk)) # Control is the reference level, which is what makes the tested coefficient # `ptsdCase` and makes its sign read as case minus control rather than the # reverse. Left to alphabetical ordering, "Case" would become the reference. ptsd <- relevel(factor(mdk$ptsd), ref = "Control") is_case <- ptsd == "Case" # `ew` is the handoff object this chapter saves as data/06_ewas.rds and that # chapter 08 reads back. The sample counts are known before any model is fit; # the design size, the results table and lambda are added below as they are # computed, and the whole list is written out in one place at the end. ew <- list(n = ncol(Mcb), ncase = sum(is_case), nctrl = sum(!is_case), k_sv = K, n_tested = nrow(Mcb), panel_excluded = sv_pieces$excluded_n) cat("EWAS matrix:", nrow(Mcb), "probes x", ncol(Mcb), "samples | k =", K, "\n") ``` ## 1. Two ways to write the model — and why the outcome matters {#sec-model} There are two ways to set up the per-CpG regression, and they answer subtly different questions. **(a) Methylation as the outcome**: $$\text{methylation}_{\text{CpG}} \sim \text{PTSD} + \text{sex} + \text{age} + \text{cell composition} (+ \text{SVs})$$ **(b) Phenotype as the outcome**: $$\text{PTSD} \sim \text{methylation}_{\text{CpG}} + \text{sex} + \text{age} + \text{cell composition} (+ \text{SVs})$$ Both are used, but **methylation-as-outcome** is the field standard for EWAS for several reasons: - **Methylation is the measured, continuous response.** Treating the β/M-value as the dependent variable matches the data-generating process: you measured methylation and ask what it varies with. - **It generalizes to any exposure type.** With methylation as outcome, the predictor of interest can be binary (PTSD), continuous (BMI, age, pack-years), or ordinal — the model form doesn't change. Phenotype-as-outcome forces a different link function for each exposure type (logistic for binary, linear for continuous), complicating a genome-wide pipeline. - **Effect sizes are biologically legible.** A coefficient in the methylation-outcome model is "change in methylation per unit exposure" — a Δβ you can read as a percent methylation difference. - **It is what meta-analysis tools expect.** METAL and comparable tools combine per-CpG **effect sizes and standard errors** on a common scale; the methylation-outcome β/SE is that common currency across cohorts. When would you flip it to phenotype-as-outcome? Mainly when methylation is conceived as a **predictor/biomarker** of the phenotype — e.g. building a methylation risk score, or asking "does methylation at this locus predict disease?" rather than "does disease associate with methylation here?". ::: {.callout-note} ## Test on M-values, report on β-values Regardless of the outcome convention, the statistical test is run on **M-values** (logit-transformed β), which are approximately homoscedastic and better satisfy the model's variance assumptions [@du2010comparison]. Effect sizes are then reported as **Δβ** because a β difference is biologically interpretable — a Δβ of 0.05 is a "5% methylation change." ::: ## 2. The covariate set {#sec-covariates} Our design includes covariates for everything we established as a driver of methylation: - **sex** and **age** — strong, ubiquitous methylation covariates ([chapter 05](05_batch_effects.qmd) confirmed both drive leading PCs). - **cell composition** (6 estimated proportions) the single largest source of variance in whole blood ([chapter 04](04_cell_composition.qmd)). - **a smoking proxy** — this cohort has no recorded smoking variable, so [chapter 05](05_batch_effects.qmd) builds one from a published panel of replicated smoking-associated CpGs. The panel probes themselves are excluded from testing. - **array position** — the chip row (`R01C01`…`R08C01`), a recorded technical factor that is estimable here and so enters as a fixed effect (7 dummy columns). - **surrogate variables** — the `r ew$k_sv` SVs retained in [chapter 05](05_batch_effects.qmd), which absorb structure that none of the measured covariates name. This is the step where the sample count changes. Chapters 01–05 worked with all **96** QC-passing arrays, because none of those steps needs the exposure. The model above does, so the nine samples with no recorded PTSD status drop out and the EWAS fits **n = `r ew$n`** (`r ew$ncase` cases, `r ew$nctrl` controls). No sample was excluded for quality — see [phenotype completeness](01_qc.qmd). ::: {.callout-note} ## Smoking Smoking is well established to cause distinct, systemic methylation changes, so it is important to include as a covariate in modeling if you have that information available. If you don't, there are methods to predict smoking status from a subset of methylation sites (these sites would then be excluded from the EWAS). Alternatively, you can exclude CpGs identified as strongly associated with smoking from previous literature prior to or post EWAS (prior is recommended if using this option). Lastly, surrogate variables can capture latent variation due to smoking if you don't have a smoking variable included in the null model, but use this option with caution as the surrogate variables may attenuate more biological signal than intended. As a last resort, you can run the EWAS without accounting for smoking, just keep this in mind while interpreting results. ::: ```{r} #| label: ewas-fit # Array position is an 8-level factor on an EPIC chip, so it costs 7 dummy # columns. Rebuilding the factor from as.character() drops any level the sample # sheet defines but no modeled sample occupies; an empty level would add an # all-zero dummy column and leave the design rank-deficient. pos <- factor(as.character(mdk$array_pos)) # Every model term goes into one data frame so that model.matrix() can name the # columns it builds: `propk` expands to the six cell proportions by their own # names, and `SV` to SV1...SV6. check.names = FALSE keeps those names intact. dd <- data.frame(ptsd = ptsd, sex = factor(mdk$sex), age = mdk$age, smoke = smoke, pos = pos, propk, SV = SV, check.names = FALSE) design <- model.matrix(~ ptsd + sex + age + smoke + pos + CD8T + CD4T + NK + Bcell + Mono + Neu + SV, data = dd) stopifnot("ptsdCase" %in% colnames(design)) ew$design_ncol <- ncol(design) ew$resid_df <- ncol(Mcb) - ncol(design) # 756,251 least-squares fits in a single call -- which is why limma is the tool # here rather than a loop over lm(). This is the expensive step in the chapter: # roughly a minute, and several GB while the per-probe coefficients and residual # variances are held in memory. fit <- lmFit(Mcb, design) # per-probe weighted least squares fit <- eBayes(fit) # empirical-Bayes moderation of the variance # sort.by = "P" fixes the row order of the saved table to increasing p-value; # topTable()'s default sorts by the log-odds column B instead. # # limma:: is required, not decoration. bacon exports its own topTable method, # and because bacon is attached after limma above, a bare topTable() call # resolves to bacon's -- whose signature is (object, number, adjust.method, # sort.by) and which fails with "unused argument (coef = ...)". Qualifying the # call is the fix that does not depend on library() order. tt <- limma::topTable(fit, coef = "ptsdCase", number = Inf, sort.by = "P") tt$probe <- rownames(tt) # Delta-beta: the case-minus-control difference on the beta scale, taken from # the same ComBat-adjusted M-values the model saw, so the effect reported and # the effect tested describe one matrix. The back-transform is a second copy of # 756,251 x 87 doubles, so it is dropped again as soon as the means are taken. bcb <- 2^Mcb / (2^Mcb + 1) db <- rowMeans(bcb[, is_case, drop = FALSE], na.rm = TRUE) - rowMeans(bcb[, !is_case, drop = FALSE], na.rm = TRUE) tt$delta_beta <- db[tt$probe] rm(bcb); invisible(gc()) # Genomic inflation factor from the p-values: the conventional definition, and # the same one chapter 05's k-selection sweep used, so the two are comparable. lam <- function(p) median(qchisq(1 - p, 1), na.rm = TRUE) / qchisq(0.5, 1) lambda <- lam(tt$P.Value) # Fix the column order of the table that gets saved and deposited. tt <- tt[, c("probe", "logFC", "delta_beta", "AveExpr", "t", "P.Value", "adj.P.Val", "B")] ``` `limma`'s `eBayes` step is what makes this appropriate for an EWAS on a modest sample size: it moderates each probe's variance toward a pooled estimate, stabilizing the t-statistics for the ~700k probes tested [@ritchie2015limma]. ```{r} #| label: model-summary cat(sprintf("EWAS sample size: n = %d (%d cases, %d controls)\n", ew$n, ew$ncase, ew$nctrl)) cat(sprintf("Design columns (covariates + intercept): %d\n", ew$design_ncol)) cat(sprintf("Residual degrees of freedom per probe: %d\n", ew$resid_df)) cat(sprintf("Surrogate variables retained: %d\n", ew$k_sv)) cat(sprintf("CpGs tested: %s (%d smoking-panel probes excluded)\n", format(ew$n_tested, big.mark = ","), ew$panel_excluded)) ``` ## 3. Is the model well-calibrated? The QQ plot and genomic inflation {#sec-calibration} Before looking at *which* CpGs are top-ranked, check whether the p-values behave. The **genomic inflation factor** λ compares the median observed test statistic to its null expectation. λ ≈ 1 means the model is well-calibrated; λ ≫ 1 signals residual confounding (unmodeled batch, cell composition, or population structure) inflating the statistics; λ ≪ 1 can signal over-correction, but may also indicate that there truly is very little differential methylation or the sample size is underpowered to detect true signals. Due to the biological differences between genetic variants and methylation and the technical design differences in how they are measured in arrays, EWAS results tend to have inflated λ compared to GWAS results. Additional bias and inflation correction is performed post-EWAS to account for this. The lambda post-bias-and-inflation correction is what should be used to judge the model fit. ## 4. Correcting residual bias and inflation with BACON {#sec-bacon} **BACON** [@vanIterson2017bacon] estimates the *empirical null distribution* of the test statistics with a Bayesian three-component mixture (one null component plus two alternative components for positive and negative effects), then rescales every statistic to that estimated null. This corrects both **bias** (a non-zero center of the null — a systematic shift) and **inflation** (a null wider than N(0,1) — the same phenomenon λ measures), and it does so without assuming that all deviation from the diagonal is confounding. Always run BACON on your EWAS results: - A λ near 1 does not mean there is *no* residual bias — bias (a shift of the null's center) barely moves λ, which is driven by the *spread* of the statistics. BACON estimates and removes both. - Applying the empirical-null correction and *inspecting its diagnostics* is what lets you distinguish "well-calibrated" from "looks calibrated by one number." BACON works from either t-statistics *or* effect sizes and their standard errors. We chose effect sizes and standard errors as we can derive these from the output of limma and these estimates are needed down stream if a meta-analysis is performed. ```{r} #| label: run-bacon # The Gibbs sampler is stochastic, so the seed is what makes the inflation, # the bias and every adjusted p-value below reproducible run to run. set.seed(42) es <- tt$logFC # effect size per CpG se <- tt$logFC / tt$t # standard error = beta / t-statistic # 5,000 iterations with 2,000 burn-in (bacon's defaults) over all 756,251 # statistics: about a minute. Passing NULL for the test statistics tells bacon # to work from the effect size / standard error pair instead. bc <- bacon(NULL, effectsizes = es, standarderrors = se) # es() and se() below are bacon's accessor functions, not the two numeric # vectors of the same name defined above: R resolves a name used in function # position to a function and skips the vector. inflation(bc) # estimated inflation (sigma.0); 1 = none bias(bc) # estimated bias (mu); 0 = none tt$bacon.p <- pval(bc) # empirical-null-adjusted p-values tt$bacon.es <- es(bc) # adjusted effect sizes tt$bacon.se <- se(bc) # adjusted standard errors # BH across the adjusted p-values. The results section thresholds on this # column, so it has to be computed from bacon.p rather than reused from limma. tt$bacon.adj.P <- p.adjust(tt$bacon.p, method = "BH") ``` ```{r} #| label: bacon-summary # The four calibration numbers, collected into the object saved at the end of # the chapter as data/06_bacon_summary.rds. lambda_raw is limma's, lambda_bacon # is the same statistic recomputed on the empirical-null-adjusted p-values. bs <- list(inflation = inflation(bc), bias = bias(bc), lambda_raw = lambda, lambda_bacon = lam(tt$bacon.p), n = ew$n) cat(sprintf("Estimated inflation (sigma.0): %.3f\n", bs$inflation)) cat(sprintf("Estimated bias (mu): %.3f\n", bs$bias)) cat(sprintf("lambda before BACON: %.3f\n", bs$lambda_raw)) cat(sprintf("lambda after BACON: %.3f\n", bs$lambda_bacon)) ``` Here the model was already well-calibrated (λ ≈ `r sprintf("%.2f", bs$lambda_raw)`), so BACON changes little — the estimated inflation is ≈ 1 and the adjusted λ stays close to 1. That is the *reassuring* outcome: the correction is available and applied, and it confirms rather than rescues the analysis. On a poorly-specified model you would instead see inflation ≫ 1 pulled back toward 1, and the diagnostics below would show it. ### Always inspect the four BACON diagnostic plots The single number can hide a mis-specified or non-converged fit. BACON's Gibbs sampler produces four diagnostic plots, and you should look at all four every time. (The interpretation below follows the guidance in the [pipeline repository](https://github.com/krferrier/EWAS) README, where these are documented in full.) ```{r} #| label: bacon-diagnostics # The four plots below are drawn from `bc` and written to disk here, so the # figures that follow display this run rather than a stored copy of an older # one. They are base and lattice graphics rather than ggplots, so each goes to # its own jpeg() device instead of through ggsave(). dir.create("data/06_bacon", showWarnings = FALSE, recursive = TRUE) # fit() is bacon's mixture-fit plot. The limma object named `fit` above does not # get in the way, for the same reason es() and se() still worked: a name in # function position resolves to a function. jpeg("data/06_bacon/fit.jpg", width = 1600, height = 1200, res = 200) print(fit(bc, n = 100)) # n = 100 histogram bins dev.off() jpeg("data/06_bacon/qq.jpg", width = 1600, height = 1200, res = 200) print(plot(bc, type = "qq")) # a ggplot, hence the print() dev.off() # traces() and posteriors() draw straight onto the device, so no print() here. # burnin = FALSE shows the chain after burn-in only, which is the part you # judge convergence on. jpeg("data/06_bacon/traces.jpg", width = 1600, height = 1400, res = 180) traces(bc, burnin = FALSE) dev.off() jpeg("data/06_bacon/posteriors.jpg", width = 1600, height = 1200, res = 200) posteriors(bc) dev.off() ``` ```{r} #| label: bacon-fit #| fig-cap: "**Fit plot.** Histogram of the observed z-scores with the fitted three-component mixture overlaid: the black curve is the overall fit, red is the estimated null component, and blue/green are the two alternative components. The black curve tracks the histogram closely — the model is well-specified. A large gap between the fitted (black) and observed (histogram) densities would signal mis-specification or non-convergence." #| echo: false knitr::include_graphics("data/06_bacon/fit.jpg") ``` ```{r} #| label: bacon-qq #| fig-cap: "**QQ plots, before (uncorrected) and after (corrected) BACON.** Points should hug the diagonal, breaking upward into a tail spike only for genuine associations. Systematic departure from the line *before* the spike is bias/inflation; after correction the bulk of points should sit closer to the diagonal. Here the two panels are nearly identical because the input was already calibrated." #| echo: false knitr::include_graphics("data/06_bacon/qq.jpg") ``` ```{r} #| label: bacon-traces #| fig-cap: "**Traces plot.** The Gibbs sampler's value for each parameter (p, mu, sigma for each of the three components) across iterations. A well-mixed, converged chain looks like a fuzzy horizontal band with a narrow range in estimates — a 'hairy caterpillar' — with no drift or trend. A sharp early sweep followed by settling is normal (the prior was far from the estimate); a chain that never settles means non-convergence." #| echo: false knitr::include_graphics("data/06_bacon/traces.jpg") ``` ```{r} #| label: bacon-posteriors #| fig-cap: "**Posteriors plot.** Scatter of the sampled posterior draws for the null inflation (sigma.0) and the null proportion (p.0), with 75/90/95% probability ellipses. A dense, roughly elliptical cloud sitting inside the contours indicates the sampler converged on a stable estimate; a sparse or lopsided cloud spilling outside the ellipses indicates it did not." #| echo: false knitr::include_graphics("data/06_bacon/posteriors.jpg") ``` In this run the fit tracks the histogram, the QQ panels are near-diagonal, the traces are hairy caterpillars, and the posteriors form a dense cloud inside the ellipses — all four say the sampler converged and the empirical null is trustworthy. ::: {.callout-tip} ## Make BACON a default, and read the diagnostics Report BACON-adjusted effect sizes and p-values as your primary EWAS results. It is not usually necessary to report the diagnostic plots, but I strongly recommend reporting the adjusted lambda. ::: ## 5. Results {#sec-results} ```{r} #| label: hits bonf <- 0.05 / nrow(tt) n_bonf <- sum(tt$bacon.p < bonf) n_fdr <- sum(tt$bacon.adj.P < 0.05) cat(sprintf("Probes tested: %s\n", format(nrow(tt), big.mark = ","))) cat(sprintf("Bonferroni-significant (p < %.2g): %d\n", bonf, n_bonf)) cat(sprintf("FDR < 0.05: %d\n", n_fdr)) knitr::kable(head(tt[order(tt$bacon.p), c("probe","logFC","delta_beta","P.Value","bacon.p","bacon.adj.P")], 6), digits = c(0, 3, 4, 8, 8, 4), row.names = FALSE, col.names = c("probe", "logFC", "delta_beta", "P (limma)", "P (BACON)", "FDR (BACON)"), caption = "Top 6 CpGs, ranked by BACON-adjusted p-value. logFC is the M-value effect; delta_beta is the case−control β difference (interpretable as percent methylation).") ``` Only `r n_bonf` CpG passes Bonferroni and `r n_fdr` pass FDR < 0.05 — a borderline result that a study of this size cannot settle either way. With n = 87, a single marginal hit is not unexpected. The next two figures are the same Manhattan plot drawn twice, once on limma's uncorrected p-values and once on the BACON-adjusted ones. Plotting both is worth the space: it shows exactly what the correction did to your results, and it makes clear that the one CpG which clears Bonferroni does so *because of* the rescaling, not independently of it. Report the adjusted version. ```{r} #| label: manhattan #| fig-cap: "Manhattan plot of the EWAS, drawn on limma's uncorrected p-values. Nothing crosses the Bonferroni line (red) on this scale — the single CpG that clears Bonferroni does so only after BACON rescaling, which shifts it just past the threshold. A near-empty Manhattan is the expected result for a subtle psychiatric phenotype in a small teaching subset (n=87); the full-scale analysis (n≈800) and meta-analysis across cohorts is what powers genome-wide discovery for PTSD [@katrinli2020ptsd]." #| echo: false knitr::include_graphics("data/06_manhattan.png") ``` ```{r} #| label: manhattan-bacon #| fig-cap: "The same Manhattan plot drawn on the BACON-adjusted p-values — the ones the chapter recommends you report. The axes, the coordinates, and the Bonferroni line are identical to the previous figure; only the y-values change. Rescaling by the empirical null lifts `cg15434749` (filled point) just across Bonferroni and leaves `cg20600436` (open ring) passing FDR < 0.05 alone. Comparing the two panels is the point: BACON does not manufacture signal, it re-centres and re-scales the whole test-statistic distribution, and a CpG that was already at the top of the uncorrected plot is the one that moves across." #| echo: false knitr::include_graphics("data/06_manhattan_bacon.png") ``` ```{r} #| label: volcano #| fig-cap: "Volcano plot. The x-axis is the EWAS effect size itself — the `ptsdCase` coefficient from the limma fit, on the M-value scale, adjusted for every covariate in the design. Colour marks the sign of that coefficient: teal for CpGs less methylated in cases, plum for those more methylated. Note that this is the adjusted effect, so its sign can differ from a raw case-minus-control β difference for probes where the covariates carry part of the signal. The three most significant probes are ringed and labelled. Even they move methylation only slightly, illustrating that EWAS effect sizes for complex phenotypes are typically small and require large samples to detect reliably." #| echo: false knitr::include_graphics("data/06_volcano.png") ``` ## 6. An alternative in common use: the M-value residual workflow {#sec-residual} The design above puts **all** covariates — biological and technical — into one large regression. A common alternative, and the convention many groups (including this tutorial's pipeline author) prefer, is a **two-stage residualization**: 1. **Regress the technical variation out first.** Fit each CpG's M-value on only the technical/nuisance terms (surrogate variables, or recorded batch such as chip and position) and keep the **residuals** — the methylation signal with technical structure removed. 2. **Run the association on the residuals** with only the biological model: $$\text{resid}(\text{M})_{\text{CpG}} \sim \text{PTSD} + \text{sex} + \text{age} + \text{cell composition}$$ ```r # Stage 1: residualize M-values on the technical terms (SVs here; or chip + position) tech <- model.matrix(~ SV) # nuisance-only design M_resid <- residuals(lmFit(M, tech), M) # limma residuals, probes x samples # Stage 2: association model on the residuals, biological covariates only bio <- model.matrix(~ ptsd + sex + age + Neu + NK + CD4T + CD8T + Bcell + Mono, data = md) fit <- eBayes(lmFit(M_resid, bio)) ``` **Why do this?** The usual motivation is **degrees of freedom**. Our single-stage model spends `r ew$design_ncol` design columns *per probe*; at n = `r ew$n` that leaves `r ew$resid_df` residual df. Moving the technical terms into a first stage — where they are estimated once — leaves the association model lean (`PTSD + sex + age + 6 cell props`) and returns df to the test that matters. It is also common to regress out the blood cell proportions in addition to the technical variation prior to running the main model. There is no *wrong* choice in which or how many covariates you adjust for prior to the EWAS, just make sure to report what exactly was done in your methods. ::: {.callout-warning} ## The caveat: residualization understates uncertainty Two-stage residualization treats the stage-1 fitted values as *known* when it runs stage 2, so it doesn't propagate the uncertainty of estimating the technical terms. In practice, with many samples relative to nuisance terms, the difference from the single-stage model is small and the df savings are worth it — which is why it is widely used. But be aware it is an approximation to the joint model, not identical to it. Report which approach you used, and ideally confirm the top hits are stable under both. ::: ## 7. When and why to stratify, then meta-analyze {#sec-stratify} The model so far adjusts for sex as a **covariate**: it lets men and women have different mean methylation, but forces the *PTSD effect* to be a single number shared across both. That assumption is often wrong for methylation. ### Why sex is a special case in methylation Sex is one of the largest sources of variation in the methylome, and it is not just an X/Y-chromosome effect — autosomal CpGs show widespread sex differences in mean level, in variance, and in how they *respond* to exposures and phenotypes. When the effect of your exposure/trait/disease genuinely differs between men and women (a sex interaction), a sex-adjusted model estimates a blend of the two and can miss a real signal that is present in only one sex, or that points in opposite directions. Because sex-specific methylation effects are common, I strongly recommend performing sex-stratified analyses. More generally, stratify when you have reason to believe the **effect itself** (not just the baseline) differs across a grouping variable: genetic ancestry, tissue or cell subtype, exposure dose, or age categories can all justify it. The cost is power — each stratum has fewer samples — so stratification trades a cleaner, assumption-light estimate within each group against reduced within-stratum sample size. However, the loss of power from sample size can be recovered with meta-analysis to re-join the strata results. ### Stratify → analyze → meta-analyze The workflow is three steps: 1. **Split** the data by the stratifying variable (e.g. sex) and drop that variable from the covariate set within each stratum — it is now constant. 2. **Run the full EWAS separately** in each stratum: the same model, the same BACON adjustment, the same diagnostics, per stratum. 3. **Meta-analyze** the per-stratum results with an inverse-variance-weighted fixed-effect meta-analysis, which combines the two effect estimates weighting each by its precision (1/SE²). This recovers power from both strata while still allowing the underlying effects to differ — and it gives you a **heterogeneity test** (Cochran's *Q*, *I²*) that formally asks whether the sexes differ at each CpG. Inverse-variance meta-analysis of a single CpG is short enough to show by hand: ```r # per-CpG effect + SE from each stratum's BACON-adjusted results # (bF, seF) from females; (bM, seM) from males wF <- 1 / seF^2 # inverse-variance weights wM <- 1 / seM^2 beta_meta <- (wF*bF + wM*bM) / (wF + wM) # combined effect se_meta <- sqrt(1 / (wF + wM)) # combined standard error z_meta <- beta_meta / se_meta p_meta <- 2 * pnorm(-abs(z_meta)) # heterogeneity: does the effect differ between the sexes? Q <- wF*(bF - beta_meta)^2 + wM*(bM - beta_meta)^2 # Cochran's Q, 1 df p_het <- pchisq(Q, df = 1, lower.tail = FALSE) ``` In practice you run this across all CpGs at once. The standard tool for it is **METAL** [@willer2010metal] — the same command-line meta-analyzer used for GWAS. It reads each stratum's results file, matches CpGs by ID, does the inverse-variance weighting and the heterogeneity test genome-wide, and writes a combined table. A concordant hit (same direction in both sexes, low heterogeneity) is far more credible than a hit seen in only one stratum. [Chapter 07](07_pipeline.qmd) runs the stratified EWAS and METAL as a managed workflow. ### You do not need METAL for this METAL's `SCHEME STDERR` *is* inverse-variance-weighted fixed-effects meta-analysis. The single-CpG code above vectorizes directly, so the whole genome-wide meta-analysis is a dozen lines of R over the two strata's BACON-adjusted result tables: ```r library(data.table) # each stratum's BACON-adjusted results: one row per CpG # Written by the stratified pipeline run in chapter 07; paths are relative to # the tutorial folder. If you have not run the pipeline, `./get_data.sh # G_pipeline_run` fetches the published run these two files come from. mF <- fread("../ewas_pipeline/run_grady/F/F_PTSD_ewas_bacon_results.csv.gz")[ , .(probe = cpgid, bF = bacon.es, seF = bacon.se)] mM <- fread("../ewas_pipeline/run_grady/M/M_PTSD_ewas_bacon_results.csv.gz")[ , .(probe = cpgid, bM = bacon.es, seM = bacon.se)] mt <- merge(mF, mM, by = "probe") # CpGs tested in both strata mt[, wF := 1 / seF^2][, wM := 1 / seM^2] # precision weights mt[, beta_meta := (wF * bF + wM * bM) / (wF + wM)] # combined effect mt[, se_meta := sqrt(1 / (wF + wM))] # combined SE mt[, z_meta := beta_meta / se_meta] mt[, p_meta := 2 * pnorm(-abs(z_meta))] # heterogeneity, per CpG mt[, Q := wF * (bF - beta_meta)^2 + wM * (bM - beta_meta)^2] # Cochran's Q, 1 df mt[, p_het := pchisq(Q, df = 1, lower.tail = FALSE)] mt[, I2 := pmax(0, (Q - 1) / Q) * 100] # direction string, the way METAL reports it mt[, direction := paste0(ifelse(bF > 0, "+", "-"), ifelse(bM > 0, "+", "-"))] setorder(mt, p_meta) ``` This is not an approximation. Run against the same two strata this tutorial meta-analyzes with METAL, across all 756,251 CpGs tested in both, it returns the **identical** answer: the same top hit (`cg22671410`, *p* = 2.027 × 10⁻¹³, direction `+-`, *I²* = 0), the same **9** CpGs below the Bonferroni threshold of 6.61 × 10⁻⁸, the same 50 at FDR < 0.05, and a rank correlation of 1 between the two sets of p-values. METAL writes its effect and standard error columns rounded to four decimal places, and once you round the R result the same way, the largest disagreement across all 756,251 CpGs is **exactly zero**. ::: {.callout-tip} ## So when is METAL worth it? Not for the arithmetic — for the bookkeeping. METAL handles allele/strand alignment (irrelevant for methylation, central for GWAS), reads many cohort files with inconsistent column names, tracks sample sizes per marker, and applies genomic-control correction on request. For a two-stratum methylation meta-analysis where you control both input files, the R above is easier to audit and keeps everything in one language. For a multi-cohort consortium meta-analysis, use METAL. ::: ## 8. When you don't have IDATs {#sec-no-idats} Everything in this notebook operates on a **methylation matrix** and a phenotype table, so the non-IDAT path is identical once you have M-values: - Convert a deposited β-matrix to M-values (`M = log2(beta/(1-beta))`; see [normalization](02_normalization.qmd)). - Obtain covariates from the series metadata: sex and age from `pData`, cell composition via a **reference-based estimator that runs on β-matrices** (e.g. `EpiDISH`/`RPC`) when you can't run `estimateCellCounts2` from IDATs ([chapter 04](04_cell_composition.qmd) covers this), and surrogate variables from `sva` run on the matrix ([chapter 05](05_batch_effects.qmd)). - Fit the identical `limma` model. The association step never touches raw intensities. ## Saving what we just computed The fit and the BACON correction are the whole cost of this chapter. Three files carry them forward: the object [chapter 08](08_annotation.qmd) reads, the calibration summary section 4 reports from, and a flat copy of the results table. ```{r} #| label: save-ewas # 1. The handoff object. `tt` is attached only now, after BACON, so that the # saved table already carries bacon.p / bacon.es / bacon.se / bacon.adj.P -- # chapter 08 reads those columns straight out of this file and never reruns # the sampler. ew$tt <- tt ew$lambda <- lambda saveRDS(ew, "data/06_ewas.rds") # 2. The four calibration numbers, on their own, for anything that needs them # without loading a 756,251-row table. saveRDS(bs, "data/06_bacon_summary.rds") # 3. A flat copy of the results table, one row per CpG. Gzipped CSV rather than # a tab-delimited file in the repository: 756,251 rows is far too large to # commit, so this one is deposited in the Zenodo record instead. The write # takes about a minute. write.csv(tt, gzfile("data/06_ewas_bacon_toptable.csv.gz"), row.names = FALSE) ``` ::: {.callout-note collapse="true"} ## Picking the chapter up from the saved fit instead Fitting 756,251 probes needs the 510 MB ComBat matrix in memory with a few GB of working space on top of it. If that is not practical on your machine, fetch the published fit and read it instead of running the two chunks above: ```{.bash filename="Terminal"} ./get_data.sh F_ewas_results ``` ```{.r filename="RStudio Console"} ew <- readRDS("data/06_ewas.rds") tt <- ew$tt # already carries the BACON columns bs <- readRDS("data/06_bacon_summary.rds") ``` Every number and table in this chapter will then run. The four diagnostic plots will not, because they are drawn from the `bc` object rather than from the saved table. ::: --- **Next:** [Running the EWAS pipeline](07_pipeline.qmd) — scale this analysis with a reproducible Snakemake workflow that models methylation as the outcome, runs a **sex-stratified** EWAS, and meta-analyzes the strata with METAL. Then [functional annotation](08_annotation.qmd) of the results. ## References {.unnumbered} ::: {#refs} :::