I want to identify differentially methylated probes(DMPs) between two groups.
After preprocessed raw data, 666805 probes remained. Differential methylation analysis was performed on the remain probes using limma. Age and gender information were used as covariates. Probes with BH adjusted P value<0.05 were considered as DMPs.
As a result, 134841 probes were identified as DMPs. It reveals nearly 1/5 probes was differentially methylated between our two groups. btw :the chip we used is Illumina Infinium MethylationEPIC BeadChip. It is kind of strange because samples from both groups were healthy people. They just grouped by different living environment.
So, we asked one technician to perform differential methylation analysis on the same 666805 probes. As a result, 8388 probes observed as DMPs,using IMA. I compared those two result, only 3037 probes were the same.
It is so strange that the numbers of two results is extremely different. I am newcomer to this research field, and i did not figure out the errors in my code. So i paste the code for DMP detection here and hope someone can help me.
library(limma)
targets<-read.table("Samplesheet_basicinfo",header=TRUE, stringsAsFactors=FALSE)
betaval<-read.delim("data_afterpreprocessed.txt",header=TRUE,stringsAsFactors=FALSE)
mval<-log2(betaval/(1-betaval))
var<-model.matrix(~Age + as.factor(Gender) + as.factor(Group),data=targets)
colnames(var)<-c("control","age","gender","group")
fit<-lmFit(mval,var)
fit2<-eBayes(fit)
probe<-topTable(fit2,adjust="BH",coef=4,num=Inf)
sig.probe<-probe[which(probe$adj.P.Val<=0.05),]
write.table(sig.probe,file="Result.DMP.txt",sep="\t",quote=TRUE)
Assuming the "data_afterpreprocessed.txt" file contains beta-values of normalised samples, in the same order as the samples in the samplesheet file, I can't see anything terribly wrong with the code. Perhaps if you do an
plotMD(fit2, column=4)
orvolcanoplot(fit2, coef=4)
you might see what's driving so much differential methylation. Or pick a probe that IMA fails to report as significant, and examine its mval profile and the coefficients and residuals from limma.Could you please help me how to make Methylation data input file like should I take CpG probes as rows or column? Many thanks
You should ask a new question, but the short answer is that features such as CpG probes are in rows.