参考来源1：https://www.jianshu.com/p/5d7e550abc1a
参考来源2：https://www.jianshu.com/p/971ff2ee533a
参考来源3：https://www.plob.org/article/10843.html
第三部分 实战篇之差异分析
1. 前言
衔接上一部分数据比对后的结果，使用R包DSS进行处理。
2. Input文件准备
我们先来复习一下上一部分中得到的数据结果： *.bismark.cov.gz 文件
$head test_data_bismark_bt2.deduplicated.bismark.cov
Chr1 975476 975476 100 1 0
Chr1 975488 975488 100 1 0
Chr1 975490 975490 100 1 0
Chr1 2224487 2224487 100 1 0
Chr1 2224489 2224489 100 1 0
Chr1 2224514 2224514 100 1 0
Chr1 2224520 2224520 100 1 0
Chr1 2313220 2313220 0 0 1
Chr1 9313902 9313902 100 1 0
Chr1 9313914 9313914 100 1 0
.bismark.cov文件记录了每个甲基化位点的覆盖度，包含发生甲基化的reads数，未发生甲基化的reads数以及甲基化频率。
# 第一列代表chromosome
# 第二，三列代表location
# 第四列代表甲基化百分比
# 第五列代表发生甲基化的reads数
# 第六列代表未发生甲基化的reads数
3. 使用Bioconductor中的DSS包进行差异分析
3.1 DSS包介绍
这里我们使用的R包为DSS，使用Bioconductor进行安装。
这个包可以对甲基化数据做两件事：
- 找出差异甲基化位点
- 根据甲基化位点找出差异化甲基化区域
DSS包对差异甲基化的检测基于β负二项分布的严格沃尔德检验。
3.2 DSS包的安装
# R version 3.5.2 (2018-12-20) -- "Eggshell Igloo"
if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
BiocManager::install("GenomicFeatures")
BiocManager::install("R.methodsS3")
BiocManager::install("R.oo")
BiocManager::install("bsseq")
BiocManager::install("DSS")
library("DSS")
3.3 DSS包的输入格式
DSS包可以对无重复数据进行处理
输入数据的格式如下:
-
第一列为染色体
-
第二列为位置
-
第三列为total reads
-
第四列为甲基化的reads
chr pos N X
chr18 3014904 26 2
chr18 3031032 33 12
chr18 3031044 33 13
chr18 3031065 48 24
# 根据这个特性，我们可以将上面的输出文件通过Linux或者R进行简单的处理得到input文件。
4. Call DML or DMR
DML是甲基化差异位点，DMR为甲基化差异区域。
使用DSS包自带的数据进行演示
4.1 载入两个包DSS和bsseq（需要该包构建obj对象）
setwd("~/GSE81436/BS-seq/DSS_data/Cr-DJ_osddm1a1b")
library(DSS)
require(bsseq)
dat1 <- read.table("Cr-DJ_DSS.txt", header=TRUE)
dat2 <- read.table("osddm1a1b_DSS.txt", header=TRUE)
BSobj <- makeBSseqData( list(dat1, dat2),c("C","N") )
BSobj
# 查看BSobj的类型
BSobj
An object of type 'BSseq' with
1000 methylation loci
4 samples
has not been smoothed
All assays are in-memory、
4.2 根据甲基化水平进行loci的差异分析
利用DMLtest函数call DML，分为以下几个步骤：
-
预测所有CpG位点的平均甲基化水平
-
对每个CpG位点估算其甲基化水平的dispersions
-
进行沃尔德检验
在第一步过程中，smoothing可以更好帮助估算甲基化（针对whole-genome BS-seq）。
而RRBS不需要smoothing。
根据甲基化水平进行loci的差异分析：
dmls <- callDML(dmlTest.sm, p.threshold=0.001)
head(dmls)
dmls2 <- callDML(dmlTest.sm, delta=0.1, p.threshold=0.001)
head(dmls2)
4.3 根据dmlTest.sm来call DML
dmls <- callDML(dmlTest.sm, p.threshold=0.001)
head(dmls)
# 用户也可以指定差异的阈值，只有差异大于阈值的才会被call出来
dmls2 <- callDML(dmlTest.sm, delta=0.1, p.threshold=0.001)
head(dmls2)
注意：smoothing前后得到的结果差异还是很大，所以针对不同的实验类型我们需要注意是否使用smoothing。同理，也可以使用的delta参数以及调整p.threshold得到合适的结果。
4.4 根据dmlTest Call DMR
dmrs <- callDMR(dmlTest.sm, p.threshold=0.01)
head(dmrs)
# 用户也可以指定差异的阈值，只有差异大于阈值的才会被call出来
dmrs2 <- callDMR(dmlTest.sm, delta=0.1, p.threshold=0.05)
head(dmrs2)
5. 可视化
DSS包提供了一个不是很美观的可视化函数，用户其实可以使用coverage结果在R里面作图。
showOneDMR(dmrs[1,], BSobj)
6. 结果输出
write.table(dmls,file="Cr-DJ_osddm1a1b.dmls.part.txt",sep="\t")
write.table(dmls2,file="Cr-DJ_osddm1a1b.dmls.part2.txt",sep="\t")
write.table(dmrs,file="Cr-DJ_osddm1a1b.dmrs.part.txt",sep="\t")
write.table(dmrs2,file="Cr-DJ_osddm1a1b.dmrs.part2.txt",sep="\t")
