limma versus DSS for methylation: a statistical comparison

statistics

Bioconductor

methylation

software

Linear-model empirical Bayes versus beta-binomial count modeling—and why raw EM-seq CpG evidence favors DSS in this study

Author

Methylation & Colitis team

Published

August 24, 2026

The conclusion first

limma and DSS both borrow strength across genomic features, but at different levels of the model. limma moderates feature-wise residual variances after fitting linear models to approximately continuous responses. DSS models methylated reads conditional on coverage with a beta-binomial distribution and shrinks CpG-wise biological dispersion. For the project’s raw EM-seq methylated/total counts, DSS preserves the sampling structure that limma would otherwise need to approximate.

Scientific illustration of molecular evidence entering statistical analysis — The appropriate statistical engine depends on whether the evidence arrives as normalized continuous measurements or methylated reads conditional on coverage.

What problem are limma and DSS solving?

In plain language, this is a group-comparison and feature-discovery analysis. It asks:

At each CpG, is methylation systematically different between the biological groups after accounting for other variables in the design?

For this study, a useful model might compare colitis, colorectal cancer, and controls while adjusting for factors such as age, sex, batch, or cell composition. The output is an effect estimate and uncertainty for every CpG. Nearby CpGs can then be joined into differentially methylated regions (DMRs).

It is easy to call the resulting loci a “fingerprint,” but that word can hide several distinct goals:

Goal	Question	What the result looks like	Do limma/DSS solve it alone?
Differential methylation	Which CpGs or regions differ between groups?	Effect sizes, p-values, FDR, DMLs and DMRs	Yes—this is their main job
Biological fingerprint discovery	Is there a reproducible pattern associated with colitis or cancer?	A candidate set of loci/regions and pathways	Partly; they nominate features, but replication is required
Classifier	Can an unknown sample be assigned to the correct group?	One prediction per sample, with sensitivity, specificity and ROC/PR performance	No; a prediction model and held-out validation are needed
Early detector	Can disease be detected before symptoms or routine diagnosis?	Performance in the intended screening population, ideally using pre-diagnostic samples	No; this is a clinical-use claim requiring a separate study

The correct label for this analysis

This limma-versus-DSS comparison is about differential methylation discovery. It can supply candidates for a colitis/cancer methylation fingerprint. It is not yet a classifier, and it cannot establish early detection from case–control samples collected after diagnosis.

How discovery could become a detector

Code

flowchart LR
  A[EM-seq reads] --> B[Filtered CpG counts]
  B --> C[DSS differential analysis]
  C --> D[Candidate CpGs and DMRs]
  D --> E[Locked feature panel]
  E --> F[Classifier trained without leakage]
  F --> G[Held-out external validation]
  G --> H[Prospective early-detection study]

  classDef now fill:#dcefee,stroke:#167d78,color:#123b5d
  classDef later fill:#f7efe4,stroke:#b8634f,color:#49352e
  class A,B,C,D now
  class E,F,G,H later

flowchart LR
  A[EM-seq reads] --> B[Filtered CpG counts]
  B --> C[DSS differential analysis]
  C --> D[Candidate CpGs and DMRs]
  D --> E[Locked feature panel]
  E --> F[Classifier trained without leakage]
  F --> G[Held-out external validation]
  G --> H[Prospective early-detection study]

  classDef now fill:#dcefee,stroke:#167d78,color:#123b5d
  classDef later fill:#f7efe4,stroke:#b8634f,color:#49352e
  class A,B,C,D now
  class E,F,G,H later

The current analysis occupies the first four boxes. To move farther right, feature selection must occur inside training resamples; otherwise, using all samples to find CpGs before cross-validation leaks outcome information and exaggerates accuracy. A convincing detector also needs independent samples representing the population and clinical moment in which the test would actually be used.

limma: feature-wise linear models with moderated variance

For feature \(g\), limma fits the familiar model

\[\mathbf y_g = X\boldsymbol\beta_g + \boldsymbol\varepsilon_g, \qquad \boldsymbol\varepsilon_g \sim N(0,\sigma_g^2 W_g^{-1}),\]

where \(X\) encodes groups and covariates, contrasts select estimands, and \(W_g\) may contain precision weights. The distinctive step is empirical-Bayes moderation of the feature-wise residual variances. A scaled inverse-chi-squared prior combines the noisy \(s_g^2\) with a shared prior variance \(s_0^2\) and prior degrees of freedom \(d_0\):

\[\widetilde s_g^2 = \frac{d_0s_0^2 + d_gs_g^2}{d_0+d_g}.\]

