Dear all,

I have 450k methylation data from 120 samples. I perform differential methylation analysis (using the m-values) for a interval scaled variable (current_phenotype with levels of 1,2,3) with limma which identifies 1336 hypo- and 635 hypermethylated CpG´s.The significant CpG´s (FDR<0.05 & logFC > 2) are then used for an analysis of differentially methylated pathways using the missMethyl package (gometh function). Although I get several differential methylated pathways I have troubles with interpreting the results. Specifically, how do I know whether higher values of "current_phenotpye" are associated with increased or decreased pathway methylation? Moreover, does it make sense to subject both significant hyper- and hypermethylated CpGs to the pathway analysis (since the gometh function does only know which CpG is significant, however not the direction i.e. whether it is hypo- or hypermethylated)?

Thanks for your help,

Philipp 

library(limma)
library(missMethyl)

design2 = model.matrix(~current_phenotype)
fit2 = lmFit(m_values, design2)  
keep <- fit2$Amean > median(fit2$Amean)
fitEb <- eBayes(fit2[keep,], robust=T, trend=T)
summary(decideTests(fitEb))

   (Intercept) current_phenotype
-1         812              1336
0        20713            204111
1       184557               635

tt <- topTable(fitEb,coef=2,sort.by="p", p.value=0.05, lfc=2, adjust.method="BH",number=Inf)
     
gst.KEGG <- gometh(sig.cpg=rownames(tt), all.cpg=rownames(m_values), collection="KEGG", prior.prob = T)
gst.KEGG <- gst.KEGG[order(gst.KEGG$FDR),]
gst.KEGG  <- gst.KEGG[gst.KEGG$FDR<0.05,]
head(gst.KEGG)