前言
今天小编为大家介绍一款差异甲基化分析软件包，methylKit用于分析和注释通过高通量亚硫酸氢盐测序获得的DNA甲基化信息。如果提供了适当的输入格式，它可以处理全基因组重亚硫酸盐测序数据。
1、数据输入
亚硫酸氢盐测序实验分为测试和对照样品。 测试样品可以来自疾病组织，而对照样品可以来自健康组织。首先读取数据，该数据集存储每个样本的甲基化信息。数据格式如下：
在以下示例中，前两个文件具有样本id“test1”和“test2”，并且由控制向量确定它们属于同一组。第三和第四个文件具有样本ID“ctrl1”和“ctrl2”，并且它们属于同一组。
##安装加载R包
source("http://bioconductor.org/biocLite.R")
biocLite("methylKit")
library(methylKit)
file.list <- list(system.file("extdata", "test1.myCpG.txt",package= "methylKit"), system.file("extdata",
"test2.myCpG.txt",package=
"methylKit"), system.file("extdata","control1.myCpG.txt",package=
"methylKit"), system.file("extdata","control2.myCpG.txt",package=
"methylKit"))
##将文件读取到一个methylRawList对象：myobj，treatment确定对照与测试组，assembly选择参考基因组，context参数控制甲基化位点类型。
myobj <- read(file.list, sample.id= list("test1","test2", "ctrl1", "ctrl2"),assembly = "hg18", treatment=
c(1,1, 0, 0), context = "CpG")
2、从已排序的Bismark比对中读取甲基化calls
read.bismark函数可以将Bismark 的SAM文件读入methylRaw或methylRawList对象。 SAM文件必须按染色体排序并读取位置列。
my.methRaw <- read.bismark(location= system.file("extdata","test.fastq_bismark.sorted.min.sam", package= "methylKit"),sample.id ="test1", assembly = "hg18", read.context= "CpG",save.folder= getwd())
3、样本的描述性统计
getMethylationStats可检查甲基化数据的基本统计数据，例如覆盖率和甲基化百分比。以下命令输出第二个样本的百分比甲基化统计信息：“test2
getMethylationStats(myobj[[2]], plot= F, both.strands = F)
getMethylationStats(myobj[[2]], plot= T, both.strands = F)
library("graphics")
getCoverageStats(myobj[[2]],plot = T, both.strands = F)
4、根据read覆盖度过滤样本
如果样品存在PCR偏差，那么需要丢弃具有非常高read覆盖率的碱基。 此外，我们还可以抛弃低read覆盖率的碱基，下面的代码过滤甲基化数据并丢弃覆盖率低于10X的碱基，并且丢弃每个样本中覆盖率超过99.9％的碱基。
filtered.myobj <- filterByCoverage(myobj,lo.count = 10,lo.perc = NULL,hi.count = NULL, hi.perc = 99.9)
5、整合样本数据
为了进一步分析，将所有样本合并到一个对象中，以便覆盖所有样本中的碱基对位置。 设置destrand = TRUE（默认值为FALSE）将合并CpG二核苷酸的两条链上的read。这提供了更好的覆盖范围，但仅在观察CpG甲基化时（对于CpH甲基化，这会在随后的分析中导致错误结果）建议。unite（）函数将返回一个methylBase对象，包含所有样本中覆盖的区域/碱基的甲基化信息。
meth <- unite(myobj, destrand= FALSE)
我们可以使用getCorrelation来检查样本之间的相关性。这个函数可以绘制散点图和相关系数，或者只是一个相关矩阵.
getCorrelation(meth, plot= T)
6、 寻找差异甲基化的碱基或区域
calculateDiffMeth（）函数是计算差异甲基化的主要函数。根据每组样本大小，它将使用Fisher精确或逻辑回归来计算P值。 使用SLIM方法将P值调整为Q值.
myDiff <- calculateDiffMeth(meth)
在q值计算之后，我们可以根据q值和甲基化差异截止百分比选择差异甲基化区域/碱基。以下位选择q值<0.01且甲基化差异百分比大于25％的碱基。如果指定type =“hyper”或type =“hypo”选项，则会出现超甲基化或低甲基化区域/碱基。
myDiff25p.hyper <- get.methylDiff(myDiff,difference = 25,qvalue = 0.01,type = "hyper")
# get hypo methylated bases
myDiff25p.hypo <- get.methylDiff(myDiff,difference = 25,qvalue = 0.01,type = "hypo")
# #
# get all differentially methylated bases
myDiff25p <- get.methylDiff(myDiff,difference = 25,qvalue = 0.01)
7、注释差异甲基化的碱基或区域
我们可以基于基因注释来注释我们的差异甲基化区域/碱基。 在这个例子中，从bed文件中读取基因注释，并用这些信息注释差异甲基化区域。在启动子/内含子/外显子/基因间区域中差异甲基化区域的百分比。 使用Bioconductor.org提供的GenomicFeatures软件包可以获取类似的基因注释
gene.obj<- read.transcript.features(system.file("extdata","refseq.hg18.bed.txt", package= "methylKit"))
annotate.WithGenicParts(myDiff25p,gene.obj)
同样，我们可以读取CpG岛注释并用他们注释我们差异甲基化的碱基/区域。
#read the shores and flanking regions and name the flanks
# as shores and CpG islands as CpGi
cpg.obj <- read.feature.flank(system.file("extdata","cpgi.hg18.bed.txt", package= "methylKit"),
feature.flank.name = c("CpGi","shores"))
#
diffCpGann <- annotate.WithFeature.Flank(myDiff25p,cpg.obj$CpGi,cpg.obj$shores, feature.name = "CpGi",flank.name= "shores")
参考文献：
Akalin, A., et al., methylKit: a comprehensive R package for theanalysis of genome-wide DNA methylation profiles. Genome Biol, 2012. 13(10): p. R87.
