参考：

https://bioconductor.riken.jp/packages/3.9/bioc/vignettes/methylKit/inst/doc/methylKit.html

//www.greatytc.com/p/5c27908ff1e3

//www.greatytc.com/p/da74b4975019

//www.greatytc.com/p/88031796f333

rm(list=ls()); gc();

library(methylKit)

file.list = list("CYWD-1.CX_report.txt_CG_methykit.txt", "CYWD-2.CX_report.txt_CG_methykit.txt", "CYWD-3.CX_report.txt_CG_methykit.txt", "CY-1.CX_report.txt_CG_methykit.txt", "CY-2.CX_report.txt_CG_methykit.txt", "CY-3.CX_report.txt_CG_methykit.txt")

#myobjDB=methRead(file.list,sample.id=list("CYWD1","CYWD2","CYWD3","CY1","CY2","CY3"),assembly="Cs01",treatment=c(1,1,1,0,0,0),context="CpG",mincov = 10)

# 以数据库模式读入大文件，以减轻服务器内存压力, 二选一读入数据

myobjDB=methRead(file.list, sample.id=list("CYWD1","CYWD2","CYWD3","CY1","CY2","CY3"),assembly="Cs01", treatment=c(1,1,1,0,0,0), context="CpG", mincov = 10, dbtype = "tabix", dbdir = "methylDB",)

print(myobjDB[[1]]@dbpath)

## [1] "/tmp/RtmpBYQMT8/Rbuild4258771d5438/methylKit/vignettes/methylDB/test1.txt.bgz"

#去除覆盖率较低的reads(count<10)

myobjDB=filterByCoverage(myobjDB,lo.count=10,lo.perc=NULL,hi.count=NULL,hi.perc=99.9)

#归一化

myobjDB <- normalizeCoverage(myobjDB)

####还可以使用reorganize重新提取样本和处理信息，构建新的对象

####myobjDB2 = reorganize(myobjDB,sample.ids=c("test1","ctrl2"),treatment=c(1,0))

#合并所有样本中共有的位点

methDB = unite(myobjDB,destrand = F,suffix="CYWD")

####选取子集

####meth2 =reorganize(methDB,sample.ids=c("CYWD1","CY3"),treatment=c(1,0) ,suffix = "output_name")

head(methDB)

##save(methDB,file = "CYWD.Rdata")

##load(file = "CYWD.Rdata")

png("cor.png")

getCorrelation(methDB,method = "pearson", plot=TRUE)

dev.off()

png("test.png") ##样品聚类

clusterSamples(methDB,dist="correlation",method="ward.D",plot=TRUE)

dev.off()

png("PCA.png")

PCASamples(methDB)

dev.off()

#得到甲基化比率

#perc.meth=percMethylation(methDB)

#计算差异甲基化位点

myDiffDB = calculateDiffMeth(methDB,mc.cores=20)

myDiff25p.hyper=getMethylDiff(myDiffDB,difference=25,qvalue=0.01,type="hyper")

myDiff25p.hypo=getMethylDiff(myDiffDB,difference=25,qvalue=0.01,type="hypo")

myDiff25p=getMethylDiff(myDiffDB,difference=25,qvalue=0.01)

##diffMethPerChr(myDiffDB,plot=T,qvalue.cutoff=0.01, meth.cutoff=25)

write.table(myDiff25p.hyper,file = "WD25p01hyper.txt",sep = "\t", row.names = TRUE, col.names = TRUE)

write.table(myDiff25p.hypo,file = "WD25p01hypo.txt",sep = "\t", row.names = TRUE, col.names = TRUE)

write.table(myDiff25p,file = "WD25p01all.txt",sep = "\t", row.names = TRUE, col.names = TRUE)

#差异甲基化区域分析

myobj_lowCov = methRead(file.list, sample.id=list("CYWD1","CYWD2","CYWD3","CY1","CY2","CY3"),assembly="Cs01",treatment=c(1,1,1,0,0,0),context="CpG",mincov = 3) ###区域比较时，对单个碱基的覆盖度要求较低

