引言
methylKit是一个R软件包，用于分析和注释通过高通量亚硫酸氢盐测序获得的DNA甲基化信息。该软件包旨在处理RRBS及其变体的测序数据。但是，如果提供正确的输入格式，它可能会处理全基因组的亚硫酸氢盐测序数据。
脊椎动物的DNA甲基化通常发生在CpG二核苷酸处，但是某些组织中（如胚胎干细胞）非CPG的位点也会被甲基化。DNA甲基化可以作为基因调控的表观遗传控制机制。甲基化可能会阻碍转录因子的结合和/或甲基化的碱基可能会被甲基结合域蛋白结合，而后者会募集染色质重塑因子。在这两种情况下，受调节基因的转录都会发生。另外，异常的DNA甲基化模式已经与许多人类恶性肿瘤相关联并且可以以可预测的方式使用。在恶性组织中，与正常组织相比，DNA被低甲基化或高甲基化。高和低甲基化位点的位置为许多疾病提供了独特的特征。
该软件主要使用
methylKit对象存储甲基化数据，这些对象在分析过程中可以与GRanges对象无缝衔接，作者提供了函数as(object,"GRanges")可以将其转化为GRanges对象进行其他分析。
读取甲基化调用文件
我们首先从具有methRead功能的亚硫酸氢盐测序中读取甲基化数据。以这种方式读取数据将返回一个methylRawList对象，该对象为每个覆盖的碱基存储每个样品的甲基化信息。甲基化调用文件基本上是文本文件，每个碱基包含甲基化分数。典型的甲基化数据文件格式如下，用户可以使用bismark的coverage文件进行简单的处理：
## 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
在大多数时候，亚硫酸氢盐测序实验都有测试样品和对照样品。测试样品可以来自疾病组织，而对照样品可以来自健康组织。用户可以读取一组甲基化调用文件，这些文件的测试/控制条件提供了treatment选项。为了进行后续分析，file.list，sample.id和处理选项应具有相同的顺序。在以下示例中，前两个文件的样本ID为“ test1”和“ test2”，并且根据处理向量确定，它们属于同一组。第三和第四文件具有样本ID“ ctrl1”和“ ctrl2”，并且它们与治疗向量所指示的属于同一组。
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") )
# 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" ，#CG岛类型，可选CPG，CHG等
#dbtype = "tabix", 指定tabix索引，在进行大型数据分析时避免内存不足
#dbdir = "methylDB"
)
同时软件支持从bam, sam文件中获取甲基化信息,比较耗时
my.methRaw=processBismarkAln( location =
system.file("extdata",
"test.fastq_bismark.sorted.min.sam",
package = "methylKit"),
sample.id="test1", assembly="hg18",
read.context="CpG", save.folder=getwd())
样本的描述统计
在进行差异分析之前，我们首先应该观察样本的分布，覆盖范围等信息，保证数据准确性，以下命令绘制了甲基化百分比分布的直方图。条形图上的数字表示该文件中包含的位置百分比。通常，甲基化百分比直方图的两端应有两个峰。大多数文件应该会产生相似的模式，在该图中我们看到很多甲基化程度较高的位置和很多甲基化程度较低的位置。
getMethylationStats(myobj[[2]],plot=TRUE,both.strands=FALSE)
#我们还可以以类似的方式绘制每个基本信息的读取覆盖率，条形图上的数字再次
#表示该文件中包含多少百分比的位置。遭受PCR复制偏差严重影响的实验将在直#方图的右侧有一个次要峰。
getCoverageStats(myobj[[2]],plot=TRUE,both.strands=FALSE)
根据覆盖率过滤样本可能很有用。尤其是，如果我们的样品遭受PCR复制偏差，则丢弃具有非常高读取覆盖率的碱基将很有用。此外，我们还希望丢弃读取覆盖率较低的碱基，足够高的读取覆盖率将提高统计检验的能力。下面的代码过滤a methylRawList并丢弃覆盖率低于10%的碱基，并且还丢弃每个样本中覆盖率超过99.9％的碱基。
filtered.myobj=filterByCoverage(myobj,lo.count=10,lo.perc=NULL,
hi.count=NULL,hi.perc=99.9)
为了进行进一步的分析，我们需要获取所有样本中涵盖的位点。以下函数会将所有样本合并到一个对象，以获取所有样本中涵盖的碱基对位置。设置destrand=TRUE（默认值为FALSE）将合并正负链。这样可以提供更好的覆盖范围，但是仅在查看CpG甲基化时才建议使用（对于CpH甲基化，这将在后续分析中导致错误的结果）。此外，该设置destrand=TRUE仅在以碱基对分辨率运行时才起作用，否则将此选项设置为TRUE无效。该unite()函数将返回一个methylBase对象，这将是我们进行所有比较分析的主要对象。该methylBase对象包含所有样本中涵盖的区域/碱基的甲基化信息。
meth=unite(myobj, destrand=FALSE)
head(meth)
## chr start end strand coverage1 numCs1 numTs1 coverage2 numCs2
## 1 chr21 9853296 9853296 + 17 10 7 333 268
## 2 chr21 9853326 9853326 + 17 12 5 329 249
## 3 chr21 9860126 9860126 + 39 38 1 83 78
## 4 chr21 9906604 9906604 + 68 42 26 111 97
## 5 chr21 9906616 9906616 + 68 52 16 111 104
## 6 chr21 9906619 9906619 + 68 59 9 111 109
## numTs2 coverage3 numCs3 numTs3 coverage4 numCs4 numTs4
## 1 65 18 16 2 395 341 54
## 2 79 16 14 2 379 284 95
## 3 5 83 83 0 41 40 1
## 4 14 23 18 5 37 33 4
## 5 7 23 14 9 37 27 10
## 6 2 22 18 4 37 29 8
默认情况下，unite函数会产生所有样本中覆盖的碱基/区域。
#查看样本相关性
getCorrelation(meth,plot=TRUE)
#样本间聚类
clusterSamples(meth, dist="correlation", method="ward", plot=TRUE)
#PCA
PCASamples(meth)
在某些情况下，可能希望在平铺窗口上汇总甲基化信息，而不是进行碱基对分辨率分析。methylKit提供进行此类分析的功能。下面的功能将窗口长度为1000bp，步长为1000bp的基因组平铺，并总结了这些平铺中的甲基化信息。在这种情况下，它返回一个methylRawList对象，该对象可以被送入unite并使用calculateDiffMeth以获取差异甲基化区域。该功能将每个覆盖的胞嘧啶的C和T计数相加，并返回每个窗口的总C和T计数。
tiles=tileMethylCounts(myobj,win.size=1000,step.size=1000)
head(tiles[[1]],3)
## chr start end strand coverage numCs numTs
## 1 chr21 9764001 9765000 * 24 3 21
## 2 chr21 9820001 9821000 * 13 0 13
## 3 chr21 9837001 9838000 * 11 0 11
差异分析
这是甲基化分析中十分重要的步骤，calculateDiffMeth()功能是计算差异甲基化的主要功能。根据每组样本的大小，将使用Fisher精确或逻辑回归来计算P值。如果有重复项，该函数将自动使用逻辑回归。用户如果合并重复项，则可以强制该函数使用费舍尔精确测试。这可以通过pool()功能来实现。
#该函数使用了上面读取甲基化文件时的分组信息，再进行分析前请确保分组信息#与甲基化文件一一对应
myDiff=calculateDiffMeth(meth，mc.cores=2， #多核并行操作)
## two groups detected:
## will calculate methylation difference as the difference of
## treatment (group: 1) - control (group: 0)
经过q值计算后，我们可以基于q值和甲基化差异阈值百分比选择差异甲基化位点。接下来选择q值<0.05且甲基化百分比差异大于20％的碱基。这是筛选差异甲基化位点的一般选择，当然也可以使用更加严格的参数。如果指定type="hyper"或type="hypo"选项，将获得高甲基化或低甲基化的位点。
# get hyper methylated bases
myDiff25p.hyper=getMethylDiff(myDiff,difference=20,qvalue=0.05,type="hyper")
# get hypo methylated bases
myDiff25p.hypo=getMethylDiff(myDiff,difference=20,qvalue=0.05,type="hypo")
# get all differentially methylated bases
myDiff25p=getMethylDiff(myDiff,difference=20,qvalue=0.05)
我们还可以使用以下函数可视化每个染色体的低/高甲基化位点的分布。在这种情况下，展示每条染色体中高/低甲基化区域的数量以及占比，由于示例集仅包含一个染色体，仅显示一条染色体的信息。
diffMethPerChr(myDiff,plot=FALSE,qvalue.cutoff=0.01, meth.cutoff=25)
## $diffMeth.per.chr
## chr number.of.hypermethylated percentage.of.hypermethylated
## 1 chr21 75 7.788162
## number.of.hypomethylated percentage.of.hypomethylated
## 1 59 6.126687
##
## $diffMeth.all
## number.of.hypermethylated percentage.of.hypermethylated
## 1 75 7.788162
## number.of.hypomethylated percentage.of.hypomethylated
## 1 59 6.126687
差异甲基化位点的注释
我们可以使用
genomation包基于基因注释来注释差异甲基化区域/碱基。在此示例中，我们从BED文件中读取了基因注释，并使用genomation函数用该信息注释了差异甲基化区域。请注意，这些函数对GRanges对象起作用，因此我们首先将methylKit对象强制为GRanges。此注释操作将告诉我们启动子/内含子/外显子/基因间区域上差异甲基化区域的百分比。在这种情况下，我们从BED文件中读取注释，可以使用GenomicFeatures可从Bioconductor.org 获得的软件包或其他软件包来获取相似的基因注释信息。我们可以在USCS等网站下载所需物种的bed文件，随后使用
genomation将其转换为可注释的对象，也可以使用TxDb系列包自己转换，需要注意，进行注释的GRangesList对象必须包括启动子/内含子/外显子/基因间区域上的信息。
# read the gene BED file
gene.obj=readTranscriptFeatures(system.file("extdata", "refseq.hg18.bed.txt",
package = "methylKit"))
# annotate differentially methylated CpGs with
# promoter/exon/intron using annotation data
annotateWithGeneParts(as(myDiff25p,"GRanges"),gene.obj)
#同样，我们可以使用CpG岛注释文件，并用它们注释差异甲基化的碱基/区域。
# read the shores and flanking regions and name the flanks as shores
# and CpG islands as CpGi
cpg.obj=readFeatureFlank(system.file("extdata", "cpgi.hg18.bed.txt",
package = "methylKit"),
feature.flank.name=c("CpGi","shores"))
# convert methylDiff object to GRanges and annotate
diffCpGann=annotateWithFeatureFlank(as(myDiff25p,"GRanges"),
cpg.obj$CpGi,cpg.obj$shores,
feature.name="CpGi",flank.name="shores")
其他功能
得到差异甲基化区域的注释后，我们可以使用getAssociationWithTSS函数来获得到每个差异甲基化位点到TSS的距离和最近的基因名。
diffAnn=annotateWithGeneParts(as(myDiff25p,"GRanges"),gene.obj,intersect.chr = TRUE）
#intersect.chr参数指定是否交替正负链以TSS为基准展示结果
# target.row is the row number in myDiff25p
head(getAssociationWithTSS(diffAnn))
#获得与内含子/外显子/启动子重叠的差异甲基化区域的百分比/数目。
getTargetAnnotationStats(diffAnn,percentage=TRUE,precedence=TRUE)
#绘制饼图查看与基因组各区域重叠的差异甲基化区域占比
plotTargetAnnotation(diffAnn,precedence=TRUE,
main="differential methylation annotation")
作者提供了便利功能可以在分析过程中任意转换对象的形式，包括methylKit对象，methylDB对象和GRanges对象，以保证数据的无缝衔接
class(meth)
## [1] "methylBase"
## attr(,"package")
## [1] "methylKit"
#转换为GRanges对象
as(meth,"GRanges")
#将methylKit对象转换为methylDB对象，反之亦然
as(myobj[[1]],"methylRaw")
