参考：
https://bioconductor.riken.jp/packages/3.9/bioc/vignettes/methylKit/inst/doc/methylKit.html
https://www.jianshu.com/p/5c27908ff1e3
https://www.jianshu.com/p/da74b4975019
https://www.jianshu.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
#链接：https://www.jianshu.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)
网友评论