DSS (Dispersion Shrinkage for Sequencing data)绣的，為基于高通量測序數(shù)據(jù)的差異分析而設(shè)計的Bioconductor包俩滥。主要應用于BS-seq（亞硫酸氫鹽測序）中計算不同組別間差異甲基化位點（DML）和差異甲基化區(qū)域（DMR）即Call DML or DMR网持。
Bisulfite Sequencing (BS-Seq)上游測序數(shù)據(jù)可以得到甲基化位點的信息尼斧，而后續(xù)DML以及DMR的確定以及可視化就需要DSS包荚恶。
DSS包的使用主要包括：輸入文件的準備 --> 利用DMLtest函數(shù)call DML --> 利用callDML函數(shù)Call DML --> 利用callDMR函數(shù)Call DMR --> 利用showOneDMR函數(shù)對DMRs可視化

1.輸入文件準備

DSS包要求輸入文件數(shù)據(jù)的格式如下:
每一行代表一個CpG site

Below shows an example of a small part of such a file:

chr pos N X
chr18 3014904 26 2
chr18 3031032 33 12
chr18 3031044 33 13
chr18 3031065 48 24

• 第一列為染色體
• 第二列為位置
• 第三列為total reads
• 第四列為甲基化的reads

拿到上游比對結(jié)果后需要把結(jié)果文件*.bismark.cov.gz改成DSS包所要求的樣子撩穿，使用Linux或者R進行簡單的處理及可得到input文件。

2. 計算不同組別間差異甲基化位點和區(qū)域—Call DML or DMR

DML：甲基化差異位點谒撼；DMR：甲基化差異區(qū)域
使用DSS包自帶的數(shù)據(jù)演示如何計算不同組別間差異甲基化位點和區(qū)域

2.1 載入DSS和bsseq包構(gòu)建BSobj對象

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)
BSobj <- makeBSseqData( list(dat1.1, dat1.2, dat2.1, dat2.2),
                        c("C1","C2", "N1", "N2") )[1:1000,]

> BSobj
An object of type 'BSseq' with
1000 methylation loci
4 samples
has not been smoothed
All assays are in-memory

2.2 利用DMLtest函數(shù)call DML

DML：甲基化差異位點食寡；DMR：甲基化差異區(qū)域

DMLtest函數(shù)主要包括以下步驟：

• 計算所有CpG位點的平均甲基化水平;
• 計算每個CpG位點的分散度dispersions；
• 進行沃爾德檢驗 conduct Wald test

在第一步過程中嗤栓，我們可以選擇是否smoothing處理甲基化水平冻河。當測序結(jié)果中CpG 位點特別密集時（比如：whole-genome BS-seq得到的數(shù)據(jù)）smoothing處理可以以更簡潔直接的方式幫助估算平均甲基化水平箍邮；當CpG 位點比較稀疏時（比如：RRBS or hydroxyl-methylation得到的數(shù)據(jù)）則不需要smoothing處理。

Call DML時不經(jīng)過smoothing處理：

# To perform DML test without smoothing, do:
dmlTest <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"))

> head(dmlTest)
    chr     pos       mu1       mu2        diff    diff.se       stat        phi1       phi2       pval       fdr
1 chr18 3014904 0.3817233 0.4624549 -0.08073162 0.24997034 -0.3229648 0.300542998 0.01706260 0.74672190 0.9985094
2 chr18 3031032 0.3380579 0.1417008  0.19635711 0.11086362  1.7711592 0.008911745 0.04783892 0.07653423 0.6792127
3 chr18 3031044 0.3432172 0.3298853  0.01333190 0.12203116  0.1092500 0.010409029 0.01994821 0.91300423 0.9985094
4 chr18 3031065 0.4369377 0.3649218  0.07201587 0.10099395  0.7130711 0.010320888 0.01603200 0.47580174 0.9985094
5 chr18 3031069 0.2933572 0.5387464 -0.24538920 0.13178800 -1.8619996 0.012537553 0.02320887 0.06260315 0.6158797
6 chr18 3031082 0.3526311 0.3905718 -0.03794068 0.07847999 -0.4834440 0.007665696 0.01145531 0.62878051 0.9985094

