R语言 -- 寻找差异甲基化区域（DMR）-- DSS 包

最好的文档其实还是官方文档。。。http://bioconductor.org/packages/release/bioc/vignettes/DSS/inst/doc/DSS.html#3_Using_DSS_for_BS-seq_differential_methylation_analysis

有一个步骤简直太慢了，所以等的时候顺便简单记录一下~

适用于各种甲基化测序（WGBS、RRBS、TBS、芯片等）

1. Data preparation

上游分析略（请参考)，我是用的是 BSMAP + MethylDackel extract 流程，最终得到甲基化bedGraph文件：

head "/sibcb2/bioinformatics2/wangjiahao/GDM/RRBS/BS_workflow/M17A020338-M_val/mCall/M17A020338-M_val_CpG.bedGraph"
track type="bedGraph" description="BS_workflow/M17A020338-M_val/mCall/M17A020338-M_val CpG methylation levels"
chr1    10496   10497   94      16      1
chr1    10524   10525   88      15      2
chr1    10525   10526   86      13      2
chr1    10541   10542   100     17      0
chr1    10542   10543   93      14      1
chr1    10562   10563   94      16      1
chr1    10563   10564   86      13      2
chr1    10570   10571   94      16      1
chr1    10571   10572   86      13      2

第五列：甲基化reads数，第六列：未甲基化reads数。

2. Start DSS ！

• 收集数据

# 自定义的函数
read.faster <- function(file, header = FALSE, sep = "\t", showProgress = TRUE, skip = 0){

    suppressPackageStartupMessages(library("data.table"))
    Tx = data.frame(fread(file = file, header = header, sep = sep, showProgress = showProgress, skip = skip))
    return(Tx)
}

read.bedgraph <- function(file, showProgress = FALSE){

    Tx = read.faster(file = file, header = FALSE, sep = "\t", skip = 1, showProgress = showProgress)
    return(Tx)
}

library("DSS")
library("bsseq") 
SRX = c(Normal_B, GDM_B)
files = list.files(sourceFolder, "CpG.bedGraph", recursive = TRUE, full.names = TRUE)
files = sapply(SRX, function(x) grep(x, files, value = TRUE))
allDat = lapply(files, function(file){
    cat("Remian", length(files) - match(file, files), "\n")
    T1 = read.bedgraph(file)
    T2 = data.frame(chr = T1[,1], pos = T1[,3], N = T1[,5] + T1[,6], X = T1[,5]) # 列名必须为这四个
    return(T2)
})

sourceFolder ：含有所有样本CpG.bedGraph文件的文件夹

Normal_B、GDM_B： 两组样本的 ID

因为我的分组有四组，因此还要把这次需要分析的两组给挑出来

下一步的构建对象需要提供包含所有样本甲基化信息的列表，因此使用lapply可方便获得

每个样本的甲基化数据格式为:

> hd(allDat[[1]])
   chr   pos N X
1 chr1 10497        12         12
2 chr1 10525        12         12
3 chr1 10526         6          6
4 chr1 10542        12         11
5 chr1 10543         6          6

N 为 总reads数，X 为甲基化的reads数，因此需要转化一下

• 构建对象、callDML、callDMR

BSobj   <- makeBSseqData(allDat, SRX) #; 时间较长长长长。。。（1长约等于5分钟，后面已放弃计时，41个RRBS样本我可能花了一个半小时，可能服务器最近很渣？
dmlTest <- DMLtest(BSobj, group1 = Normal_B, group2 = GDM_B, smoothing = FALSE)
dmls    <- callDML(dmlTest, delta = 0.1, p.threshold = 0.05)
dmrs    <- callDMR(dmlTest, delta = 0.1, p.threshold = 0.05)

WGBS数据需要smoothing = TRUE，RRBS等测序方式官方建议也最好加上，虽然对结果影响不大：

User has the option to smooth the methylation levels or not. A simple moving average algorithm is implemented for smoothing. For WGBS data, smoothing is always recommended so that information from nearby CpG sites can be combined to improve the estimation of methylation levels. For data with sparse CpG coverage such as RRBS or hydroxymethylation, smoothing might not alter the results much, but is still recommended. In RRBS, CpG sites are likely to be clustered locally within small genomic regions, so smoothing can potentially help the methylation estimation.

If smoothing is requested, smoothing span is an important parameter which has non-trivial impact on DMR calling. We use 500 bp as default, because it performs well in real data tests according to our experience.

如果提供的数据列名不是chr pos N X则会报以下错误：

 Error in `[.data.frame`(alldat, , "N") : undefined columns selected

• 注释DMRs

write.bed <- function(Tx, file){

    write.table(Tx, file = file, sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE)
}

library("ChIPseeker")
write.bed(dmrs[1:3], file="txt/DMRs_placenta.bed")
dmrs_gr   <- readPeakFile("txt/DMRs_placenta.bed")
peak_anno <- data.frame(annotatePeak(dmrs_gr, annoDb = "org.Hs.eg.db"))

注释结果可直接用于各种分析~

大力感谢：使用DSS包多种方式检验差异甲基化信号区域