The moderated \(t\) statistic uses \(\widetilde s_g\) and augmented degrees of freedom. This stabilizes inference when each feature has few residual degrees of freedom, while preserving arbitrary design matrices and contrasts. Robust and trend options relax the single-prior assumption; voom estimates a mean–variance trend and observation-level weights for RNA-seq log-counts.

For methylation arrays, the common response is a normalized M-value, \(\log_2(M/U)\), which is generally more suitable for linear-model inference than the bounded beta value \(M/(M+U)\). minfi handles array-specific preprocessing; limma handles the design and differential analysis.

DSS: beta-binomial evidence with dispersion shrinkage

For CpG \(j\) in sample \(i\), DSS begins with the actual sequencing evidence:

\[M_{ij}\mid N_{ij},p_{ij},\phi_j \sim \operatorname{BetaBinomial}(N_{ij},p_{ij},\phi_j),\]

where \(M_{ij}\) is methylated reads, \(N_{ij}\) is informative coverage, \(p_{ij}\) is underlying methylation, and \(\phi_j\) represents extra-binomial biological dispersion. One useful parameterization gives

\[\operatorname{Var}(M_{ij}\mid N_{ij}) = N_{ij}p_{ij}(1-p_{ij})\{1+(N_{ij}-1)\phi_j\}.\]

The first component is read-sampling variance; the multiplier admits participant-to-participant heterogeneity beyond binomial sampling. DSS estimates locus-specific dispersions with empirical-Bayes shrinkage, then tests group or multifactor effects. Nearby locus-level evidence can be aggregated into DMRs using spatial and effect-size criteria.

Similar philosophy, different likelihood

limma

Response: continuous transformed feature

Borrowing: residual variances across features

Precision: residual variance and optional weights

Natural home: arrays and derived summaries

DSS

Response: methylated count given total coverage

Borrowing: beta-binomial dispersion across CpGs

Precision: coverage plus overdispersion

Natural home: BS-seq / EM-seq counts

Statistical issue	limma	DSS	Consequence here
Unequal CpG coverage	Must be represented through a transformed response and defensible precision weights	Enters directly as \(N_{ij}\)	DSS retains the observation mechanism
Bounded methylation proportion	Usually handled by M-value/logit-style transformation	Modeled through a count probability	DSS avoids treating raw percentages as homoscedastic Gaussian data
Extra biological variation	Moderated residual variance	Beta-binomial dispersion with shrinkage	Both borrow strength, but DSS does so in the count likelihood
Complex design	Excellent design/contrast machinery	Supports multi-factor designs, subject to sample size	Neither creates degrees of freedom in an 18-person cohort
Regional signal	Requires explicit aggregation or region-level inputs	DML statistics feed DMR calling	DSS offers a coherent CpG-to-region route
Array validation	Established choice with `minfi` preprocessing	Not designed for array intensities	limma remains preferred for 450K/EPIC external data

A real-data toy comparison in R

The example below uses the RRBS dataset bundled with DSS: 5,000 CpGs from 16 real reduced-representation bisulfite-sequencing samples. Its \(2\times2\) design contains healthy-control (HC) and systemic-lupus (SLE) samples across two cell states (aN and rN). The biological question is not colitis; the dataset is useful because it is small, public, count-based, and has a genuine nuisance factor to adjust.

Why this is relevant to EM-seq

RRBS and EM-seq use different laboratory chemistries and have different coverage profiles. After conversion-aware alignment and methylation extraction, however, both can provide the same statistical unit used here: methylated count \(M\) and total informative coverage \(N\) at each CpG. Replace RRBS with the project’s filtered EM-seq BSseq object and keep the design explicit; do not assume the toy result transfers biologically.

Load the counts and design

Code

suppressPackageStartupMessages({
  library(DSS)
  library(bsseq)
  library(limma)
})

data(RRBS, package = "DSS")

RRBS

An object of type 'BSseq' with
  5000 methylation loci
  16 samples
has not been smoothed
All assays are in-memory

Code

design

   case cell
1    HC   rN
2    HC   rN
3    HC   rN
4   SLE   aN
5   SLE   rN
6   SLE   aN
7   SLE   rN
8   SLE   aN
9   SLE   rN
10  SLE   aN
11  SLE   rN
12   HC   aN
13   HC   aN
14   HC   aN
15   HC   aN
16   HC   rN

Code

table(design$case, design$cell)

     
      aN rN
  HC   4  4
  SLE  4  4

