BS-Seq(4) –差异分析

前言

衔接上一篇数据比对后的结果，使用R包DSS进行处理。

Input文件准备

我们先来复习一下上一节课得到的数据结果：
*.bismark.cov.gz 文件

$head test_data_bismark_bt2.deduplicated.bismark.cov
chr1    10542   10542   0   0   1
chr1    10563   10563   100 2   0
chr1    10571   10571   100 2   0
chr1    10577   10577   100 2   0
chr1    10579   10579   100 2   0
# 第一列代表chromosome
# 第二，三列代表location
# 第四列代表甲基化百分比
# 第五列代表甲基化数目
# 第六列代表未甲基化数目

这里我们使用的R包为DSS，使用Bioconductor进行安装。
这个包可以对甲基化数据做两件事：

• 找出差异甲基化位点
• 根据甲基化位点找出差异化甲基化区域
  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文件。
这里我是这样转换的

zcat test_data_bismark_bt2.deduplicated.bismark.cov.gz | awk '{print $1"/t"$2"/t"$5+$6"/t"$5}' > test_data.dss.input.txt
head test_data.dss.input.txt
chr1    10542   1   0
chr1    10563   2   2
chr1    10571   2   2
chr1    10577   2   2
chr1    10579   2   2
然后手动在第一列添加chr pos N X，变成
chr pos N X
chr1    10542   1   0
chr1    10563   2   2

Call DML or DMR

DML是甲基化差异位点，DMR为甲基化差异区域。
使用DSS包自带的数据进行演示

1. 载入两个包DSS和bsseq（需要该包构建obj对象）

BiocManager::install("GenomicFeatures")
BiocManager::install("R.methodsS3")
BiocManager::install("R.oo")
BiocManager::install("bsseq")
BiocManager::install("DSS")


library(DSS)
library(stats)
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对象
BSobj <- makeBSseqData( list(dat1.1, dat1.2, dat2.1, dat2.2),
                           c("C1","C2", "N1", "N2") )
# 查看BSobj的类型
BSobj 
An object of type 'BSseq' with
  1000 methylation loci
  4 samples
has not been smoothed
All assays are in-memory
可以看出我们的数据包括4个样本，34739个甲基化位点。

2.利用DMLtest函数call DML

预测所有CpG位点的平均甲基化水平
对每个CpG位点估算其甲基化水平的dispersions
进行沃尔德检验
在第一步过程中，smoothing可以更好帮助估算甲基化（针对whole-genome BS-seq,而RRBS不需要smoothing。
smoothing有点类似滑动窗口的意思。我们也可以选择smoothing=False
根据甲基化水平进行loci的差异分析

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

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

3. 利用callDML函数可以找出差异甲基化位点

dmls <- callDML(dmlTest.sm, p.threshold=0.001)
dmls
      chr     pos       mu1       mu2      diff    diff.se     stat       phi1       phi2         pval          fdr postprob.overThreshold
447 chr18 3973699 0.8530694 0.3432547 0.5098147 0.05726793 8.902272 0.03662109 0.98947516 5.471481e-19 4.962634e-16              1.0000000
709 chr18 4564190 0.7773858 0.1977036 0.5796822 0.10294267 5.631116 0.21725415 0.02952656 1.790467e-08 5.413180e-06              0.9999984
710 chr18 4564237 0.7773858 0.1977036 0.5796822 0.15431085 3.756587 0.02931516 0.26238558 1.722462e-04 7.439396e-03              0.9990652
dmls <- callDML(dmlTest, p.threshold=0.001)
dmls
      chr     pos        mu1         mu2       diff    diff.se       stat        phi1       phi2         pval          fdr postprob.overThreshold
