最好的文档其实还是官方文档。。。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"))
注释结果可直接用于各种分析~
