Differential Methylation Testing at Individual CpGs

Why Test Methylation at an Individual CpG?

A differential methylation test asks whether methylation at a particular CpG site differs between groups, after accounting for specified covariates. The target is a group contrast for one site, not a claim that an entire genomic region changes or that group membership causes a change. Since an array can measure hundreds of thousands of sites, each site-specific result is also one part of a larger testing problem.

The choice of response scale affects how an estimated difference should be read. A beta value is on a 0-to-1 scale and supports direct statements in percentage points. An M value is a log-base-2 transformation of the methylation odds; its linear-model coefficients are often convenient for modelling, but are not constant differences in beta values. The earlier material on beta values and M values defines these scales; the essential formulas are given here so that the interpretation is explicit.

Model, Contrast, and Assumptions

For each CpG, ordinary least squares estimates the model coefficients by minimizing the sum of squared residuals. Testing no adjusted group difference at site \(j\) means testing \(H_0:\beta_{1j}=0\) against, commonly, the two-sided alternative \(H_a:\beta_{1j}\ne 0\). More generally, the scientific comparison may be a linear contrast \(c^\top\boldsymbol{\beta}_j\), where \(\boldsymbol{\beta}_j\) is the coefficient vector and \(c\) specifies the comparison. For example, with treatment coding and no interaction, the group coefficient itself is the contrast.

The interpretation depends on the model and design. The mean response must be adequately described by the included terms; the design matrix must have full column rank so the contrast is estimable; and, for the usual small-sample \(t\) test and confidence interval to be exact, errors must be independent, normally distributed, and have a common variance conditional on the predictors. Approximate inference can be useful with larger samples, but serious departures, influential observations, or an inadequate mean model can still undermine it. Samples, rather than repeated technical measurements of the same sample, are the independent observational units in this basic model.

Covariates can include measured characteristics or estimated cell-type fractions. If cell fractions sum to one, including every fraction together with an intercept creates exact collinearity; use an identifiable parameterization, such as omitting a reference fraction. The fitted group contrast is conditional on the covariates included. It is not automatically a causal effect, and the model does not repair a poorly chosen or non-identifiable adjustment set.

For an estimable contrast \(\theta_j=c^\top\boldsymbol{\beta}_j\), its estimated standard error is obtained from the coefficient covariance matrix \(\widehat{\mathbf V}_j\):

Under the null and the model assumptions, this statistic follows a \(t\) distribution with residual degrees of freedom \(df\). A two-sided \(p\)-value is \(2P(T_{df}\ge |t|)\). A \(100(1-\alpha)\%\) confidence interval is \(\widehat{\theta}_j\pm t_{1-\alpha/2,df}SE(\widehat{\theta}_j)\). The interval describes the range of contrast values compatible with the data and model at that confidence level; it is not a probability statement about the fixed parameter after observing the data.

Worked Examples on the Two Scales

Contrasts When the Group Effect Varies

If the model includes a group-by-covariate interaction, the group coefficient alone is not the group difference at every covariate value. The relevant contrast must be evaluated at the covariate profile of interest, and its uncertainty must include covariance among the estimated coefficients.

Genome-Wide Testing, Diagnostics, and Alternatives

A test at one CpG has a sitewise \(p\)-value. When tests are performed at many CpGs, some small \(p\)-values can arise by chance even if all null hypotheses hold. A common approach is to apply the Benjamini–Hochberg procedure to the full set of sitewise \(p\)-values: order them from smallest to largest, compare the \(k\)th ordered value with \(k\alpha/m\), where \(m\) is the number of tests, and identify the largest \(k\) that meets the comparison. This controls the false discovery rate under its standard dependence conditions. Report the adjustment method and the set of tests to which it was applied. An adjusted result is not the probability that a particular site is truly differential.

Inspect diagnostics at representative sites and, where feasible, at sites driving important conclusions. Residual-versus-fitted plots can reveal nonconstant variance or a misspecified mean; normal quantile plots can reveal severe tail departures; leverage and influence checks can identify samples that dominate a contrast. Beta values are bounded, so a linear model can in principle predict outside 0 to 1, especially near the boundaries. M values remove that bound but do not guarantee normal, homoskedastic residuals. The scale should be chosen with both inference and biological interpretation in mind.

An unadjusted two-group comparison can be simpler, but does not estimate a covariate-adjusted contrast. A rank-based test can reduce reliance on a normal model, but generally targets a different comparison and does not automatically provide the same covariate-adjusted effect. A transformation or alternative response model may be worth considering when the linear-model diagnostics are poor; no method makes design limitations disappear. For the present model, a clear report gives the response scale, contrast, effect estimate and units, uncertainty, sitewise and adjusted evidence, and the covariates used.