450 chr18 3976129 0.01027497 0.939033927 -0.9287590 0.06544340 -14.191789 0.052591567 0.02428826 1.029974e-45 2.499403e-43              1.0000000
451 chr18 3976138 0.01027497 0.939033927 -0.9287590 0.06544340 -14.191789 0.052591567 0.02428826 1.029974e-45 2.499403e-43              1.0000000
638 chr18 4431501 0.01331553 0.943056638 -0.9297411 0.09273779 -10.025483 0.053172411 0.07746835 1.177826e-23 1.429096e-21              1.0000000
639 chr18 4431511 0.01327049 0.943056638 -0.9297862 0.09270080 -10.029969 0.053121697 0.07746835 1.125518e-23 1.429096e-21              1.0000000
710 chr18 4564237 0.91454619 0.011930005  0.9026162 0.05260037  17.159883 0.009528898 0.04942849 5.302004e-66 3.859859e-63              1.0000000
782 chr18 4657576 0.98257334 0.067835497  0.9147378 0.06815000  13.422418 0.010424723 0.06755651 4.468885e-41 8.133371e-39              1.0000000
582 chr18 4340682 0.95398081 0.030390730  0.9235901 0.10935874   8.445508 0.085494283 0.04540643 3.027264e-17 2.754810e-15              1.0000000
583 chr18 4340709 0.95398081 0.030390730  0.9235901 0.10935874   8.445508 0.085494283 0.04540643 3.027264e-17 2.754810e-15              1.0000000
340 chr18 3542732 0.95023554 0.034383112  0.9158524 0.11937407   7.672122 0.089137013 0.04474741 1.691739e-14 1.368429e-12              1.0000000
395 chr18 3723448 0.06570765 0.751990744 -0.6862831 0.09825286  -6.984866 0.011958092 0.01646418 2.851278e-12 2.075730e-10              1.0000000
188 chr18 3370113 0.01488553 0.787980174 -0.7730946 0.11380172  -6.793347 0.054190769 0.02752024 1.095611e-11 7.250956e-10              1.0000000
400 chr18 3785543 0.79337804 0.118353679  0.6750244 0.12183251   5.540593 0.017720841 0.02442007 3.014493e-08 1.688116e-06              0.9999988
683 chr18 4494490 0.77615275 0.104235359  0.6719174 0.12407577   5.415380 0.009219023 0.08210742 6.115879e-08 3.180257e-06              0.9999980
783 chr18 4657592 0.96858371 0.266894590  0.7016891 0.13077255   5.365722 0.064228633 0.07721117 8.062618e-08 3.668491e-06              0.9999979
642 chr18 4431618 0.43034706 0.965539700 -0.5351926 0.09578178  -5.587625 0.012064922 0.07497569 2.301966e-08 1.396526e-06              0.9999972
189 chr18 3370141 0.01488553 0.676293642 -0.6614081 0.12554112  -5.268458 0.054190769 0.02138845 1.375746e-07 5.891428e-06              0.9999961
738 chr18 4635185 0.44671015 0.007342128  0.4393680 0.08148745   5.391849 0.010069744 0.05154034 6.973627e-08 3.384534e-06              0.9999844
330 chr18 3542403 0.02875655 0.714168320 -0.6854118 0.15294560  -4.481409 0.041892077 0.03368033 7.415192e-06 2.841190e-04              0.9999354
185 chr18 3347936 0.49560648 0.017291187  0.4783153 0.10184167   4.696656 0.012173893 0.04533502 2.644552e-06 1.069574e-04              0.9998983
92  chr18 3217241 0.03876653 0.767643956 -0.7288774 0.17547605  -4.153714 0.038095316 0.03753383 3.271214e-05 9.602621e-04              0.9998319
396 chr18 3723468 0.23043341 0.951779714 -0.7213463 0.18220738  -3.958930 0.173291593 0.07965750 7.528626e-05 1.838513e-03              0.9996786
40  chr18 3047980 0.07206183 0.658535570 -0.5864737 0.14396467  -4.073734 0.012517547 0.02615351 4.626536e-05 1.247451e-03              0.9996373
614 chr18 4353399 0.01355584 0.547298573 -0.5337427 0.12855458  -4.151876 0.053422013 0.02055735 3.297603e-05 9.602621e-04              0.9996300
99  chr18 3226251 0.97443783 0.217343279  0.7570945 0.19807591   3.822244 0.012911947 0.27481269 1.322425e-04 3.105566e-03              0.9995532
613 chr18 4353379 0.01355584 0.529341634 -0.5157858 0.13033377  -3.957422 0.053422013 0.02080043 7.576289e-05 1.838513e-03              0.9992902
894 chr18 4921874 0.35815864 0.009172558  0.3489861 0.07955244   4.386868 0.009774839 0.05010151 1.149944e-05 3.986472e-04              0.9991255
895 chr18 4921881 0.35815864 0.009172558  0.3489861 0.07955244   4.386868 0.009774839 0.05010151 1.149944e-05 3.986472e-04              0.9991255

delta:定义DML的阈值。在DML检验程序中，在每个CpG位点进行两组均值相等的假设检验。这里如果指定了delta，函数将计算均值差大于delta的后验概率，然后基于此调用DML。

p.threshold: 当未指定delta时，这是定义DML的p值阈值，例如p值小于该阈值的位点将被视为DML。当指定delta时，后验概率大于1-p阈值的CpG位点被视为DML。
当然，用户也可以指定差异的阈值，只有差异大于阈值的才会被call出来。

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

4. 利用callDMR函数可以找出差异甲基化区域

dmrs <- callDMR(dmlTest.sm, p.threshold=0.01)
dmrs
     chr   start     end length nCG meanMethy1 meanMethy2 diff.Methy areaStat
23 chr18 4921637 4922059    423  27 0.11002940 0.01809674 0.09193266 86.29588
7  chr18 3507919 3508022    104  10 0.07524915 0.03294316 0.04230600 30.55943
15 chr18 4340682 4340753     72   4 0.89237955 0.35052968 0.54184987 10.83526
dmrs <- callDMR(dmlTest, p.threshold=0.01)
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

我们可以发现，smoothing前后得到的结果差异还是很大，所以针对不同的实验类型我们需要注意是否使用smoothing。
同理，也可以使用的delta参数以及调整p.threshold得到合适的结果。

dmrs2 <- callDMR(dmlTest, delta=0.1, p.threshold=0.05)
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

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

write.csv(dmrs,file= "dmrs.csv")
write.csv(dmls,file= "dmls.csv")

5.使用showOneDMR函数可视化DMR

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

showOneDMR(dmrs[1,], BSobj)

image.png

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