inconsistent calculateDiffMeth results when re-ordering methylBase columns and swapping treatment 0/1 annotations

[hd2326]

Greetings!

Differences between code chunks were highlighted.

Case 1, new methylBase object, re-order columns:

dat1 <- getData(methylBase.obj)
dat1 <- dat1[, c(1:4, methylBase.obj@coverage.index, methylBase.obj@numCs.index, methylBase.obj@numTs.index)]
dat1 <- new("methylBase", as.data.frame(dat1),
sample.ids=c(paste("C", 1:2, sep=""), paste("N", 1:2, sep="")),
treatment=c(rep(1, 2), rep(0, 2)),
coverage.index=5:8,
numCs.index=9:12,
numTs.index=13:16,
assembly="simulation", context="simulation", destranded=F, resolution="base")
dat1 <- calculateDiffMeth(dat1, adjust="fdr")

Case 2, new methylBase object, without re-ordering:

dat2 <- getData(methylBase.obj)
dat2 <- new("methylBase", as.data.frame(dat2),
sample.ids=c(paste("C", 1:2, sep=""), paste("N", 1:2, sep="")),
treatment=c(rep(1, 2), rep(0, 2)),
coverage.index=methylBase.obj@coverage.index,
numCs.index=methylBase.obj@numCs.index,
numTs.index=methylBase.obj@numTs.index,
assembly="simulation", context="simulation", destranded=F, resolution="base")
dat2 <- calculateDiffMeth(dat2, adjust="fdr")

Case 3, new methylBase object, swap "treatment" 0/1 annotations:

dat3 <- getData(methylBase.obj)
dat3 <- new("methylBase", as.data.frame(dat3),
sample.ids=c(paste("C", 1:2, sep=""), paste("N", 1:2, sep="")),
treatment=c(rep(0, 2), rep(1, 2)),
coverage.index=methylBase.obj@coverage.index,
numCs.index=methylBase.obj@numCs.index,
numTs.index=methylBase.obj@numTs.index,
assembly="simulation", context="simulation", destranded=F, resolution="base")
dat3 <- calculateDiffMeth(dat3, adjust="fdr")

Case 4, control:

dat4 <- calculateDiffMeth(methylBase.obj, adjust="fdr")

Results:

Case 1 != Case 4

sum(dat1$pvalue != dat4$pvalue)
[1] 963
sum(dat1$qvalue != dat4$qvalue)
[1] 963

Case 2 == Case 4

sum(dat2$pvalue != dat4$pvalue)
[1] 0
sum(dat2$qvalue != dat4$qvalue)
[1] 0

Case 3 != Case 4

sum(dat3$pvalue != dat4$pvalue)
[1] 505
sum(dat3$qvalue != dat4$qvalue)
[1] 462

Thoughts?

[alexg9010]

Hi @hd2326, Interesting tests. Thanks for the details! First of all, I would like to mention that we have a function `reorganize()` that allows changing the order and selection of samples. But given your test cases I can confirm this is a real bug in `.calculateDiffMeth()`. It never reads ***@***.******@***.******@***.***`, but hardcodes column positions instead. This assumes the canonical interleaved layout (`coverage,numCs,numTs` per sample, in that order). So: - **Case 1** fails because you regrouped columns — the index slots were set correctly, but the function ignores them and pulls the wrong columns. Real bug, no warning given. - **Case 2** passes only because the layout happened to stay canonical. - **Case 3** (treatment swap) — mathematically the GLM deviance/p-value should be invariant to recoding 0/1 labels. I'd guess your diffs there are floating-point noise from `glm.fit`'s IRLS convergence rather than a logic bug — worth checking magnitude with `all.equal()` instead of exact `!=` before treating it as separate from Case 1. Fix should either use the stored index slots directly or normalize/reorder columns before computing `Ccols`/`Tcols`. I'd be happy to review a pull request that adds your checks as a unit test and fixes the issue. Best, Alex hd2326 ***@***.***> schrieb am Sa., 29. Aug. 2026, 00:52:

…

*hd2326* created an issue (al2na/methylKit#370) <#370> Greetings! Differences between code chunks were highlighted. *Case 1, new methylBase object, re-order columns:* dat1 <- getData(methylBase.obj) dat1 <- dat1[, c(1:4, ***@***.***, ***@***.***, ***@***.***)] dat1 <- new("methylBase", as.data.frame(dat1), sample.ids=c(paste("C", 1:2, sep=""), paste("N", 1:2, sep="")), treatment=c(rep(1, 2), rep(0, 2)), *coverage.index=5:8,* *numCs.index=9:12,* *numTs.index=13:16,* assembly="simulation", context="simulation", destranded=F, resolution="base") dat1 <- calculateDiffMeth(dat1, adjust="fdr") *Case 2, new methylBase object, without re-ordering:* dat2 <- getData(methylBase.obj) dat2 <- new("methylBase", as.data.frame(dat2), sample.ids=c(paste("C", 1:2, sep=""), paste("N", 1:2, sep="")), treatment=c(rep(1, 2), rep(0, 2)), ***@***.***, ***@***.***, ***@***.***, assembly="simulation", context="simulation", destranded=F, resolution="base") dat2 <- calculateDiffMeth(dat2, adjust="fdr") *Case 3, new methylBase object, swap "treatment" 0/1 annotations:* dat3 <- getData(methylBase.obj) dat3 <- new("methylBase", as.data.frame(dat3), sample.ids=c(paste("C", 1:2, sep=""), paste("N", 1:2, sep="")), *treatment=c(rep(0, 2), rep(1, 2)),* ***@***.***, ***@***.***, ***@***.***, assembly="simulation", context="simulation", destranded=F, resolution="base") dat3 <- calculateDiffMeth(dat3, adjust="fdr") *Case 4, control:* dat4 <- calculateDiffMeth(methylBase.obj, adjust="fdr") *Results:* Case 1 != Case 4 sum(dat1$pvalue != dat4$pvalue) [1] 963 sum(dat1$qvalue != dat4$qvalue) [1] 963 Case 2 == Case 4 sum(dat2$pvalue != dat4$pvalue) [1] 0 sum(dat2$qvalue != dat4$qvalue) [1] 0 Case 3 != Case 4 sum(dat3$pvalue != dat4$pvalue) [1] 505 sum(dat3$qvalue != dat4$qvalue) [1] 462 Thoughts? — Reply to this email directly, view it on GitHub <#370?email_source=notifications&email_token=ADK7JD6U7KT5ZBJ7IN24IX35MIEJJA5CNFSL4Z3JMQ5C6L3HNF2C22DVMIXUS43TOVSS6NJSHAZDSMZWGA3TFJTSMVQXG33OVJZXKYTTMNZGSYTFMSSWK5TFNZ2KYZTPN52GK4S7MNWGSY3L>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADK7JDYGH5IA4P24NJVRTH35MIEJJAVCNFSNUABEKJSXA33TNF2G64TZHMYTKNBSGQYTIMJ3JFZXG5LFHM2TEOBSHEZTMMBXGKQXMAQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[hd2326]

Hey Alex thank you for the quick response!

Yes I confirm all.equal solves the treatment swap problem.

all.equal(dat2$pvalue, dat3$pvalue)
[1] TRUE
all.equal(dat2$qvalue, dat3$qvalue)
[1] TRUE