The object preserves two matrices: methylated reads (M) and total coverage (Cov). This is the structure expected from an EM-seq caller such as MethylDackel after genomic coordinates and samples are assembled into a count-aware object.

Fit DSS to the counts

Code

invisible(capture.output(
  dss_fit <- DMLfit.multiFactor(
    RRBS,
    design = design,
    formula = ~ cell + case
  )
))

dss_result <- DMLtest.multiFactor(dss_fit, coef = "caseSLE")
head(dss_result)

   chr     pos       stat     pvals      fdrs
1 chr1 2771453  1.5123894 0.1304348 0.9837621
2 chr1 2771503  0.5213406 0.6021295 0.9937994
3 chr1 2771512 -0.1505977 0.8802931 0.9940470
4 chr1 2771609  1.5368674 0.1243258 0.9698905
5 chr1 2771635  0.6941054 0.4876161 0.9937994
6 chr1 2771655  0.2321348 0.8164333 0.9937994

DSS tests the disease coefficient while retaining coverage and estimating beta-binomial dispersion. The cell term adjusts for the second factor.

Fit limma to transformed proportions

To make a deliberately simple comparison, the same counts are converted to M-values with a 0.5 pseudocount. This avoids infinite logits at fully methylated or unmethylated loci, but it demonstrates the extra choice limma requires. This toy path is not presented as the optimal way to weight sequencing methylation data.

Code

methylated <- getCoverage(RRBS, type = "M")
coverage   <- getCoverage(RRBS, type = "Cov")

m_value <- log2(
  (methylated + 0.5) /
  (coverage - methylated + 0.5)
)

X <- model.matrix(~ cell + case, data = design)

limma_fit <- lmFit(m_value, X)
limma_fit <- eBayes(limma_fit, trend = TRUE, robust = TRUE)

limma_result <- topTable(
  limma_fit,
  coef = "caseSLE",
  number = Inf,
  sort.by = "none"
)

Put both answers side by side

Code

gr <- granges(RRBS)
limma_result$chr <- as.character(seqnames(gr))
limma_result$pos <- start(gr)

# A reader-friendly biological effect: mean methylation in SLE minus HC.
# The design is balanced across cell state, so this is useful for visualization;
# inference still comes from the adjusted models above.
beta <- methylated / coverage
delta_beta <- rowMeans(beta[, design$case == "SLE"], na.rm = TRUE) -
  rowMeans(beta[, design$case == "HC"], na.rm = TRUE)
limma_result$delta_beta <- delta_beta

comparison <- merge(
  dss_result[, c("chr", "pos", "stat", "pvals", "fdrs")],
  limma_result[, c("chr", "pos", "delta_beta", "logFC", "t", "P.Value", "adj.P.Val")],
  by = c("chr", "pos")
)

names(comparison)[3:10] <- c(
  "DSS_stat", "DSS_p", "DSS_FDR",
  "delta_beta", "limma_logFC", "limma_t", "limma_p", "limma_FDR"
)

comparison <- comparison[
  order(comparison$DSS_FDR, comparison$limma_FDR),
]

knitr::kable(
  head(comparison, 6),
  digits = 3,
  caption = "Six loci ranked first by DSS FDR; limma statistics are shown for the same CpGs."
)

Six loci ranked first by DSS FDR; limma statistics are shown for the same CpGs.
	chr	pos	DSS_stat	DSS_FDR	delta_beta	limma_logFC	limma_t	limma_p	limma_FDR
4053	chr1	3230622	-5.608	0.000	-0.372	-2.411	-4.434	0.000	0.084
4051	chr1	3230609	-4.705	0.006	-0.263	-1.625	-3.155	0.003	0.434
1334	chr1	2938265	4.406	0.018	0.374	4.182	5.475	0.000	0.041
4054	chr1	3230635	-4.079	0.057	-0.293	-2.242	-3.887	0.000	0.298
4553	chr1	3305518	3.956	0.064	0.361	2.931	4.689	0.000	0.053
201	chr1	2785745	3.995	0.064	0.268	1.684	3.205	0.002	0.430

Code

hit <- ifelse(
  comparison$DSS_FDR < 0.05 & comparison$limma_FDR < 0.05, "both",
  ifelse(comparison$DSS_FDR < 0.05, "DSS only",
         ifelse(comparison$limma_FDR < 0.05, "limma only", "neither"))
)

palette <- c(
  "neither" = "#a8b8ba",
  "DSS only" = "#167d78",
  "limma only" = "#d39a38",
  "both" = "#b8634f"
)

