差异peaks分析
1.学习目标
学习用DiffBind流程评估两个样本间的差异结合区域
用PCA评估样本间的关系
评估不同工具得到的差异peaks的一致性
2.DiffBind介绍
DiffBind是鉴定两个样本间差异结合位点的一个R包。主要用于peak数据集，包括对peaks的重叠和合并的处理，计算peaks重复间隔的测序reads数，并基于结合亲和力鉴定具有统计显著性的差异结合位点。适用的统计模型有DESeq、DESeq2、edgeR。详细内容可参考DiffBind的文档：http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf。
该包的主要重点是确定样本之间有差异的位点。它包括支持峰集处理的功能，包括重叠和合并峰集，在峰集里进行重叠区间的测序read计数，以及基于结合亲和的证据(通过read密度的差异测定)识别统计上显著的差异结合位点。为此，它使用了在RNA-Seq环境中开发的统计包(edgeR和DESeq2)。此外，这个包构建在Rgraphics的基础上，提供了一组标准化的图来帮助进行结合分析。
DiffBind主要对峰集(peaksets)进行分析，峰集是一组代表候选蛋白质结合位点的基因组区间（也适用于ATAC-seq的峰集）。每个区间包括一个染色体，一个开始和结束位置，通常还有一个表示对峰得confidence或强度的分数。与每个峰集相关联的是与产生峰集的实验相关的metadata。此外，包含比对上的测序read文件(.bam文件)可以与每个峰集关联。
通常的，DiffBind处理数据一般是五个部分：
读取峰集：第一步是读取一组峰集和相关的metadata。峰集要么来自ChIP-Seq的peak callers，如MACS，要么使用其他标准(如基因组窗口，或基因组中
的所有启动子区域)。读取峰集的最简单方法是使用逗号分隔值(csv)列表，每个峰集一行(也接受带有.xls或.xlsx后缀的Excel格式的电子表格)。单个样本可以有多个相关峰集。例如，如果使用多个peak callers进行比较，每个样本将在列表中多出一行。
Occupancy分析：峰集，特别是由peak callers产生的，提供了对特定基因组位点上某个蛋白质的潜在 occupancy的结果。peaksets已加载后，可以执行一些探索性绘图，确定这些occupancy图谱在实验重复之间的一致性。如实验之间的重复，在同一实验中使用不同peak callers。DiffBind提供了能够检查重叠的功能，以及确定相似样本聚集在一起的函数。此外，峰可能基于blacklists过滤，以及自定义greylists来源于control track特定的实验(见第6节)。除了质量控制，Occupancy的分析可以是一个共识峰集，代表一组整体的候选结合位点用于进一步分析。
reads计数：一旦一个共识的峰集已经被推导出来，DiffBind可以使用提供的测序read文件来计算每个样本的每个区间有多少reads重叠。默认情况下，为了提供更多标准化的峰值区间，共识峰集中的峰会根据其峰值(最大读重叠点)重新调整中心点和trimmed。计数的最终结果是一个结合亲和矩阵，其中包含每个样本在每个共识结合位点的read count，无论它是否被识别为该样本中的一个峰。有了这个矩阵，样本可以使用affinity重新聚类，而不是occupancy数据。结合亲和矩阵用于QC绘图以及随后的差异分析。
差异结合亲和分析：DiffBind的核心功能是差异结合亲和分析，它可以识别样本间显著的差异结合位点。这一步骤包括将实验数据标准化，建立模型设计和对比(或contrasts)。接下来执行底层的核心分析，默认情况下使用DESeq2。这将为每个候选结合位点分配一个p值和FDR，表明它们的差异结合置信度。
画图和reporting：一旦运行了一个或多个contrasts，DiffBind就提供了许多用于报告和绘制结果的函数。MA和火山图给出了分析结果的概述，而相关性热图和PCA图显示了如何基于不同的结合位点的组聚类。箱形图显示了reads在差异结合位点内的分布，对应于它们在两个样本组之间是否获得或失去亲和力。reporting机制能够提取差异结合位点进行进一步处理，如注释、motif和通路分析。
数据分析
1.DiffBind的安装
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DiffBind")
library(DiffBind)
2.准备sample sheet
sample sheet是一个列表，需要包括以下几列："SampleID"， "Tissue"， "Factor"， "Condition" ，"Treatment"，"Replicate"， "bamReads" ，"ControlID" ，"bamControl" ，"Peaks"和"PeakCaller"。
这个列表可以是dataframe，或者是从csv文件读取，比如：
NOTE: 我用的数据是ATAC-seq数据，如果是ChIP-seq数据的话上面这个表里需要填写的内容会有些不同。具体请参考官网：
https://www.bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf
3.读取峰集（peaksets）
DBA_object <- dba(sampleSheet = samples)
DBA_object #这里得到的就是上面提到的共识峰
###我用的是
#导入数据，dba函数会将会将bam文件一同导入，因此csv中的路径一定要准确
tamoxifen <- dba(sampleSheet="atDiff.csv")
#tamoxifen <- dba.count(tamoxifen, minOverlap=3, bParallel=FALSE)
tamoxifen <- dba.count(tamoxifen, minOverlap=3, bParallel=FALSE, bUseSummarizeOverlaps=TRUE)
NOTE: 这一步有很多可选参数，我就直接用的都是默认值了，如果有特殊需要的，请参考DiffBind官网说明书。
使用peak calls的数据，可以生成一个相关热图，通过结合矩阵的每一行的相互关系给出样本的初始聚类:
plot(DBA_object)
###我用的是
plot(tamoxifen) ###查看样本间的相关性，得到相关热图。
dba.plotPCA(tamoxifen, attributes=DBA_FACTOR, label=DBA_ID) ###pca图
4.计算reads
下一步是根据每个样本的read计数(亲和力分数)计算一个带有分数的结合矩阵，而不是仅针对特定样本中call的那些峰(occupancy分数)计算置信分数。这些reads是使用 dba.count 功能获得的:
DBA_object <- dba.count(DBA_object,bUseSummarizeOverlaps=TRUE)
DBA_object
#计算每个样本中的reads数
tamoxifen <- dba.count(tamoxifen)
#简单查看计算结果
info <- dba.show(tamoxifen)
libsizes <- cbind(LibReads=info$Reads, FRiP=info$FRiP, PeakReads=round(info$Reads * info$FRiP))
rownames(libsizes) <- info$ID
libsizes
write.csv(libsizes, file="libsizes.csv")
###我用的是
tamoxifen <- dba.count(tamoxifen, minOverlap=3, bParallel=FALSE, bUseSummarizeOverlaps=TRUE)
这里同样有很多参数可以选择，这里我只选择一个参数：bUseSummarizeOverlaps 。
这个参数会使得运行比较缓慢，但是是一个更标准的计算功能。如果你把它设置为TRUE，所有的read文件必须是bam，并且必须有其自己的索引文件（.bam.bai）。另外fragmentSize参数必须是缺省值。
结果如下：
这里添加了两个新列。第一个显示了每个样本的比对的reads的总数(“完整”文库大小)。第二个是标记的FRiP，它代表在峰的Fraction of Reads。这是共识峰集中与这个样品里某一重叠的峰所占reads的比例，可以用来表示总体上哪个样品富集程度更高。
这时可以画个PCA plot:
这时候的相关性热图就不一样了：
plot(DBA_object)
5.构建design和contrast模型
DBA_object <- dba.contrast(DBA_object,categories=DBA_TREATMENT,minMembers = 2)
6.差异分析
差异分析是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)
DBA_object <- dba.analyze(DBA_object, method=DBA_DESEQ2)
dba.show(DBA_object, bContrasts=T) #这里我有4组不同的处理，那么两两比较就有6种组合，这里还要注意的是Group2是baseline。
Group1 Members1 Group2 Members2 DB.DESeq2
1 two 2 one 2 16216
2 two 2 four 2 21686
3 two 2 three 2 25708
4 one 2 four 2 2827
5 one 2 three 2 1382
6 four 2 three 2 1231
plot(DBA_object, contrast=6) #这里的contrast就是上面6行里最开始的一列“1-6”的数字，这里设置6，是针对第六组组合比较：我的“four”处理组合“three”处理组的相关性热图
还可以画个MA-plot看一下：
dba.plotMA(DBA_object,contrast = 6)
7.提取差异分析结果
这一部分你可以设置提取哪些组比较的差异结果，如果你只有两个组比较，就很容易，用默认值就好了。如果你的数据和我一样，我是4个组比较，上面提到就有6组比较组合，这时我可以设置提取哪两个实验条件的组合结果：
comp1.edgeR <- dba.report(dbObj, method=DBA_EDGER, contrast = 1, th=1)
comp1.deseq <- dba.report(dbObj, method=DBA_DESEQ2, contrast = 1, th=1)
> DBA_object.DB <- dba.report(DBA_object,contrast = 6)
> DBA_object.DB
GRanges object with 1231 ranges and 6 metadata columns:
seqnames ranges strand | Conc Conc_four Conc_three Fold p-value FDR
<Rle> <IRanges> <Rle> | <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
302 chr1 24612300-24613350 * | 13.35 13.02 13.62 -0.60 2.86e-13 1.11e-08
6232 chr11 57091050-57091850 * | 5.91 2.05 6.86 -4.81 5.55e-13 1.11e-08
15938 chr16 67435750-67437100 * | 7.36 8.02 6.11 1.91 1.15e-09 9.66e-06
10844 chr13 73455700-73456600 * | 5.74 6.60 3.31 3.29 1.15e-09 9.66e-06
25740 chr3 151351500-151352350 * | 7.48 6.42 8.08 -1.66 1.21e-09 9.66e-06
... ... ... ... . ... ... ... ... ... ...
12086 chr14 20765350-20766300 * | 6.76 6.13 7.20 -1.08 0.00153 0.0497
5266 chr10 116739900-116740500 * | 5.19 4.11 5.80 -1.69 0.00153 0.0497
29648 chr5 129974800-129975450 * | 5.20 3.79 5.90 -2.10 0.00154 0.0498
30836 chr6 47920250-47921000 * | 5.92 5.10 6.44 -1.34 0.00154 0.0498
2950 chr1 183905150-183905850 * | 6.41 5.70 6.88 -1.18 0.00154 0.0499
-------
seqinfo: 20 sequences from an unspecified genome; no seqlengths
结果文件包含所有位点的基因组坐标，以及差异富集的统计数据包括fold change、p值和FDR。其中Conc的表示read的平均浓度，即对peak 的read counts进行log2标准化
The value columns show the mean read concentration over all the samples (the default calculation uses log2 normalized ChIP read counts with control read counts subtracted) and the mean concentration over the first (Resistant) group and second (Responsive) group.
保存文件
# EdgeR
out <- as.data.frame(comp1.edgeR)
write.table(out, file="results/Nanog_vs_Pou5f1_edgeR.txt", sep="t", quote=F, col.names = NA)
# DESeq2
out <- as.data.frame(comp1.deseq)
write.table(out, file="results/Nanog_vs_Pou5f1_deseq2.txt", sep="t", quote=F, col.names = NA)
以bed格式保存显著性的差异结果
# Create bed files for each keeping only significant peaks (p < 0.05)
# EdgeR
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)
# DESeq2
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)
#看一下升高和降低的差异峰的数量
> sum(DBA_object.DB$Fold>0)
[1] 155
> sum(DBA_object.DB$Fold<0)
[1] 1076
#这时可以针对four组和three组的样品画PCA plot：
> dba.plotPCA(DBA_object, contrast= 6, attributes=DBA_TREATMENT, label=DBA_ID)
#火山图
dba.plotVolcano(DBA_object,contrast = 6)
#把某两个组所有的差异峰用热图表示出来：
hmap <- colorRampPalette(c("red", "black", "green"))(n = 13)
readscores <- dba.plotHeatmap(DBA_object, contrast=6, correlations=FALSE,scale="row", colScheme = hmap)
#当然你也可以把所有的样品都画出来：
> readscores <- dba.plotHeatmap(DBA_object,correlations=FALSE,scale="row",colScheme = hmap)
但是这里有一个问题，它把我的同一个Treatment的样品分隔开了。我尝试过调整列的顺序，但是貌似在 dba.plotHeatmap 这个函数里没有类似的参数可以调整。于是我试着用pheatmap将这个热图进行可视化。
使用pheatmap进行调整：
#reorder sample sequences(using pheatmap)
> reorder_matrix = readscores@elementMetadata[c(3,2,1,4,5,6,7,8)]
> color <- colorRampPalette(c("green","black", "red"))(n = 50)
> annotation_col = data.frame(group = c("one","one","three","three","two","two","four","four"), replicates = c("1","2"))
> rownames(annotation_col) = colnames(reorder_matrix)
# set the range of scale legend bar
> n=t(scale(t(reorder_matrix)))
> n[n> 1.5]=1.5 #限定上限
> n[n< -1.5]= -1.5 #限定下限
# plot
> pheatmap(n,
color = color,
cluster_cols=FALSE,
annotation_col = annotation_col)
