第9篇:差异peaks分析——DiffBind
ATAC系列连载:
第1篇:ATAC-seq的背景介绍以及与ChIP-Seq的异同
第4篇:对ATAC-Seq/ChIP-seq的质量评估(一)——phantompeakqualtools
第5篇:对ATAC-Seq/ChIP-seq的质量评估(二)——ChIPQC
第7篇:用Y叔的ChIPseeker对peaks进行注释和可视化
学习目标
学习用DiffBind流程评估两个样本间的差异结合区域
用PCA评估样本间的关系
评估不同工具得到的差异peaks的一致性
评估差异peaks富集的工具
ATAC-Seq下游分析的另一个重点是差异peaks分析。如分析不同的实验条件、多个时间节点、不同的发育时期等的差异区域。鉴定这些差异peaks区域在生物医学研究中也具有重要意义,目前也有多种相关的工具被开发:
选择合适的工具需考虑以下几个因素:
所用的软件不需要大量的代码移植,该工具是否被维护和频繁更新。
用户需要提供什么样的输入文件,是否与peaks calling所用工具可以衔接。
信号分布的基础统计模型是什么?是基于泊松分布或更灵活的负二项分布。
一些工具是专门针对特定的ATAC-seq或ChIP-seq 数据(信号类型)设计的,如组蛋白修饰或转录因子(TF)结合。
具体选择哪种工具决定于实验设计,下面的决策树可以帮助缩小你的选择范围。
DiffBind
DiffBind是鉴定两个样本间差异结合位点的一个R包。主要用于peak数据集,包括对peaks的重叠和合并的处理,计算peaks重复间隔的测序reads数,并基于结合亲和力鉴定具有统计显著性的差异结合位点。适用的统计模型有DESeq、DESeq2、edgeR。详细内容可参考DiffBind的文档:http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf
使用方法:
1. 下载安装DiffBind
source("https://bioconductor.org/biocLite.R")
biocLite("DiffBind")
# view documentation
browseVignettes("DiffBind")
2. 准备输入文件
需要准备一个SampleSheet文件,与ChIPQC的方法一样。SampleSheet文件是根据实验设计和数据存储地址等信息创建的一个csv格式文件,包含的表头信息有"SampleID"、 "Tissue"、 "Factor"、 "Condition" 、"Treatment"、"Replicate" 、"bamReads" 、"ControlID"、 "bamControl" "Peaks"、 "PeakCaller"(bam,peak文件分别在比对和call peak的步骤产生)。
3. 差异peaks分析
读入文件
一旦读入了peaksets,合并函数就找到所有重叠的peaks,并导出一致性的peaksets。
library(DiffBind)
dbObj <- dba(sampleSheet="SampleSheet.csv")
亲和结合矩阵
计算每个peaks/regions的count信息。先对一致性的peaks数据集进行标准化,然后根据他们的峰值(point of greatest read overlap)再次中心化并修剪一些peaks,最终得到更加标准的peak间隔。使用函数dba.count(),参数
bUseSummarizeOverlaps
可以得到更加标准的计算流程。dbObj <- dba.count(dbObj, bUseSummarizeOverlaps=TRUE)
dba.plotPCA(dbObj, attributes=DBA_FACTOR, label=DBA_ID)
plot(dbObj)
结果中同时计算了FRiP (质量评估的一个标准,可以参考第5篇:对ATAC-Seq/ChIP-seq的质量评估(二)——ChIPQC)。
样本间的聚类:
Venn图展示样本间peaks的重合
差异分析
差异分析是DiffBind的核心功能,默认情况是基于DEseq2, 可以设置参数
method=DBA_EDGER
选择edgeR,或者设置method=DBA_ALL_METHODS
。每种方法都会评估差异结果的p-vaue和FDR。# Establishing a contrast
dbObj <- dba.contrast(dbObj, categories=DBA_FACTOR,minMembers = 2)
dbObj <- dba.analyze(dbObj, method=DBA_ALL_METHODS)
# summary of results
dba.show(dbObj, bContrasts=T)
# overlapping peaks identified by the two different tools (DESeq2 and edgeR)
dba.plotVenn(dbObj,contrast=1,method=DBA_ALL_METHODS)
提取结果
comp1.edgeR <- dba.report(dbObj, method=DBA_EDGER, contrast = 1, th=1)
comp1.deseq <- dba.report(dbObj, method=DBA_DESEQ2, contrast = 1, th=1)
结果文件包含所有位点的基因组坐标,以及差异富集的统计数据包括fold change、p值和FDR。
保存文件
out <- as.data.frame(comp1.edgeR)
write.table(out, file="results/Nanog_vs_Pou5f1_edgeR.txt", sep="\t", quote=F, col.names = NA)
out <- as.data.frame(comp1.deseq)
write.table(out, file="results/Nanog_vs_Pou5f1_deseq2.txt", sep="\t", quote=F, col.names = NA)
以bed格式保存显著性的差异结果
out <- as.data.frame(comp1.edgeR)
edge.bed <- out[ which(out$FDR < 0.05),
c("seqnames", "start", "end", "strand", "Fold")]
write.table(edge.bed, file="results/Nanog_vs_Pou5f1_edgeR_sig.bed", sep="\t", quote=F, row.names=F, col.names=F)
out <- as.data.frame(comp1.deseq)
deseq.bed <- out[ which(out$FDR < 0.05),
c("seqnames", "start", "end", "strand", "Fold")]
write.table(deseq.bed, file="results/Nanog_vs_Pou5f1_deseq2_sig.bed", sep="\t", quote=F, row.names=F, col.names=F)
参考资料
DiffBind文档
哈佛深度NGS数据分析课程, 09- Differential Peak calling using DiffBind
独家福利
如果需要组装自己的服务器;代办生物信息学服务器
如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等),点我?
如果需要线下辅导及培训,看招学徒
如果需要个人电脑:个人计算机推荐
如果需要置办生物信息学书籍,看:生信人必备书单
如果需要实习岗位:实习职位发布
如果需要售后:点我