group1.mean and group2.mean do not fit with mean meth values per DMR (in-script)

[RicardoMartins-Ferreira]

Greetings,

I am analysing wgbs data from a total of 10 samples divided per 5 groups (n=2 per group). I calculated dmrs with the standart t.statistics for all pairwise comparisons.

For visualization in a heatmap, I generated a script (below) in which I generate a matrix with the mean methylation values of all CpGs within a dmr per sample. However, I noticed that the mean methylation values generated in my script do not match with the group.means in the t.stat output. In some cases, I see opposite behavior, a DMR that appears to be hypermethylated in one specific group in the heatmap, is actually hypomethylated in the same comparison in the t.stat results. Furthermore, I noticed that this behavior occurs often in DMRs with high intra-group variability.

Can you weight on this? Is this possible to happen with the t.stat method or am I missing some mistake in my workflow?

In-house script:

Make genomic ranges for bsseq object

ranges <- as.data.frame(bsseq.fit@rowRanges)
ranges$seqnames <- gsub("\.[0-9]+$", "", ranges$seqnames)
row.names(ranges) <- paste0(ranges$seqnames, "", ranges$start, "", ranges$end)
cpg_gr <- makeGRangesFromDataFrame(ranges,
keep.extra.columns = FALSE,
ignore.strand = TRUE,
seqnames.field = "seqnames",
start.field = "start",
end.field = "end",
strand.field = "strand",
starts.in.df.are.0based = FALSE)

Extract methylation matrix

Meth <- getMeth(bsseq.fit, type = "smooth") %>% as.data.frame()
colnames(Meth) <- bsseq.fit@colData$Sample
row.names(Meth) <- rownames(ranges)

Create unique DMR names

dmr_results$Name <- paste0(dmr_results$chr, "", dmr_results$start, "", dmr_results$end)
dmr_dict <- split(dmr_results, dmr_results$Name)
dmr_list <- names(dmr_dict)

Precompute overlaps

cpg_ids <- rownames(ranges)
cpg_gr$ID <- cpg_ids
dmr_gr <- makeGRangesFromDataFrame(dmr_results,
seqnames.field = "chr",
start.field = "start",
end.field = "end",
keep.extra.columns = TRUE,
ignore.strand = TRUE)

dmr_gr$Name <- dmr_results$Name

overlaps <- findOverlaps(dmr_gr, cpg_gr, ignore.strand = TRUE)
overlap_list <- split(subjectHits(overlaps), queryHits(overlaps))

Function using precomputed overlaps

calculate_dmr_means <- function(dmr_idx, dmr_list, Meth, overlap_list, cpg_ids) {
ol <- overlap_list[[as.character(dmr_idx)]]
if (!is.null(ol) && length(ol) > 0) {
cpg_subset <- cpg_ids[ol]
Bvalues_subset <- Meth[cpg_subset, , drop = FALSE]
mean_beta <- colMeans(Bvalues_subset, na.rm = TRUE)
} else {
mean_beta <- rep(NA, ncol(Meth))
}
dmr_name <- dmr_list[dmr_idx]
dmr_means <- as.data.frame(t(mean_beta))
rownames(dmr_means) <- dmr_name
colnames(dmr_means) <- colnames(Meth)
return(dmr_means)
}

Setup parallelism and run

message("Starting DMR Beta calculation...")
plan(multisession)
DMR_beta_list <- future_lapply(seq_along(dmr_list), calculate_dmr_means,
dmr_list = dmr_list,
Meth = Meth,
overlap_list = overlap_list,
cpg_ids = cpg_ids)

Combine and write result

DMR_beta <- do.call(rbind, DMR_beta_list)

[PeteHaitch]

Please provide a minimal reproducible example, i.e. some code someone can copy+paste to reproduce the issue.

[RicardoMartins-Ferreira]

Hi, thank you for the reply. I was able to overcome this by adapting the above script to account for the weight of the coverage in methylation calculation (script bellow). So, right now the mean meth values of all cpgs in a dmr coincides with the values in the t-statistics results.
This was more of a conceptual question, how are the mean methylation in the t-statists calculated? Does it consider coverage weight, CpG density, etc?

Filter CpGs based on coverage (>=10 in at least 3 samples)

bsseq.cov <- getCoverage(bsseq.fit)
keepLoci <- rowSums(bsseq.cov >= 10) >= 3
bsseq.fit <- bsseq.fit[keepLoci, ]

Set sample names

sampleNames(bsseq.fit) <- bsseq.fit@colData$Sample

Load and filter DMRs

