Methylation: scales, mixtures, and regional targets

What a CpG summary does and does not identify

Before reading: The Understanding an omics dataset page, proportions, logarithms, and elementary regression.

1 What does a methylation feature represent?

DNA methylation assays measure chemical modifications of DNA, often summarized at CpG sites, positions where a cytosine is followed by a guanine along a DNA strand. A site or region is the feature; the specimen supplies DNA from many molecules and cells. A methylation measurement is not a transcript count. Mapping a site to a nearby gene identifies an annotation relationship, not evidence that methylation caused a change in that gene’s expression.

On an array, methylated and unmethylated channel intensities contribute to a site-level summary. Bisulfite sequencing instead provides evidence from sequenced DNA molecules at covered positions. Standard bisulfite treatment does not distinguish 5-methylcytosine from 5-hydroxymethylcytosine, so distinguishing those modifications needs additional chemistry or information. Record assay chemistry, genome build, probe/site IDs, material, preprocessing, and coverage or detection-quality information.

Even an ideal molecular fraction is not automatically a fraction of cells. A value of one half could reflect methylation on one of two alleles in each cell, a mixture of differently methylated cells, or other mixtures of molecular states. Array response and preprocessing add further measurement considerations.

The Bioconductor methylation workflow covers array QC, probe filtering, normalization, cell-composition assessment, and site/region analyses. This primer concentrates on the statistical quantities produced by that process.

2 Beta and M are different outcome scales

Ignoring offsets for the moment, write the methylated and unmethylated intensities as I_M,I_U>0:

\begin{aligned} \beta&=\frac{I_M}{I_M+I_U},\\ M&=\log_2(I_M/I_U)\\ &=\log_2\{\beta/(1-\beta)\}. \end{aligned}

Real software can use offsets or additional preprocessing; record the actual formula. Beta values provide a bounded assay summary. They are not guaranteed to be exact fractions of methylated cells or molecules without a measurement argument. M values use a log-ratio scale. Du et al. compare the scales; a transformation’s variance properties do not make its scientific target interchangeable with a beta-scale target.

b <- seq(0.005, 0.995, length.out = 300)
ToM <- function(b) log2(b/(1-b))
par(mar = c(4, 4, 1, 1))
plot(b, ToM(b), type = "l", lwd = 2, xlab = "Beta summary", ylab = "M summary")
for (s in list(c(0.1, 0.2), c(0.5, 0.6))) {
  segments(s[1], ToM(s[1]), s[2], ToM(s[2]), col = "#D55E00", lwd = 4)
  points(s, ToM(s), pch = 16, col = "#D55E00")
}
data.frame(beta0 = c(0.1, 0.5), beta1 = c(0.2, 0.6),
           delta_beta = 0.1, delta_M = ToM(c(0.2, 0.6)) - ToM(c(0.1, 0.5)))

  beta0 beta1 delta_beta   delta_M
1   0.1   0.2        0.1 1.1699250
2   0.5   0.6        0.1 0.5849625

Figure 1: The idealized beta-to-M relation without offsets. Equal beta changes correspond to different M changes near the boundary and in the middle. No sample or significance test is represented.

A difference in E[\beta] has beta-scale units; 0.10 is ten percentage points, rather than a universal 10% relative increase. A change from 0.1 to 0.2 doubles the beta summary, while 0.5 to 0.6 increases it by 20% relative to baseline. On the M scale, the corresponding changes are \log_2(2.25)\approx1.17 and \log_2(1.5)\approx0.585.

The derivative dM/d\beta=1/\{\ln(2)\beta(1-\beta)\} grows near zero and one. Thus equal beta changes correspond to different M changes depending on the starting value. An M-scale coefficient has no single conversion to a beta difference independent of baseline and covariates.

For a beta-scale standardized target, model the beta outcome or justify a model for its conditional mean. If modeling the transformed distribution instead, its inverse is b(m)=2^m/(1+2^m). The needed mean is E[b(M)\mid A,Z], generally different from b(E[M\mid A,Z]). The former averages transformed individual values; the latter transforms one average. Flexible nuisance regression does not remove this distinction.

beta_values <- c(0.1, 0.2, 0.8)
m_values <- ToM(beta_values)
c(mean_beta = mean(beta_values), inverse_mean_M = plogis(mean(m_values) * log(2)))

     mean_beta inverse_mean_M 
     0.3666667      0.3246665 

3 Bulk tissue combines cells as well as molecular states

A bulk measurement pools DNA contributions from cell types. In an ideal two-component model, let w be the fraction of assayed DNA contributed by cell type 1, with fixed methylation fractions b_1,b_2. Then bulk fraction is w b_1+(1-w)b_2. These are DNA-contribution weights, which need not equal cell-number proportions if DNA content or recovery differs. The linear mixture illustrates a molecular fraction; real array responses need an additional measurement argument.

w <- seq(0, 1, length.out = 100)
bulk <- function(w) w * 0.8 + (1 - w) * 0.2
plot(w, bulk(w), type = "l", lwd = 2, col = "#2F80C8",
     xlab = "Cell type 1 fraction of assayed DNA", ylab = "Ideal bulk fraction", ylim = c(0, 1))
points(c(0.25, 0.75), bulk(c(0.25, 0.75)), pch = 16, col = "#D55E00", cex = 1.4)
text(c(0.25, 0.75), bulk(c(0.25, 0.75)) + 0.08, c("Group 0", "Group 1"))

