使用DSS包多种方式检验差异甲基化信号区域
因为一个有趣的原因，我们来一起学一个R包吧
一个背景
哺乳动物基因组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.txt
CHR POS A3R_d6_250
chr1 10525 21/23
chr1 10542 21/23
chr1 10609 1/23
chr1 10617 0/24
chr1 10620 0/24
chr1 10631 0/24
chr1 10633 0/24
chr1 10636 0/24
chr1 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包的使用主要包括：
输入文件的准备
利用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 function
showOneDMR(dmrs[1,], BSobj)
很明显，参数都是可以调整的，统计学显著性的阈值自己把握。
作业
分析文章 Enduring epigenetic landmarks define the cancer microenvironment ，拿到 患者间差异甲基化区域（DMRs）与 DMR 相关差异表达基因（DE-DMRs）