bsseq 进行差异甲基化分析
bsseq
主要用来分析WGBS的数据, 安装过程如下
1 2 | source(“http://bioconductor.org/biocLite.R“) biocLite(“bsseq”) |
bsseq
的分析主要包括以下4步:
- 读取原始数据
- BSmooth
- t-test检验
- DMR
1. 读取原始数据
bsseq
要求的原始数据格式如下:
共6列数据,制表符分隔,每一行代表一个甲基化位点,前5列很好理解,描述甲基化位点的染色体位置和类别,默认情况下bbseq
用于分析CpG
类型的甲基化位点。当然其他类型的数据,比如CHG
, CHH
也支持,但是需要调整参数。Cov
代表覆盖到这个位点的reads数,M
代表其中发生了甲基化的reads数目。
每个样本一个这样的原始数据,用来表示该样本methylation calling
的结果,这样的数据我们从bismark
的结果中也可以得到。当原始数据准备好之后,首选需要读取所有样本的原始数据,然后导入到R中,生成一个bbseq
定义的对象。在bbseq
安装的路径下,提供了一个名为get_BS.chr22.R
的脚本,展示了如何从读取所有样本原始数据的过程。
代码如下
这里以mc_imr90_r1_22
和mc_imr90_r2_22
两个样本的原始数据为例,详细展示了读取过程。我们只需要根据自己的数据,适当修改上述代码就可以了。主要注意sampleNames
和pData
数据就可以了。
2 Smooth
已测试数据BS.chr22
为例,smooth
的过程如下
在实际分析中,由于甲基化位点很多,所以这一步时间特别久,为了提高速度,可以添加mc.cores
参数,这个参数指定了CPU个数,用于并行执行。
1 | BS.chr22.1 <- BSmooth(BS.chr22, mc.cores = 2, verbose = TRUE) |
3. T-test
在分析之前,有必要过滤掉覆盖度较低的甲基化位点。通常保留在所有样本中覆盖度大于2的甲基化位点,但是也可以修改这个条件。过滤之后,直接通过BSmooth.tstat
进行分析
下面的代码基于bsseqData
包中的数据,这个数据包含了6个样本,分为normal
和cancer
两组
group1
指定属于treatment
组的样本,group2
指定属于control
组的样本。
4. DMR
通过dmrFinder
函数进行差异甲基化分析, 代码如下:
cutoff
指定DMR
的阈值,这个阈值根据t-test
的结果进行调整。subset
对差异甲基化的结果进行筛选,筛选包含甲基化位点个数大于3而且meanDiff
大于0.1的甲基化区域。
相关推荐:
最新创建圈子
-
原料药研发及国内外注册申报
2019-01-25 10:41圈主:caolianhui 帖子:33 -
制药工程交流
2019-01-25 10:40圈主:polysciences 帖子:30 -
健康管理
2019-01-25 10:40圈主:neuromics 帖子:20 -
发酵技术
2019-01-25 10:39圈主:fitzgerald 帖子:17 -
医学肿瘤学临床试验
2019-01-25 10:39圈主:bma 帖子:58