DSS包寻找DMS/DMR的统计学原理
参考语雀用户笔记：DSS检测差异甲基化区域
对于WGBS甲基化数据，做差异分析通常会用到DSS来完成，那么DSS是如何寻找差异甲基化位点的呢？在文献《A Bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data》中介绍了如何利用 Bayesian hierarchical model 来寻找差异甲基化位点。
以下内容主要来源于公众号：生信技能树，其中内容有稍作修改。
背景
哺乳动物基因组 CpG 位点通常集中在称为 CpG 岛（CpG island，CGI）的区域中，并且已知人基因启动子~ 60％含有 CpG 岛。CpG 岛上下游不超过 2000 个碱基对（2kb）的基因组区域称为 CpG“岛岸”（shores），其中 CpG shelves 指位于 CpG shores 上下游 2kb 以内的区域，open sea 指 CpG islands、CpG shores 和 CpG shelves 之外的其他区域。这 4 种情况形成了 CpG resort。CpG 位点的密度从 island 到 open sea 递减。
3 个技术
主要是 甲基化测序的 WGBS 和 RRBS，还有 芯片：
全基因组 DNA 甲基化测序（Whole Genome Bisulfite Sequencing，WGBS）是 DNA 甲基化研究的金标准，它通过 Bisulfite 处理和全基因组 DNA 测序结合的方式，对整个基因组上的甲基化情况进行分析，具有单碱基分辨率，可精确评估单个 C 碱基的甲基化水平，构建全基因组精细甲基化图谱。数据量非常大。
简化甲基化测序 (Reduced representation bisulfite sequencing, RRBS) 是一种准确、高效、经济的 DNA 甲基化研究方法，通过酶切 (Msp I) 富集启动子及 CpG 岛区域，并进行 Bisulfite 测序，同时实现 DNA 甲基化状态检测的高分辨率和测序数据的高利用率。作为一种高性价比的甲基化研究方法，简化甲基化测序在大规模临床样本的研究中具有广泛的应用前景。
Illumina 的 Infinium BeadChip 芯片，包括 HumanMethyation450（450K）和 MethylationEPIC（850K）。Infinium 芯片存在染料偏差、不同探针化学和位置效应的问题，已知这些问题会影响结果，必须在数据处理过程中进行校正。Infinium 450K 探针交叉反应和模糊比对到人类基因组中的多个位置影响了 485,000 个探测器中的约 140,000 个探针（29％），将可用探针的数量减少到约 345,000 个。这个问题在新发布 850K 仍然存在，其包括 > 90％的 450K 探针。
有文章比较这 3 个技术：Empirical comparison of reduced representation bisulfite sequencing and Infinium BeadChip reproducibility and coverage of DNA methylation in humans
主要是分析 差异甲基化区域（DMRs）与 DMR 相关差异表达基因
数据介绍
这里选择 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52140
提供每个样本的信号值矩阵下载
下载并且了解数据:
查看压缩包内容，如下：
head GSM1084240_A3R_d6_250.cpgs.txtCHR POS A3R_d6_250chr1 10525 21/23chr1 10542 21/23chr1 10609 1/23chr1 10617 0/24chr1 10620 0/24chr1 10631 0/24chr1 10633 0/24chr1 10636 0/24chr1 10638 0/24
DSS 包要求输入文件数据:每一行代表一个 CpG site, 格式如下:
- 第一列为染色体
- 第二列为位置
- 第三列为 total reads
- 第四列为甲基化的 reads
所以我们下载的数据需要进行拆分，然后导入到 R 里面才能被 DSS 包使用。
DSS 包介绍
主要是把上面项目的数据文件下载，然后导入到 R 里面，是有 DSS 包进行分析。
DSS (Dispersion Shrinkage for Sequencing data)，为基于高通量测序数据的差异分析而设计的 Bioconductor 包。主要应用于 BS-seq（亚硫酸氢盐测序）中计算不同组别间差异甲基化位点（DML）和差异甲基化区域（DMR）即 Call DML or DMR。
输入文件格式
DSS要求的输入文件格式如下：
| Chr | Pos | N | X |
|---|---|---|---|
| chr18 | 3014904 | 26 | 2 |
| chr18 | 3031032 | 33 | 12 |
| chr18 | 3031044 | 33 | 13 |
其中：
1. Chr代表染色体号，染色体编号需要有”chr”字符，如”chr1”，而不能是”1”；
2. Pos代表位点的基因组坐标/位置
3. N代表该位点一共被多少条reads覆盖
4. X代表有多少条reads显示该位点为甲基化；
注：这种格式可以从bismark_methylation_extractor命令的输出文件转换而来，一种是”.bismark.cov”，另一种是”.CpG_report.txt”。
如何进行转换：使用shell下的AWK，把列的顺序重新编排一下，没有的列可以通过简单的加法获得。
DSS 包的使用主要包括：
- 输入文件的准备
- 利用 DMLtest 函数检验所有的位点
- 利用 callDML 函数挑选统计学显著的位点
- 利用 callDMR 函数 Call DMR
- 利用 showOneDMR 函数对 DMRs 可视化
首先我们导入上面 GSE52140 数据集的文件：
library(data.table)library(stringr)library(tidyverse) allDat <- lapply(list.files('GSE52140_RAW/',pattern='*cpgs.txt.gz'),function(f){ # f="GSM1251242_H2R_d0.cpgs.txt.gz"; print(f); tmp=fread(file.path('GSE52140_RAW/',f)) chr=as.character(tmp$CHR) pos=as.character(tmp$POS) newTmp=separate(tmp,col =3,into = c("methy", "unmethy"), sep = "/") newTmp$all=as.numeric(newTmp$methy)+as.numeric(newTmp$unmethy) newTmp=as.data.frame(newTmp[,c(1,2,5,3)]) colnames(newTmp)=c('chr', 'pos' ,'N' ,'X') return(newTmp)})## 值得注意的是每个样本的位点数量不一致哦do.call(rbind,lapply(allDat,dim))tmp=do.call(cbind,lapply(allDat,head))sn=gsub('.cpgs.txt.gz','',list.files('GSE52140_RAW/',pattern='*cpgs.txt.gz'))sn=gsub('GSM.*?_','',sn)sn
也就是说把这 17 个文件读入了，样本名字是：
> sn [1] "A0R_d0_rep1" "A3R_d0_rep1" "A3R_d6_250" "A3R_d6_1000" [5] "A3R_d13_250" "A3R_d13_1000" "H0R_d0_rep1" "H3R_d0_rep1" [9] "H3R_d6_250" "H3R_d13_250" "A0R_d0_rep2" "A3R_d0_rep2" [13] "H0R_d0_rep2" "H3R_d0_rep2" "A1R_d0" "A2R_d0" [17] "H2R_d0"
然后我们用下面的 3 个例子来说明这个 DSS 包的用法，需要掌握上面样本的命名：
# lung cancer cell lines A549 (A) and HTB56 (H)# normal cell lines (0R) # a highly metastatic phenotype (3R)# 5-Azacytidine treatment at low concentrations (250 nM & 1000 nM) # for 6 days, additional 7 days in regular medium
单样本 VS 单样本
代码如下，重要就是构建对象和做统计检验
library(DSS)require(bsseq) if(T){ BSobj <- makeBSseqData(allDat[1:2], c("A0R", "A3R") )[1:1000,] BSobj save(BSobj,file = 'single-BSobj.Rdata') # There is no biological replicates in at least one condition. dmlTest <- DMLtest(BSobj, group1=c("A0R"), group2=c("A3R"),smoothing=TRUE) head(dmlTest) }
多样本的组 VS 另一个组
代码如下，重要就是构建对象和做统计检验，这里比较 “A0R_d0”和”A3R_d0” 组别的 2 个样本：
if(T){ sn[c(1,11,2,12)] BSobj <- makeBSseqData(allDat[c(1,11,2,12)], sn[c(1,11,2,12)] )[1:1000,] BSobj save(BSobj,file = 'group-BSobj.Rdata') dmlTest <- DMLtest(BSobj, group1=c("A0R_d0_rep1","A0R_d0_rep2"), group2=c("A3R_d0_rep1","A3R_d0_rep2"),smoothing=F) head(dmlTest) }
多种比较方式
代码如下，重要仍然是构建对象和做统计检验，但是需要构建属性矩阵，而且增加了 DMLfit 步骤。
sn sn[grep('rep',sn)] cellline=substring(sn[grep('rep',sn)],1,1) type=substring(sn[grep('rep',sn)],2,2) design=data.frame(cellline=cellline,type=type) design
得到的属性矩阵如下：
sn sn[grep('rep',sn)] cellline=substring(sn[grep('rep',sn)],1,1) type=substring(sn[grep('rep',sn)],2,2) design=data.frame(cellline=cellline,type=type) design # 构建对象特别耗时； BSobj <- makeBSseqData(allDat[c(grep('rep',sn))], sn[grep('rep',sn)]) BSobj save(BSobj,file = 'multi-BSobj.Rdata') load(file = 'multi-BSobj.Rdata') DMLfit=DMLfit.multiFactor(BSobj,design, formula = ~cellline+type+cellline:type) colnames(DMLfit$X) # 这里可以使用 ‘coef’, ‘term’, or ‘Contrast’我们仅仅是演示 coef dmlTest=DMLtest.multiFactor(DMLfit,coef=2) head(dmlTest)
值得注意的是DMLtest.multiFactor结果不需要callDML，只需要callDMR即可!
结果介绍
不管是哪种比较，最后都得到dmlTest变量走后面的流程，包括确定显著差异甲基化区域及基因，以及可视化展现，代码如下：
# 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 functionshowOneDMR(dmrs[1,], BSobj)
很明显，参数都是可以调整的，统计学显著性的阈值自己把握。
若有收获，就点个赞吧
0 人点赞