Call DML時經(jīng)過smoothing處理代碼：

# To perform statistical test for DML with smoothing, do:
dmlTest.sm <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"), smoothing=TRUE)

> head(dmlTest.sm)
    chr     pos       mu1       mu2        diff    diff.se       stat       phi1       phi2      pval       fdr
1 chr18 3014904 0.3693669 0.4566563 -0.08728939 0.29967322 -0.2912819 0.30054300 0.01706260 0.7708357 0.9656515
2 chr18 3031032 0.3433882 0.3679732 -0.02458503 0.03970109 -0.6192533 0.03177894 0.28323422 0.5357495 0.8639036
3 chr18 3031044 0.3412867 0.3678807 -0.02659404 0.04032823 -0.6594397 0.02536938 0.02080295 0.5096134 0.8596522
4 chr18 3031065 0.3358830 0.3511983 -0.01531533 0.04799161 -0.3191252 0.01123412 0.01621926 0.7496316 0.9652417
5 chr18 3031069 0.3358830 0.3511983 -0.01531533 0.03205500 -0.4777830 0.02832889 0.05857316 0.6328047 0.8968029
6 chr18 3031082 0.3358830 0.3511983 -0.01531533 0.05846593 -0.2619531 0.01682981 0.01368466 0.7933576 0.9745116

2.3 利用callDML函數(shù)call DML

使用callDML函數(shù)call DML叨叙，結(jié)果可以按顯著性排序：

dmls <- callDML(dmlTest, p.threshold=0.001)
> head(dmls)
      chr     pos        mu1       mu2       diff    diff.se      stat        phi1       phi2         pval          fdr
450 chr18 3976129 0.01027497 0.9390339 -0.9287590 0.06544340 -14.19179 0.052591567 0.02428826 1.029974e-45 2.499403e-43
451 chr18 3976138 0.01027497 0.9390339 -0.9287590 0.06544340 -14.19179 0.052591567 0.02428826 1.029974e-45 2.499403e-43
638 chr18 4431501 0.01331553 0.9430566 -0.9297411 0.09273779 -10.02548 0.053172411 0.07746835 1.177826e-23 1.429096e-21
639 chr18 4431511 0.01327049 0.9430566 -0.9297862 0.09270080 -10.02997 0.053121697 0.07746835 1.125518e-23 1.429096e-21
710 chr18 4564237 0.91454619 0.0119300  0.9026162 0.05260037  17.15988 0.009528898 0.04942849 5.302004e-66 3.859859e-63
782 chr18 4657576 0.98257334 0.0678355  0.9147378 0.06815000  13.42242 0.010424723 0.06755651 4.468885e-41 8.133371e-39
    postprob.overThreshold
450                      1
451                      1
638                      1
639                      1
710                      1
782                      1

默認情況下锭弊，計算基于零假設(shè)，即默認甲基化水平的差異為0擂错。當然味滞，我們可以指定差異的閾值，只有差異大于閾值（0.1）的才會被call出來：

# To detect loci with difference greater than 0.1, do:
> dmls2 <- callDML(dmlTest, delta=0.1, p.threshold=0.001)
> head(dmls2)
      chr     pos        mu1       mu2       diff    diff.se      stat        phi1       phi2         pval
450 chr18 3976129 0.01027497 0.9390339 -0.9287590 0.06544340 -14.19179 0.052591567 0.02428826 1.029974e-45
451 chr18 3976138 0.01027497 0.9390339 -0.9287590 0.06544340 -14.19179 0.052591567 0.02428826 1.029974e-45
638 chr18 4431501 0.01331553 0.9430566 -0.9297411 0.09273779 -10.02548 0.053172411 0.07746835 1.177826e-23
639 chr18 4431511 0.01327049 0.9430566 -0.9297862 0.09270080 -10.02997 0.053121697 0.07746835 1.125518e-23
710 chr18 4564237 0.91454619 0.0119300  0.9026162 0.05260037  17.15988 0.009528898 0.04942849 5.302004e-66
782 chr18 4657576 0.98257334 0.0678355  0.9147378 0.06815000  13.42242 0.010424723 0.06755651 4.468885e-41
             fdr postprob.overThreshold