tiles = tileMethylCounts(myobj_lowCov,win.size=1000,step.size=1000,cov.bases = 10,suffix = "CYWD")

head(tiles[[1]],3)

meth_tiles = unite(tiles)

head(meth_tiles)

myDiff_tiles = calculateDiffMeth(meth_tiles,mc.cores=20)

myDiff25pregions.hyper=getMethylDiff(myDiff_tiles,difference=25,qvalue=0.01,type="hyper")

myDiff25pregions.hypo=getMethylDiff(myDiff_tiles,difference=25,qvalue=0.01,type="hypo")

myDiff25pregions=getMethylDiff(myDiff_tiles,difference=25,qvalue=0.01)

write.table(myDiff25pregions.hyper,file = "WDregion25p01hyper.txt",sep = "\t", row.names = TRUE, col.names = TRUE)

write.table(myDiff25pregions.hypo,file = "WDregion25p01hypo.txt",sep = "\t", row.names = TRUE, col.names = TRUE)

write.table(myDiff25pregions,file = "WDregion25p01all.txt",sep = "\t", row.names = TRUE, col.names = TRUE)

#差异甲基化位点/区域注释

library(genomation)

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

#我们可以在USCS等网站下载所需物种的bed文件，随后使用genomation将其转换为可注释的对象，也可以使用TxDb系列包自己转换，需要注意，进行注释的GRangesList对象必须包括启动子/内含子/外显子/基因间区域上的信息。

#作者：nnlrl

#链接：//www.greatytc.com/p/da74b4975019

#来源：简书

#著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

#读入个人数据,读取基因注释文件(bed格式),需要提前准备bed12文件(参考文献为https://www.xknote.com/ask/60f2f34e89da1.html)，使用了gff3ToGenePred，然后使用了UCSC实用程序中的genePredToBed工具。这将输出一个12列的.bed；

#文件转换用软件安装

#mkdir UCSC_tools

#cd UCSC_tools

#rsync -avzPL rsync://hgdownload.soe.ucsc.edu/genome/admin/exe/linux.x86_64/gff3ToGenePred .

#rsync -avzPL rsync://hgdownload.soe.ucsc.edu/genome/admin/exe/linux.x86_64/genePredToBed .

#rsync -avzPL rsync://hgdownload.soe.ucsc.edu/genome/admin/exe/linux.x86_64/* .##幕布所有工具

#gff3ToGenePred CsTGY.gff CsTGY.pred

#genePredToBed CsTGY.pred CsTGY.bed

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

#参数up.flank和down.flank用于定义启动子区域，默认TSS上下游各1000bp为启动子，这里定义TSS上游2000bp为启动子区。

gene.obj = readTranscriptFeatures("/public/home/tanger/Kongweilong/test/extdata/CsTGY.bed",up.flank=2000,down.flank=0)

#genomation的函数是基于GRanges对象，首先将myDiff25pregions转变为GRanges格式

annotateWithGeneParts(as(myDiff25pregions,"GRanges"),gene.obj)

#4.1 区域分析

#用于汇总启动子等区域的DNA甲基化信息

promoters = regionCounts(myobj,gene.obj$promoters) #####好像不识别压缩形式的输入myobj

promoters = regionCounts(myobj_lowCov,gene.obj$promoters)

promoters_meth_tiles = regionCounts(meth_tiles,gene.obj$promoters)

head(promoters[[1]])

#4.2 注释对象的相关操作函数

#当得到差异甲基化区域的注册信息后，可以通过getAssociationWithTSS函数获得其和TSS之间的距离，以及最近的基因。

diffAnn=annotateWithGeneParts(as(myDiff25pregions,"GRanges"),gene.obj)

head(getAssociationWithTSS(diffAnn))

#还可以得到与内含子/外显子/启动子重叠的差异甲基化区域的百分比/数量

getTargetAnnotationStats(diffAnn,percentage=TRUE,precedence=TRUE)