methylKit是一个用于分析DNA甲基化测序数据的R包，其核心功能是差异甲基化分析和差异甲基化位点/区域的注释。
1. 加载原始数据
DNA甲基化文件格式如下:
## chrBase chr base strand coverage freqC freqT
## 1 chr21.9764539 chr21 9764539 R 12 25.00 75.00
## 2 chr21.9764513 chr21 9764513 R 12 0.00 100.00
## 3 chr21.9820622 chr21 9820622 F 13 0.00 100.00
## 4 chr21.9837545 chr21 9837545 F 11 0.00 100.00
## 5 chr21.9849022 chr21 9849022 F 124 72.58 27.42
每一行是一个甲基化位点，其中，coverage代表覆盖这个位点的reads数，freqC和freqT分别代表发生甲基化和未发生甲基化胞嘧啶的比例。这种纯文本格式内容直观，文件小，读取速度快。
此文件可以通过perl脚本转化bismark输出文件CX_report.txt得到。
perl BismarkCX2methykit.pl <CX_report.txt>
#BismarkCX2methykit.pl代码见末尾。
文件加载过程如下
library(methylKit)
file.list = list("/public/home/FengSL/anaconda2/lib/R/library/methylKit/extdata/test1.myCpG.txt",\
"/public/home/FengSL/anaconda2/lib/R/library/methylKit/extdata/test2.myCpG.txt",\
"/public/home/FengSL/anaconda2/lib/R/library/methylKit/extdata/control1.myCpG.txt",\
"/public/home/FengSL/anaconda2/lib/R/library/methylKit/extdata/control2.myCpG.txt")
# read the files to a methylRawList object: myobj
myobj=methRead(file.list,
sample.id=list("test1","test2","ctrl1","ctrl2"),
assembly="hg18",
treatment=c(1,1,0,0),
context="CpG",
mincov = 10
)
# 以数据库模式读入大文件，以减轻服务器内存压力
myobjDB=methRead(file.list,
sample.id=list("test1","test2","ctrl1","ctrl2"),
assembly="hg18",
treatment=c(1,1,0,0),
context="CpG",
mincov = 10,
dbtype = "tabix",
dbdir = "methylDB",
)
print(myobjDB[[1]]@dbpath)
## [1] "/tmp/RtmpBYQMT8/Rbuild4258771d5438/methylKit/vignettes/methylDB/test1.txt.bgz"
#还可以使用reorganize重新提取样本和处理信息，构建新的对象
myobj2=reorganize(myobj,sample.ids=c("test1","ctrl2"),treatment=c(1,0) )
myobjDB2 = reorganize(myobjDB,sample.ids=c("test1","ctrl2"),treatment=c(1,0))
2. 合并所有的样本
meth = unite(myobj,destrand = F,suffix="output_name")
methDB = unite(myobjDB,destrand = F,suffix="output_name")
##选取子集
meth2 =reorganize(meth,sample.ids=c("test1","ctrl2"),treatment=c(1,0) ,suffix = "output_name")
head(meth)
## chr start end strand coverage1 numCs1 numTs1 coverage2 numCs2 numTs2
## 1 chr21 9853296 9853296 + 17 10 7 333 268 65
## 2 chr21 9853326 9853326 + 17 12 5 329 249 79
## 3 chr21 9860126 9860126 + 39 38 1 83 78 5
## 4 chr21 9906604 9906604 + 68 42 26 111 97 14
## 5 chr21 9906616 9906616 + 68 52 16 111 104 7
## 6 chr21 9906619 9906619 + 68 59 9 111 109 2
## coverage3 numCs3 numTs3 coverage4 numCs4 numTs4
## 1 18 16 2 395 341 54
## 2 16 14 2 379 284 95
## 3 83 83 0 41 40 1
## 4 23 18 5 37 33 4
## 5 23 14 9 37 27 10
## 6 22 18 4 37 29 8
在合并的过程中，默认情况下，只有所有的样本都包含该位点时，才会保留，本质就是取的所有样本的交集，如果你想要取并集，可以修改min.per.group参数的值，该参数的值代表每组中至少有多少个样本覆盖该位点时才保留，如果设置为1，就是取并集。
meth.min=unite(myobj,min.per.group=1L)
2.1样品相关性分析
getCorrelation(meth2,method = "pearson", plot=TRUE, nrow = 2e+06)
# method, default: "pearson",other options are "kendall" and "spearman"
#nrow, a numeric giving the number of lines to read in of methylBaseDB object, defalts to 2e6.
2.2样品聚类
clusterSamples(meth2,dist="correlation",method="ward",plot=TRUE)
#dist, the distance measure to be used. This must be one of "correlation"(default), "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski".
#method, the agglomeration method to be used. This should be (an unambiguous abbreviation of) one of "ward.D"(default), "ward.D2", "single", "complete", "average", "mcquitty", "median" or "centroid".
PCASamples(meth2)
3 差异分析
3.1 单碱基差异分析
在methylKit的差异分析针对的是合并后的甲基化表达谱，上面的合并文件中每一行是一个甲基化位点，那么差异分析的结果就是差异甲基化位点，使用函数calculateDiffMeth执行差异分析。
myDiff=calculateDiffMeth(meth,mc.cores=2)
myDiffDB = calculateDiffMeth(methDB,mc.cores=2)
#calculateDiffMeth支持多线程，这里使用mc.cores参数调用2个核心参与运算
# 提取超甲基化位点(hyper)
myDiff25p.hyper=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hyper")
#
# 提取去甲基化位点(hypo)
myDiff25p.hypo=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hypo")
#
#
# 提取全部差异甲基化位点(all)
myDiff25p=getMethylDiff(myDiff,difference=25,qvalue=0.01)
3.2 区域差异分析
执行区域分析之前，首先要合并甲基化区域文件，区域比较时，对单个碱基的覆盖度要求较低。
设定窗口和步长，比较区域内甲基化差异
##区域比较时，对单个碱基的覆盖度要求较低
myobj_lowCov = methRead(file.list,
sample.id=list("test1","test2","ctrl1","ctrl2"),
assembly="hg18",
treatment=c(1,1,0,0),
context="CpG",
mincov = 3
)
##窗口大小和步长均可调节 软件默认为1000
tiles = tileMethylCounts(myobj_lowCov,win.size=1000,step.size=1000,cov.bases = 10,suffix = "output_name")
head(tiles[[1]],3)
meth_tiles = unite(tiles)
head(meth_tiles)
#methylBase object with 3 rows
#--------------
# chr start end strand coverage1 numCs1 numTs1 coverage2 numCs2
#1 chr21 9906001 9907000 * 970 604 366 1878 1421
#2 chr21 10011001 10012000 * 2461 2169 292 493 444
#3 chr21 10012001 10013000 * 4932 4175 756 736 622
#numTs2 coverage3 numCs3 numTs3 coverage4 numCs4 numTs4
#1 457 1011 621 390 724 425 299
#2 49 704 421 283 252 178 74
#3 114 1737 1185 552 417 285 132
#--------------
#sample.ids: test1 test2 ctrl1 ctrl2
#destranded FALSE
#assembly: hg18
#context: CpG
#treament: 1 1 0 0
#resolution: region
合并完成以后，可以执行区域差异分析。
读取文件
meth_tiles = readMethylDB(mydbpath)
重新选择样品
myobjDB2 = reorganize(myobjDB,sample.ids=c("test1","ctrl2"),treatment=c(1,0),suffix = "output_name")
myDiff_tiles = calculateDiffMeth(meth_tiles,mc.cores=2)
提取差异区域
myDiff25p.hyper=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hyper")
myDiff25p.hypo=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hypo")
myDiff_tiles25p=getMethylDiff(myDiff_tiles,difference=25,qvalue=0.01)
4 差异甲基化位点/区域注释
注释首先要使用R包genomation获取基因/区域信息
library(genomation)
#读取基因注释文件(bed格式)
bed.file = system.file("extdata","refseq.hg18.bed.txt",package="methylKit")
bed.file
[1] "/public/home/FengSL/anaconda2/lib/R/library/methylKit/extdata/refseq.hg18.bed.txt"
gene.obj = readTranscriptFeatures(bed.file,up.flank=2000,down.flank=0)
参数up.flank和down.flank用于定义启动子区域，默认TSS上下游各1000bp为启动子，这里定义TSS上游2000bp为启动子区。
#genomation的函数是基于GRanges对象，首先将myDiff25p转变为GRanges格式
annotateWithGeneParts(as(myDiff25p,"GRanges"),gene.obj)
##这里bed文件同样可以是 转座子或CG岛的注释
4.1 区域分析
用于汇总启动子等区域的DNA甲基化信息
promoters = regionCounts(myobj,gene.obj$promoters)
promoters_meth = regionCounts(meth,gene.obj$promoters)
head(promoters[[1]])
## chr start end strand coverage numCs numTs
## 1 chr21 10011791 10013791 - 7953 6662 1290
## 2 chr21 10119796 10121796 - 1725 1171 554
## 3 chr21 10119808 10121808 - 1725 1171 554
## 4 chr21 13903368 13905368 + 10 10 0
## 5 chr21 14273636 14275636 - 282 220 62
## 6 chr21 14509336 14511336 + 1058 55 1003
4.2 注释对象的相关操作函数
当得到差异甲基化区域的注册信息后，可以通过getAssociationWithTSS函数获得其和TSS之间的距离，以及最近的基因。
diffAnn=annotateWithGeneParts(as(myDiff25p,"GRanges"),gene.obj)
# target.row is the row number in myDiff25p
head(getAssociationWithTSS(diffAnn))
# target.row dist.to.feature feature.name feature.strand
#60 1 792 NM_199260 -
#170 2 86634 NM_174981 +
#170.1 3 94634 NM_174981 +
#188 4 0 NM_024944 +
还可以得到与内含子/外显子/启动子重叠的差异甲基化区域的百分比/数量
getTargetAnnotationStats(diffAnn,percentage=TRUE,precedence=TRUE)
附
# BismarkCX2methykit.pl
#NC_000932.1 3 - 0 0 CHH CAT
#NC_000932.1 4 - 0 0 CHH CCA
#NC_000932.1 5 - 0 0 CHH CCC
#NC_000932.1 6 + 0 0 CG CGA
# chrbase chr base strnd coverage freqC freqT
open (IN,"$ARGV[0]");
open (OUT1,">${ARGV[0]}_CG_methykit.txt");
open (OUT2,">${ARGV[0]}_CHG_methykit.txt");
open (OUT3,">${ARGV[0]}_CHH_methykit.txt");
print OUT1 "chrbase\tchr\tbase\tstrand\tcoverage\tfreqC\tfreqT\n";
print OUT2 "chrBase\tchr\tbase\tstrand\tcoverage\tfreqC\tfreqT\n";
print OUT3 "chrBase\tchr\tbase\tstrand\tcoverage\tfreqC\tfreqT\n";
while (<IN>) {
chomp;
@a=split;
$chrbase=$a[0].".".$a[1];
$chr=$a[0];
$base=$a[1];
if ($a[2]=~/-/) { $strand= "R";
}
else { $strand ="F";}
$coverage=$a[3]+$a[4];
if ($coverage < 5) {next;}
$freqC= $a[3]/($a[3]+$a[4]) * 100;
$freqT= $a[4]/($a[3]+$a[4]) * 100;
$str=$chrbase."\t".$chr."\t".$base."\t".$strand."\t".$coverage."\t".$freqC."\t".$freqT."\n";
if ($a[5]=~ /CG/) { print OUT1 $str;
}
elsif($a[5]=~/CHG/) {print OUT2 $str;}
elsif($a[5]=~/CHH/){ print OUT3 $str;}
else{print $_."\n";}
}