plot_data <- transform(
  comparison,
  delta_percent = 100 * delta_beta,
  DSS_evidence = -log10(pmax(DSS_p, .Machine$double.xmin)),
  call = factor(hit, levels = names(palette)),
  hover = sprintf(
    paste(
      "<b>%s:%s</b>",
      "Call: %s",
      "Mean methylation difference: %+.1f percentage points",
      "DSS statistic: %.2f | p: %.3g | FDR: %.3g",
      "limma logFC: %.2f | t: %.2f | p: %.3g | FDR: %.3g",
      sep = "<br>"
    ),
    chr, pos, hit,
    100 * delta_beta,
    DSS_stat, DSS_p, DSS_FDR,
    limma_logFC, limma_t, limma_p, limma_FDR
  )
)

plot_points <- function(p, x, y) {
  for (group in levels(plot_data$call)) {
    d <- plot_data[plot_data$call == group, ]
    emphasized <- group != "neither"
    p <- plotly::add_trace(
      p, data = d, x = x, y = y, name = group,
      text = ~hover, hoverinfo = "text",
      type = "scattergl", mode = "markers",
      marker = list(
        color = unname(palette[group]),
        size = if (emphasized) 10 else 5,
        opacity = if (emphasized) 0.9 else 0.22,
        line = list(width = if (emphasized) 1 else 0, color = "white")
      )
    )
  }
  p
}

volcano <- plot_points(
  plotly::plot_ly(), ~delta_percent, ~DSS_evidence
) |>
  plotly::layout(
    title = list(text = "DSS: biological effect versus evidence", x = 0.03),
    xaxis = list(title = "SLE − HC mean methylation (percentage points)", zeroline = TRUE),
    yaxis = list(title = "DSS −log10(p)"),
    legend = list(orientation = "h", x = 0, y = 1.12),
    shapes = list(list(
      type = "line", x0 = min(plot_data$delta_percent), x1 = max(plot_data$delta_percent),
      y0 = -log10(0.05), y1 = -log10(0.05),
      line = list(color = "#718386", dash = "dot")
    ))
  )

agreement <- plot_points(
  plotly::plot_ly(), ~DSS_stat, ~limma_t
) |>
  plotly::layout(
    title = list(text = "Do the methods agree on direction and strength?", x = 0.03),
    xaxis = list(title = "DSS test statistic", zeroline = TRUE),
    yaxis = list(title = "limma moderated t", zeroline = TRUE),
    showlegend = FALSE,
    shapes = list(
      list(type = "line", x0 = 0, x1 = 0,
           y0 = min(plot_data$limma_t), y1 = max(plot_data$limma_t),
           line = list(color = "#718386", dash = "dot")),
      list(type = "line", y0 = 0, y1 = 0,
           x0 = min(plot_data$DSS_stat), x1 = max(plot_data$DSS_stat),
           line = list(color = "#718386", dash = "dot"))
    )
  )

widget_config <- list(
  displaylogo = FALSE, responsive = TRUE,
  modeBarButtonsToRemove = c("lasso2d", "select2d")
)

Code

do.call(plotly::config, c(list(volcano), widget_config))

Figure 1

The x-axis is now interpretable as percentage-point methylation change. Colored, larger points survive 5% FDR in at least one method; quiet grey points provide context without hiding the signal.

Code

do.call(plotly::config, c(list(agreement), widget_config))

Figure 2

Points in the upper-right and lower-left quadrants agree in direction. Distance from the origin represents stronger standardized evidence—not methylation magnitude. The highlighted points make disagreements visible without forcing 5,000 loci into equally prominent marks.

Code

data.frame(
  quantity = c(
    "CpGs tested",
    "DSS FDR < 0.05",
    "limma FDR < 0.05",
    "Called by both",
    "Same effect direction",
    "Spearman correlation of raw p-values"
  ),
  value = c(
    nrow(comparison),
    sum(comparison$DSS_FDR < 0.05),
    sum(comparison$limma_FDR < 0.05),
    sum(comparison$DSS_FDR < 0.05 & comparison$limma_FDR < 0.05),
    sprintf("%.1f%%", 100 * mean(sign(comparison$DSS_stat) == sign(comparison$limma_t))),
    sprintf("%.3f", cor(comparison$DSS_p, comparison$limma_p, method = "spearman"))
  )
) |>
  knitr::kable(caption = "What happened in this one illustrative analysis.")

What happened in this one illustrative analysis.
quantity	value
CpGs tested	5000
DSS FDR < 0.05	3
limma FDR < 0.05	1
Called by both	1
Same effect direction	83.7%
Spearman correlation of raw p-values	0.632