dmr_results <- fread("/ijc/LABS/BALLESTAR/DATA/RICARDO/MeiBao_progesterone/WGBS/Analysis/All_DMRs_combined_qcutoff0.95_cpg5_beta20.txt")
dmr_results$chr <- gsub("\.[0-9]+$", "", dmr_results$chr)
dmr_results <- dmr_results[!dmr_results$comparison %in% c("MO_vs_iDC_pro", "MO_vs_iDC", "MO_vs_mDC_pro", "MO_vs_mDC"), ]

Add unique DMR names

dmr_results$Name <- paste0(dmr_results$chr, "", dmr_results$start, "", dmr_results$end)
dmr_list <- dmr_results$Name

Create CpG GRanges

ranges <- as.data.frame(bsseq.fit@rowRanges)
ranges$seqnames <- gsub("\.[0-9]+$", "", ranges$seqnames)
row.names(ranges) <- paste0(ranges$seqnames, "", ranges$start, "", ranges$end)
cpg_ids <- rownames(ranges)

cpg_gr <- makeGRangesFromDataFrame(ranges,
keep.extra.columns = FALSE,
ignore.strand = TRUE,
seqnames.field = "seqnames",
start.field = "start",
end.field = "end")

cpg_gr$ID <- cpg_ids

Create DMR GRanges

dmr_gr <- makeGRangesFromDataFrame(dmr_results,
seqnames.field = "chr",
start.field = "start",
end.field = "end",
keep.extra.columns = TRUE,
ignore.strand = TRUE)
dmr_gr$Name <- dmr_results$Name

Get overlaps between DMRs and CpGs

overlaps <- findOverlaps(dmr_gr, cpg_gr, ignore.strand = TRUE)
overlap_list <- split(subjectHits(overlaps), queryHits(overlaps))

Extract smoothed methylation and raw coverage matrices

Meth <- getMeth(bsseq.fit, type = "smooth") %>% as.data.frame()
Cov <- getCoverage(bsseq.fit, type = "Cov") %>% as.data.frame()

colnames(Meth) <- colnames(Cov) <- bsseq.fit@colData$Sample
rownames(Meth) <- rownames(Cov) <- rownames(ranges)

Function to calculate coverage-weighted mean methylation for each DMR

calculate_dmr_means <- function(dmr_idx, dmr_list, Meth, Cov, overlap_list, cpg_ids) {
ol <- overlap_list[[as.character(dmr_idx)]]
if (!is.null(ol) && length(ol) > 0) {
cpg_subset <- cpg_ids[ol]
meth_subset <- Meth[cpg_subset, , drop = FALSE]
cov_subset <- Cov[cpg_subset, , drop = FALSE]

# Coverage-weighted methylation
weighted_meth <- colSums(meth_subset * cov_subset, na.rm = TRUE) / colSums(cov_subset, na.rm = TRUE)

} else {
weighted_meth <- rep(NA, ncol(Meth))
}

dmr_name <- dmr_list[dmr_idx]
dmr_means <- as.data.frame(t(weighted_meth))
rownames(dmr_means) <- dmr_name
colnames(dmr_means) <- colnames(Meth)
return(dmr_means)
}

Run calculation in parallel

message("Starting DMR Beta calculation using coverage-weighted smoothed methylation values...")
plan(multisession, workers =2)

DMR_beta_list <- future_lapply(seq_along(dmr_list), calculate_dmr_means,
dmr_list = dmr_list,
Meth = Meth,
Cov = Cov,
overlap_list = overlap_list,
cpg_ids = cpg_ids)

Combine result into one data frame

DMR_beta <- do.call(rbind, DMR_beta_list)

[PeteHaitch]

This was more of a conceptual question, how are the mean methylation in the t-statists calculated? Does it consider coverage weight, CpG density, etc?

I'm afraid it's been several years since I thought about or used bsseq, so I don't recall off the top of my head.
But the documentation ?bsseq::BSmooth.tstat and reading the code (https://github.com/hansenlab/bsseq/blob/devel/R/BSmooth.tstat.R) are where I would start if I had these questions.
@kasperdanielhansen might be able to answer specific questions

[zetian-jia]

I’ve also encountered this issue before. The t-statistics in bsseq are calculated after BSmooth, and the package generally focuses on the smoothed methylation values—so the t-test is performed on the smoothed data (i.e., after the BSmooth fit). As a result, it’s not appropriate to evaluate methylation differences based on individual CpG sites directly. Instead, the results should be interpreted at the DMR level after merging.
From my own comparisons, I’ve noticed that the smoothing process tends to adjust methylation levels across neighboring CpGs, which can significantly affect the interpretation. I’ve compared results before and after smoothing and observed this correction effect clearly.
I hope this helps clarify the discrepancy. That said, in my current work, I’ve implemented some custom modifications to bsseq that allow me to access single-CpG statistics more directly. However, this is still a work in progress and not yet ready for broader use.