差異peaks分析

1.學習目標

學習用DiffBind流程評估兩個樣本間的差異結合區(qū)域
用PCA評估樣本間的關系
評估不同工具得到的差異peaks的一致性

2.DiffBind介紹

DiffBind是鑒定兩個樣本間差異結合位點的一個R包溪掀。主要用于peak數據集，包括對peaks的重疊和合并的處理，計算peaks重復間隔的測序reads數，并基于結合親和力鑒定具有統(tǒng)計顯著性的差異結合位點绰上。適用的統(tǒng)計模型有DESeq、DESeq2纯衍、edgeR竖共。詳細內容可參考DiffBind的文檔：http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf。
該包的主要重點是確定樣本之間有差異的位點乳绕。它包括支持峰集處理的功能绞惦，包括重疊和合并峰集，在峰集里進行重疊區(qū)間的測序read計數洋措，以及基于結合親和的證據(通過read密度的差異測定)識別統(tǒng)計上顯著的差異結合位點济蝉。為此，它使用了在RNA-Seq環(huán)境中開發(fā)的統(tǒng)計包(edgeR和DESeq2)菠发。此外王滤，這個包構建在Rgraphics的基礎上，提供了一組標準化的圖來幫助進行結合分析滓鸠。
DiffBind主要對峰集(peaksets)進行分析雁乡，峰集是一組代表候選蛋白質結合位點的基因組區(qū)間（也適用于ATAC-seq的峰集）。每個區(qū)間包括一個染色體糜俗，一個開始和結束位置踱稍，通常還有一個表示對峰得confidence或強度的分數曲饱。與每個峰集相關聯(lián)的是與產生峰集的實驗相關的metadata。此外珠月，包含比對上的測序read文件(.bam文件)可以與每個峰集關聯(lián)扩淀。
通常的，DiffBind處理數據一般是五個部分：
讀取峰集：第一步是讀取一組峰集和相關的metadata桥温。峰集要么來自ChIP-Seq的peak callers引矩，如MACS梁丘，要么使用其他標準(如基因組窗口侵浸，或基因組中
的所有啟動子區(qū)域)。讀取峰集的最簡單方法是使用逗號分隔值(csv)列表氛谜，每個峰集一行(也接受帶有.xls或.xlsx后綴的Excel格式的電子表格)掏觉。單個樣本可以有多個相關峰集。例如值漫，如果使用多個peak callers進行比較澳腹，每個樣本將在列表中多出一行。
Occupancy分析：峰集杨何，特別是由peak callers產生的酱塔，提供了對特定基因組位點上某個蛋白質的潛在 occupancy的結果。peaksets已加載后危虱，可以執(zhí)行一些探索性繪圖羊娃，確定這些occupancy圖譜在實驗重復之間的一致性。如實驗之間的重復埃跷，在同一實驗中使用不同peak callers蕊玷。DiffBind提供了能夠檢查重疊的功能，以及確定相似樣本聚集在一起的函數弥雹。此外垃帅，峰可能基于blacklists過濾，以及自定義greylists來源于control track特定的實驗(見第6節(jié))剪勿。除了質量控制贸诚，Occupancy的分析可以是一個共識峰集，代表一組整體的候選結合位點用于進一步分析厕吉。
reads計數：一旦一個共識的峰集已經被推導出來赦颇，DiffBind可以使用提供的測序read文件來計算每個樣本的每個區(qū)間有多少reads重疊。默認情況下赴涵，為了提供更多標準化的峰值區(qū)間媒怯，共識峰集中的峰會根據其峰值(最大讀重疊點)重新調整中心點和trimmed。計數的最終結果是一個結合親和矩陣髓窜，其中包含每個樣本在每個共識結合位點的read count扇苞，無論它是否被識別為該樣本中的一個峰欺殿。有了這個矩陣，樣本可以使用affinity重新聚類鳖敷，而不是occupancy數據脖苏。結合親和矩陣用于QC繪圖以及隨后的差異分析。
差異結合親和分析：DiffBind的核心功能是差異結合親和分析定踱，它可以識別樣本間顯著的差異結合位點棍潘。這一步驟包括將實驗數據標準化，建立模型設計和對比(或contrasts)崖媚。接下來執(zhí)行底層的核心分析亦歉，默認情況下使用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的别瞭，所以這里factor不是很重要窿祥，如果用ChIP-seq的數據，這里就是你的特定蛋白名稱蝙寨∩柜茫“Treatment”就是分組嗤瞎，對照或者不同處理，也可以是對照和過表達/KO等 等听系”雌妫“bamReads”是bam文件的絕對路徑】渴ぃ“Peaks”是call peak之后得到的peak文件的文件夾掉瞳。

atdiff.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圖

相關熱圖

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)

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)

plotMA

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

   結果文件包含所有位點的基因組坐標瞳购，以及差異富集的統(tǒng)計數據包括fold change、p值和FDR亏推。其中Conc的表示read的平均濃度学赛，即對peak 的read counts進行l(wèi)og2標準化

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)

PCA

火山圖

#把某兩個組所有的差異峰用熱圖表示出來：
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)

所有樣品熱圖

#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)

image.png