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)

Your limma code looks fine, assuming Gender and Group only have two levels and Age is a real-valued covariate. The only extra things would be to set trend=TRUE and robust=TRUE in eBayes, but I don't think it would make much of a dent in your 130000 DMPs. If you think that number's too high, investigate whether there are confounding factors (e.g., batch) that might be driving systematic differences between groups.

I can't speak for why the limma results disagree with those from IMA, you'll have to ask the IMA authors. Also see comments about testing for differential methylation: C: What is logFC and "B" in methylation analysis ??.

It is not entirely clear that there is any problem here.

IMA uses limma as the testing method in their vignette:

   http://rforge.net/IMA/meth450k.pdf

however the default in IMA is to do a Wilcoxon test for each site. Obviously a Wilcoxon test will be far less powerful than the limma test when you don't have a lot of replicates and will result in fewer DMPs. If you tell your "technician" to use the "limma" option when running IMA, you may find that it gives the same result as your limma analysis. If not, then the pre-processing must be different in some way.

Have you considered trying missMethyl (which also uses limma under the hood):

   https://bioconductor.org/packages/release/bioc/html/missMethyl.html