Dear all,

currently Im working on methylation analysis of Illumina 450k beadchip data using the minfi workflow for identifying differentially methylated CpG sites. The design consist of 55 samples and we have different groups of covariates and one group for example has 3 levels. When applying dmp finder we obtain different significantly methylated CpG sites compared to limfit by limma. I read that the minfi dmpfinder is a wrapper around lmfit. But why are the results so different then? Appreciate any help and hints for a solution.

Phenotype data

Command dmp finder:

dmps.cat <- dmpFinder(m, pheno=phenotype_data$GroupSurgery, type="categorical")

Command limma

group1 <- factor(phenotype_data$GroupSurgery, levels=c("1","2","3"))

design <- model.matrix(~group1)
fit <- lmFit(m,design)
fit <- eBayes(fit)
topTable(fit, coef=group1)

Results dmp finder

Results lmfit

Two reasons. First, dmpFinder is smart enough to do an F-test that doesn't include the intercept (which estimates the mean of surgery group 1), whereas you did an F-test that does include this intercept. That doesn't make any sense; the F-test is testing the null hypothesis that all coefficients are equal to zero, and there is no expectation that the mean M-value for the surgery group 1 should be zero. The correct test would be

topTable(fit, 2:3)

Second, dmpFinder by default doesn't do the empirical Bayes step, whereas you are doing that when using limma directly. With 50 samples or so, I would be surprised to see any materially different results between the two, but there may be some small differences. But the main problem is that you are doing a nonsensical F-test when you use limma directly.