软件17 —— DiffBind

一、 基本介绍

DiffBind是鉴定两个样本间差异结合位点的一个R包。主要用于peak数据集，包括对peaks的重叠和合并的处理，计算peaks重复间隔的测序reads数，并基于结合亲和力鉴定具有统计显著性的差异结合位点。

该R包采用了RNA-seq中差异基因表达的思路来进行peak的差异分析，和macs2的差异功能不同，DiffBind需要依赖已有的peak calling结果，将peak区域当做RNA-seq中的基因区域，然后对这些区域进行定量和差异分析，其核心的差异分析通过调用RNA-seq中常用的R包来实现。适用的统计模型有DESeq、DESeq2、edgeR。

详细内容可参考DiffBind的文档：
http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf

二、 使用方法

(1) 安装

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("DiffBind")

library(DiffBind)

(2) 准备输入文件

DiffBind需要提供样本比对的bam文件以及peak calling得到的peak区域结果文件。为了方便导入，DiffBind提供了一个接口，将导入文件的相关信息保存在一个文件中，该文件内容包含的表头信息有"SampleID"、 "Tissue"、 "Factor"、 "Condition" 、"Treatment"、"Replicate" 、"bamReads" 、"ControlID"、 "bamControl"、"Peaks"、 "PeakCaller"。在实际分析中，可能有很多列没有对应信息，直接空值即可。

SampleID：样本的标识符字符串。每个样本必须是唯一的。
Tissue：Identifier string for tissue type.
Factor：Identifier string for factor.
Condition：Identifier string for condition.
Treatment：Identifier string for treatment.
Replicate：第几次重复。
bamReads：包含ChIP样本aligned reads的bam文件的文件路径。
ControlID：Call peak时使用的input数据的ID.
bamControl：包含control样本aligned reads的bam文件的文件路径。
Spikein：包含spike-in样本aligned reads的bam文件的文件路径。
Peaks：包含样品峰的文件路径。由PeakCaller字段或caller参数确定的格式。这里有多种数据格式可作为输入：1. macs2。输出的.narrowPeak等峰文件。2. 包括所有call peak 得到的peak位置信息的.bed 文件。3. 以上两种格式得到的.gz文件。
PeakCaller：使用的peak caller的标识符字符串。如果Peaks不是bed文件，这将决定如何解析Peaks文件。如果缺少，将使用peak caller参数中指定的默认peak caller。可能的值为：

“raw”：text file file; peak score is in fourth column
“bed”：.bed file; peak score is in fifth column
“narrow”：default peak.format：narrowPeaks file
“macs”：MACS .xls file
“swembl”：SWEMBL .peaks file
“bayes”：bayesPeak file
“peakset”：peakset written out using pv.writepeakset
“fp4”：FindPeaks v4

PeakFormat：指示峰值文件格式的字符串。
ScoreCol：峰文件中包含峰值分数的列。
LowerBetter：逻辑上表示较低的分数意味着较好的峰值。
Counts：外部计算的读取计数的文件路径。

(3) 进行分析

Diffbind进行了高度封装，所有的函数都围绕一个自定义的DBA对象为中心，整个过程分为以下4步：

• count，计算peak区域的表达量，由于不同的peak数据集会存在overlap，所以首先合并peak区域，当导入的peak数据集越多，理论上合并后的peak平均宽度就会越宽，overlap的peak越多，合并后的peak机会越宽。正是由于merge机制的存在，最终定量结果中的peak无论是个数还是宽度都和输入的不太一致。
• contrast，构建比较的分组，指定哪些分组进行比较。DiffBind要求必须有生物学重复，每组至少有两个样本，否则会报错。
• analyze，根据定量结果，调用DESeq等R包进行差异分析。
• report，提取差异分析结果。

(4) 举个例子

# 读取数据
dbObj <- dba(sampleSheet="13_DiffBind/SampleSheet.csv")

> dbObj
10 Samples, 471 sites in matrix (760 total):
    ID   Tissue   Condition   Replicate   Intervals
1   WTF4   larvae    WTF      1       381
2   WTF5   larvae    WTF      2       490
# 471表示至少两个样本中有重叠的峰的数量，760表示所有的重叠峰放在一起并去重后的数量。

# 质控
dba.plotPCA(dbObj, attributes=DBA_CONDITION, label=DBA_ID)
plot(dbObj)

# 计数
dbObj <- dba.count(dbObj, bUseSummarizeOverlaps=TRUE) 

> dbObj
10 Samples, 450 sites in matrix:
    ID   Tissue   Condition   Replicate   Reads    FRiP
1   WTF4   larvae    WTF     1   2252267   0.06
2   WTF5   larvae    WTF     2   4769374   0.04
# dbObj中多了两列，Reads表示每个样品中所有比对上的reads数量；FRiP表示Fraction of Reads in Peaks，即该样本peaks上的reads占所有reads的百分比，表示样本富集的效果。

# 基于测序深度标准化
dbObj <- dba.normalize(dbObj)

# 构建比较的分组 Establishing a contrast 
dbObj <- dba.contrast(dbObj, categories=DBA_CONDITION, minMembers=2)

# 进行差异分析
dbObj <- dba.analyze(dbObj, method=DBA_DESEQ2)
#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)

# 保存文件
write.csv(as.data.frame(comp1.deseq), file="13_DiffBind/WTF_vs_WTM.csv", row.names = F)

seqnames   start   end   width   strand
2L   2755822   2756222   401   *
2L   14743325   14743725   401   *
Conc   Conc_WTF   Conc_WTM   Fold   p.value   FDR
7.009378342   0   8.009378342   -20.27812504   2.60E-16   9.93E-14
7.675999254   0   8.675999254   -10.85136276   2.95E-13   5.64E-11
# group2被设置为了对照
# Conc指的是Mean read concentration over all the samples (进行了log2的处理)
# Conc_WTF指的是在所有WTF的样本中的统计
# Fold这里指的是WTF/WTM ，为正数则表示 Increased binding affinity in the WTF group.

# 以bed格式保存显著性的差异结果
# Create bed files for each keeping only significant peaks (p < 0.05)
out <- as.data.frame(comp1.deseq)
deseq.bed <- out[which(out$p.value < 0.05), 
                 c("seqnames", "start", "end", "strand", "Fold")]
write.table(deseq.bed, file="13_DiffBind/WTF_vs_WTM.bed", sep="\t", 
            quote=F, row.names=F, col.names=F)

# 火山图
dba.plotVolcano(dbObj, contrast = 2, bUsePval = T, th = 0.05)

# 热图
readscores <- dba.plotHeatmap(dbObj, contrast=1, bUsePval = T, th = 0.05, 
                              correlations=FALSE, scale="row", 
                              colScheme=colorRampPalette(c("red", "black", "green"))(n = 13))
readscores <- as.data.frame(readscores)