In this run, DSS calls three CpGs at 5% FDR, limma calls one, and that locus overlaps. The signs agree for about 84% of CpGs, while the raw p-value rank correlation is about 0.63. These numbers are not a benchmark of universal power: they depend on the response transformation, pseudocount, weighting decision, dispersion model, design, coverage distribution, and this particular dataset. The lesson is that similar biological direction does not imply identical uncertainty.

Pros and cons in practice

Method	Pros	Cons	Best role in this project
DSS	Uses methylated and total counts directly; beta-binomial overdispersion; empirical-Bayes dispersion shrinkage; direct DML-to-DMR route	More specialized design/output; low coverage and tiny groups still limit inference; DMR settings can change results; standard workflow does not solve pairing or cell mixture automatically	Primary EM-seq CpG and regional inference, with a pre-specified multifactor design and coverage filters
limma	Excellent model matrices and contrasts; fast; mature diagnostics; robust empirical-Bayes variance moderation; ideal for array M-values	Raw sequencing precision is not intrinsic after transforming to a continuous response; pseudocounts and weights must be justified; no native CpG-count likelihood or automatic DMR construction	450K/EPIC validation, defensible region summaries, or a pre-specified sensitivity analysis

Do not choose the winner by p-value

The method should be selected from the assay and estimand before looking at which produces more discoveries. For this study, DSS matches raw EM-seq counts. limma remains valuable when the input is genuinely continuous or when a transparent secondary analysis has been specified in advance.

Why not simply run limma on methylation percentages?

Suppose two samples both have estimated methylation \(0.5\). One is based on \(M/N=1/2\) and the other on \(50/100\). A percentage-only model initially treats them as equal observations; their binomial standard errors are roughly \(0.354\) and \(0.050\). A transformation improves boundaries but cannot reconstruct uncertainty if \(N\) has been discarded.

One could construct coverage-informed weights and use limma on a transformed response. That can be reasonable for suitable summaries, and limma’s framework is flexible. It nevertheless adds modeling choices—transformation, pseudocounts, mean–variance estimation, and weights—when DSS already expresses the raw data-generating unit directly. With sparse or uneven EM-seq coverage, direct modeling is easier to defend.

Why DSS was the better primary choice here

The assay output matches its likelihood. MethylDackel produces methylated and total counts at CpGs.
Coverage varies. Plasma cfDNA and low-input libraries do not provide uniform depth across loci and participants.
Biological replicates are overdispersed. The beta-binomial separates read sampling from extra-binomial heterogeneity.
The target includes regions. The study needs interpretable DMRs, not only isolated CpG tests.
The audit trail is direct. Every estimate can be traced to methylated reads, total reads, dispersion, effect size, and regional criteria.

This preference is conditional, not a claim that DSS is universally more powerful. With six participants per group, dispersion and covariate estimates remain uncertain; smoothing and DMR thresholds can materially change results; cell mixture and batch can still confound effects; and asymptotic p-values do not substitute for independent replication.

Where limma still belongs in the project

External 450K/EPIC data: preprocess with minfi, analyze M-values with limma, and translate results across platforms cautiously.
Low-dimensional derived outcomes: use limma for approximately Gaussian region or pathway summaries when their construction and precision are defensible.
Factorial contrasts: exploit its transparent coefficient/contrast framework when the response suits a linear model.
Sensitivity analysis: agreement between a pre-specified continuous-summary analysis and the primary count model can test robustness, but should not be used to select the more favorable result.

Short conclusion: which should we use?

Starting data or goal	Use	Why
Raw EM-seq methylated and total reads	DSS	Coverage and beta-binomial biological variation belong in the likelihood
CpGs must be assembled into DMRs	DSS	It provides a coherent locus-to-region workflow
Normalized 450K/EPIC array M-values	limma	The transformed continuous response matches its linear-model strengths
Region/pathway summaries that are approximately Gaussian	limma	Its contrasts and empirical-Bayes variance moderation are excellent
Robustness check for the EM-seq result	Both, pre-specified	DSS remains primary; limma can test whether direction survives a different representation

Conclusion

For this study’s raw EM-seq counts, use DSS as the primary model. Use limma for compatible array data or justified continuous summaries, not as a contest to find smaller p-values. The toy data show why: the methods often agree on direction, but coverage and likelihood choices change uncertainty and therefore the final calls.

Read the broader software history → · See the beta-binomial explainer → · Follow the pipeline →