450 2.499403e-43                      1
451 2.499403e-43                      1
638 1.429096e-21                      1
639 1.429096e-21                      1
710 3.859859e-63                      1
782 8.133371e-39                      1

2.4 利用callDMR函數(shù)Call DMR

DML：甲基化差異位點钮呀；DMR：甲基化差異區(qū)域
甲基化差異區(qū)域檢測也是基于差異位點的結(jié)果剑鞍，同樣使用callDML函數(shù)。當不同組別間CpG位點區(qū)域具有顯著的統(tǒng)計學差異時這段差異區(qū)域被定義為DMRs爽醋。

# Call DMR by using callDMR function
##Regions with many statistically significant CpG sites are identified as DMRs.
dmrs <- callDMR(dmlTest, p.threshold=0.01)
> head(dmrs)
     chr   start     end length nCG meanMethy1 meanMethy2 diff.Methy areaStat
27 chr18 4657576 4657639     64   4   0.506453   0.318348   0.188105 14.34236

同理蚁署，這里我們也可以使用delta參數(shù)以及調(diào)整p.threshold指定差異的閾值：

# To detect regions with difference greater than 0.1, do:
dmrs2 <- callDMR(dmlTest, delta=0.1, p.threshold=0.05)
> head(dmrs2)
     chr   start     end length nCG meanMethy1 meanMethy2 diff.Methy areaStat
31 chr18 4657576 4657639     64   4  0.5064530  0.3183480   0.188105 14.34236
19 chr18 4222533 4222608     76   4  0.7880276  0.3614195   0.426608 12.91667

這里我們需要注意，選擇一個合理的閾值來定義DMRs是非常困難的蚂四，所以建議嘗試不同的閾值光戈，以獲得滿意的結(jié)果。

2.5 可視化

使用showOneDMR函數(shù)可視化甲基化差異區(qū)域DML遂赠，該函數(shù)不僅可以繪制甲基化所占百分比還可以繪制每個CpG位點的覆蓋深度久妆。

showOneDMR(dmrs[1,], BSobj)

From RRBS experiment so the DMRs are much shorter

我們的示例數(shù)據(jù)來自RRBS實驗結(jié)果，所以甲基化差異區(qū)域DML很短跷睦。一般whole-genome BS-seq數(shù)據(jù)中DML會長一些：

whole-genome BS-seq

代碼純享版：

# 1. Load in library. Read in text files and create an object of BSseq class
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)
BSobj <- makeBSseqData( list(dat1.1, dat1.2, dat2.1, dat2.2),
                        c("C1","C2", "N1", "N2") )[1:1000,]
BSobj

# 2.Perform statistical test for DML by calling DMLtest function.
## To perform DML test without smoothing, do:
dmlTest <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"))
head(dmlTest)
## To perform statistical test for DML with smoothing, do:
dmlTest.sm <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"), smoothing=TRUE)
head(dmlTest.sm)

# 3.Call DML by using callDML function. The results DMLs are sorted by the significance.
dmls <- callDML(dmlTest, p.threshold=0.001)
head(dmls)
##To detect loci with difference greater than 0.1, do:
dmls2 <- callDML(dmlTest, delta=0.1, p.threshold=0.001)
head(dmls2)

# 4.Call DMR by using callDML function
##Regions with many statistically significant CpG sites are identified as DMRs.
dmrs <- callDMR(dmlTest, p.threshold=0.01)
head(dmrs)
##To detect regions with difference greater than 0.1, do:
dmrs2 <- callDMR(dmlTest, delta=0.1, p.threshold=0.05)
head(dmrs2)

# 5.The DMRs can be visualized using showOneDMR function
showOneDMR(dmrs[1,], BSobj)