DSS是一个R包，对基于计数的测序数据进行差异分析。 它可以检测来自RNA-seq的差异表达基因（DEG），以及来自亚硫酸氢盐测序（BS-seq）的差异甲基化位点或区域（DML / DMR）。
DSS (Dispersion Shrinkage for Sequencing data)，为基于高通量测序数据的差异分析而设计的Bioconductor包。主要应用于BS-seq（亚硫酸氢盐测序）中计算不同组别间差异甲基化位点（DML）和差异甲基化区域（DMR）即Call DML or DMR。

使用R包DSS多种方式检验差异甲基化信号区域。
这里我使用的R版本为3.6.0,安装“DSS”包

BiocManager::install("DSS")

1.准备输入文件
DSS包要求输入数据格式：每一行代表一个CpG位点，格式如下：
第一列为染色体
第二列为位置
第三列为总reads数
第四列为甲基化的reads数

我们看一下上一步我们得到的数据test_R1_bismark_bt2_pe.deduplicated.bismark.cov.gz

less  test_R1_bismark_bt2_pe.deduplicated.bismark.cov.gz

image.jpeg

在file.cov.gz中
第一列数据为染色体
第二、三列数据为甲基化位置
第四列为甲基化百分比
第五列为甲基化数目
第六列为未甲基化数目

在这里,需要对我们的数据利用R，linux或者python整理成DSS包所需输入文件格式。这里我使用python整理的输入文件格式。

image.jpeg

2.对不同组别之间DNA甲基化进行差异分析
这里我们使用DSS包自带的数据
2.1 加载DSS包，并读取甲基化数据

library(DSS)
require(bsseq)

###读取甲基化数据
path <- file.path(system.file(package="DSS"), "extdata")
dat1.1 <- read.table(file.path(path, "cond1_1.txt"), header=TRUE)
dat1.2 <- read.table(file.path(path, "cond1_2.txt"), header=TRUE)
dat2.1 <- read.table(file.path(path, "cond2_1.txt"), header=TRUE)
dat2.2 <- read.table(file.path(path, "cond2_2.txt"), header=TRUE)

2.2构建BSobj对象

BSobj <- makeBSseqData(list(dat1.1, dat1.2, dat2.1, dat2.2),
                        c("C1","C2", "N1", "N2"))

BSobj

image.jpeg

可以看出我们的数据包括4个样本，34739个甲基化位点。

2.3利用DMLtest函数call DML
For whole-genome BS-seq data, perform DML test with smoothing
smoothing: 用于指示是否在估计平均甲基化水平时应用平滑的标志,这里我理解的smoothing有点类似滑动窗口的意思。我们也可以选择smoothing=False

dmlTest.sm <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"), smoothing=TRUE)
head(dmlTest.sm)

image.jpeg

2.4利用callDML函数可以找出差异甲基化位点
delta:定义DML的阈值。在DML检验程序中，在每个CpG位点进行两组均值相等的假设检验。这里如果指定了delta，函数将计算均值差大于delta的后验概率，然后基于此调用DML。
p.threshold: 当未指定delta时，这是定义DML的p值阈值，例如p值小于该阈值的位点将被视为DML。当指定delta时，后验概率大于1-p阈值的CpG位点被视为DML。

dmls <- callDML(dmlTest.sm,delta=0.1, p.threshold=0.001)
head(dmls)

image.jpeg

2.5利用callDMR函数可以找出差异甲基化区域
当不同组别间CpG位点区域具有显著的统计学差异时这段差异区域被定义为DMRs。DMR也是基于DML被检测出来的。

dmrs <- callDMR(dmlTest.sm, delta=0.1, p.threshold=0.01)
head(dmrs)

image.jpeg

如果想把甲基化位点和甲基化区域信息输出，可通过write.csv()输出。

write.csv(dmrs,file= “sili.csv”)

2.6 #使用showOneDMR函数可视化DMR
给定一个DMR和一个BSseq对象，此函数将生成一个多面板图，每个图对应一个示例文件，以可视化甲基化水平。每个CpG上都有一个条，灰色线条表示总覆盖率，蓝色线条表示甲基化水平。

showOneDMR(dmrs[1,], BSobj)

image.jpeg

showOneDMR唯一缺陷只能展示一个区域甲基化水平，有点丑，感兴趣的同学可以自己利用脚本画一下全基因组甲基化水平。

参考：
1.https://www.rdocumentation.org/packages/DSS/versions/2.12.0

2. 计算不同组别间差异甲基化位点和区域的R包—DSS
   https://www.jianshu.com/p/203ac75c0c32?utm_campaign=maleskine&utm_content=note&utm_medium=seo_notes&utm_source=recommendation