差异甲基化分析eDMR软件分析说明

#############################################
 
library(edmr)
 
library(methylKit)
 
library(GenomicRanges)
 
library(mix tools)
 
#首先使用bismark，然后使用methylkit(具体可参见：https://www.plob.org/2014/10/22/8831.html)
 
file.list=list( "1h_CpG.txt","HL_CpG.txt")
 
myobj=read(file.list,sample.id=list("h0","L0"),assembly="hg19",treatment=c(0,1))
 
meth=unite(myobj, destrand=FALSE)
 
myDiff=calculateDiffMeth(meth,num.cores=8)
 
# fitting the bimodal normal distribution to CpGs distribution
 
tiff("h0_VS_L0.tiff")
 
myMixmdl=myDiff.to.mixmdl(myDiff, plot=T, main="h0_VS_L0")
 
dev.off()
 
# plot cost function and the determined distance cutoff
 
tiff("h0_VS_L0_cost_function.tiff")
 
plotCost(myMixmdl, main="cost function")
 
dev.off()
 
# calculate all DMRs candidate
 
mydmr=edmr(myDiff, mode=1, ACF=TRUE)
 
# further filtering the DMRs
 
mysigdmr=filter.dmr(mydmr)
 
## annotation加载注释文件这两个文件可以在edmr上下载
 
# get genebody annotation GRangesList object
 
genebody=genebody.anno(file="/DG/home/fyc/methylation/eDMR/hg19_refseq_all_types.bed")
 
# get CpG islands and shores annotation
 
cpgi=cpgi.anno(file="/DG/home/fyc/methylation/eDMR/hg19_cpgisland_all.bed")
 
# plot the eDMR genebody annotation
 
tiff("eDMR_genebody_annotation")
 
plot.dmr.distr(mysigdmr, genebody, main="eDMR genebody annotation", xlab="DMR count")
 
dev.off()
 
# plot the eDMR CpG islands and shores annotation
 
tiff("eDMR_CpG_annotation")
 
plot.dmr.distr(mysigdmr, cpgi, main="eDMR CpG islands and shores annotation", xlab="DMR count")
 
dev.off()
 
# prepare genes for pathway analysis with significant DMRs at its promoter regions
 
dmr.genes=get.dmr.genes(myDMR=mysigdmr, subject=genebody$promoter, id.type="gene.symbol")
 
dmr.genes
 
###########