使用DSS包多种方式检验差异甲基化信号区域

Posted on by

一个背景

哺乳动物基因组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）

​

Comments are closed.