Figure 2: An ideal two-cell-type DNA mixture with fixed fractions 0.8 and 0.2. Changing the mixture alone changes the bulk fraction; neither cell type changes its methylation.

With b_1=0.8 and b_2=0.2, the bulk fraction is 0.2+0.6w. Weights 0.25 and 0.75 give 0.35 and 0.65, so the group contrast is 0.30 despite zero within-cell-type contrasts. The difference comes entirely from mixture weights.

A total bulk comparison, a contrast standardized to common composition, and a cell-type-specific contrast are different estimands. Adjusting composition can remove part of a bulk difference the researcher intended to describe. Estimated proportions are learned quantities, not exact covariates by default; cell-type-specific identification needs information and assumptions separating the components, beyond adding a covariate to bulk regression.

4 Coverage and the unit of replication

In an ideal sequencing model, a specimen has methylated-read count K out of coverage N>0, with K\mid N,b\sim\operatorname{Binomial}(N,b). Here b is its molecular methylation probability. Binomial moments give

\begin{aligned} E[K/N\mid N,b]&=b,\\ \operatorname{Var}(K/N\mid N,b)&=\frac{b(1-b)}{N}. \end{aligned}

At b=1/2, coverage 5, 20, and 100 gives conditional measurement SDs about 0.224, 0.112, and 0.050. These describe repeated measurement of a specimen with fixed b, not the SE of a group mean across people. With zero coverage the proportion is undefined, rather than zero methylation. Duplicate reads or dependence between molecules can also invalidate the simple binomial variance.

If people have different b values but the same fixed coverage N, total variance gives

\operatorname{Var}(K/N) =\operatorname{Var}(b)+\frac{E[b(1-b)]}{N}.

Increasing coverage shrinks the second term; the between-person term remains. More sequencing of existing people therefore has a different role from sampling additional independent people.

coverage <- c(5, 20, 100)
data.frame(coverage, conditional_SE_at_half = sqrt(0.5 * 0.5 / coverage))

  coverage conditional_SE_at_half
1        5              0.2236068
2       20              0.1118034
3      100              0.0500000

A read-weighted pooled proportion and an equally person-weighted mean can differ. For two people with b=0.2,0.8 and fixed coverages 10 and 90, the person-average probability is 0.5. The expected pooled fraction is (10\times0.2+90\times0.8)/100=0.74. Pooling weights the second person nine times as heavily. Define which weights answer the question. In a longitudinal study, keep all specimens from a person together in folds and resampling; see Dependent observations and cluster-aware inference.

5 Site-level and region-level questions need separate definitions

Before a region analysis, define boundaries, sites, weights, and summary. A fixed equal-site average gives each CpG the same contribution; pooling read counts weights sites by coverage, potentially differently in every person. These define different outcomes. An omnibus region test can ask whether any site changes: two sites with changes +0.1 and -0.1 have zero mean shift but a nonzero vector of changes. Its p-value is not an estimate of a common regional shift.

Selecting a region because it contains small site p-values and then treating its boundaries as prespecified can change calibration. Record the site/region family, boundary-selection rule, and multiplicity guarantee.

A QC and provenance record should include detection/coverage filters, problematic or ambiguous probes, genetic variation affecting probes, genome/annotation versions, batch and position, normalization, and sample exclusions. Check the group–batch design and the role of cell composition before calling a site a disease mechanism. The shared design-and-contrast framework in Differential expression and abundance applies to an explicitly defined outcome and design; assay-specific noise and missingness still need attention.

6 Exercises

The checkpoint is this assay specialization’s completion target. Additional practice is optional. Try the checkpoint before opening its hint or solution.

Checkpoint · METH-01 — Distinguish scale and mixture effects

Task type: Computation and interpretation.

Cell types have fixed methylation fractions 0.8 and 0.2. Their assayed-DNA mixture weights are 0.25 in group 0 and 0.75 in group 1. Compute the bulk contrast and the within-cell-type contrasts. Compare the M changes for beta 0.1→0.2 and 0.5→0.6. Explain why neither a bulk association nor an M coefficient automatically gives a cell-specific beta effect.

Ready to continue when: you can calculate scale and composition effects and state the extra information needed for a stronger biological interpretation.

Transfer: If high-methylation specimens receive more sequencing reads, how could pooled-read and person-weighted regional summaries differ?

Hint

Apply the weighted average first; then use the log2 odds transform to compare scales.

Solution outline — open after trying

Bulk values are 0.35 and 0.65, so the contrast is 0.30. Both cell-type-specific contrasts are zero. The two beta changes are both 0.10, but the M changes are log2(2.25)≈1.170 and log2(1.5)≈0.585. Bulk composition changes can explain the first difference; nonlinear transformation explains the second. Cell-specific and beta-mean targets need their own measurement and identification/modeling arguments.

Additional practice (optional)

1. Using the three beta values in the example, compare mean beta with the inverse of mean M and explain the discrepancy.
2. Define a fixed-region target with equal site weights and another with coverage weights. Describe which changes if coverage changes.
3. Draw the metadata table needed to distinguish person, specimen, array position, batch, and estimated cell composition.

7 Reading

• Bioconductor methylation-array workflow: assay-aware QC and analysis.
• Du et al., beta- and M-value comparison: scale and variance behavior.
• The behaviour of 5-hydroxymethylcytosine in bisulfite sequencing: an experimental explanation of the chemical distinction standard bisulfite sequencing does not resolve.