ChIPQC——对ChIP-seq的质量评估
第4篇：對ChIP-seq的質量評估
學習目標：
文章目錄
- 第4篇：對ChIP-seq的質量評估
- 一、ChIP-Seq質量評估
- 二、鏈交叉相關（Strand cross-correlation）
- 1. phantompeakqualtools
- 1.1 軟件安裝
- 1.2 使用phantompeakqualtools計算交叉相關性和其他相關的質控度量值
- 三、ChIPQC對數據進行質量評估
- 1. 實際操作
- 1.1 設置
- 1.2 軟件安裝
- 1.3 Running ChIPQC
- 1.4 運行 ChIPQC時報錯及解決
- 2. ChIPQC: quality metrics report
- 2.1 ChIPQC report
- 2.1.1 Metrics: Read characteristics
- 2.1.2 Metrics: Enrichment of reads in peaks
- 2.1.3 Metrics: Peak profiles
- 2.1.4 Metrics: Peak signal strength
- 2.2 Final takehome from ChIPQC
- 2.3 低質量數據的可能來源是什么?
一、ChIP-Seq質量評估
在下游分析前，最好是先對peak calling 后的ChIP-Seq數據進行質量評估。
二、鏈交叉相關（Strand cross-correlation）
1. phantompeakqualtools
1.1 軟件安裝
安裝方法一：https://anaconda.org/bioconda/phantompeakqualtools
安裝方法二：https://github.com/kundajelab/phantompeakqualtools
這里采用方法一進行安裝：
conda install -c bioconda phantompeakqualtools1.2 使用phantompeakqualtools計算交叉相關性和其他相關的質控度量值
參考：第4篇：對ATAC-Seq/ChIP-seq的質量評估（一）——phantompeakqualtools
三、ChIPQC對數據進行質量評估
學習目標
參考一：https://github.com/hbctraining/Intro-to-ChIPseq/blob/master/lessons/06_combine_chipQC_and_metrics.md
參考二：https://www.jianshu.com/p/a11147808d14
Prior to performing any downstream analyses with the results from a peak caller, it is best practice to assess the quality of your ChIP-seq data. What we are looking for is good quality ChIP enrichment over background.
1. 實際操作
1.1 設置
- mkdir一個新文件夾進行操作，將需要的文件轉移到新文件夾中；需要的文件有：
- BAM files
首先對比對過濾后的bam數據建索引，然后將bam和index文件轉移到新建的工作目錄中 - peak files
將macs2 peak calling 后的narrowPeak 文件移動到工作目錄中 - sampleSheet file
sampleSheet file是需要根據實驗設計和數據存儲地址等信息創建的一個csv格式文件（bam,peak文件分別在比對和call peak的步驟產生）。sampleSheet具體需要包含的信息如下：
SampleID: Identifier string for sample
Tissue, Factor, Condition: Identifier strings for up to three different factors (You will need to have all columns listed. If you don’t have infomation, then set values to NA)
Replicate: Replicate number of sample
bamReads: file path for BAM file containing aligned reads for ChIP sample
ControlID: an identifier string for the control sample
bamControl: file path for bam file containing aligned reads for control sample
Peaks: path for file containing peaks for sample
PeakCaller: Identifier string for peak caller used. Possible values include “raw”, “bed”, “narrow”, “macs”
1.2 軟件安裝
if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager") BiocManager::install("ChIPQC") #一開始我是進行上述安裝步驟，但后面頻頻出現報錯，后面找到原因是該程序包出現bug，shengqh/ChIPQC修復該bug，于是改為 BiocManager::install("shengqh/ChIPQC") #或者 library(devtools) install_github("shengqh/ChIPQC")安裝時最好設置CRAN使用國內鏡像，如清華大學鏡像：
options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")))1.3 Running ChIPQC
ChIPQC只需要三步就可以完成質量評估和報告生成。
首先載入包和sampleSheet信息
>library(ChIPQC) >#載入ChIPQC時需要先退出conda環境 #Load sample data samples <- read.csv('meta/samplesheet_chr12.csv') View(samples) #創建ChIPQC對象 #利用sampleSheet的信息讀取每個樣本的bam和narrowpeak文件，并計算質量評估值，結果存在一個對象里。 #Create ChIPQC object chipObj <- ChIPQC(samples, annotation="mm10") ## Create ChIPQC report ChIPQCreport(chipObj, reportName="ChIP QC report: Nanog and Pou5f1", reportFolder="ChIPQCreport")1.4 運行 ChIPQC時報錯及解決
我的數據分析時出現以下問題：
#問題一： [W::hts_idx_load2] The index file is older than the data file: #問題二： ```c Error in `combine_vars()`:! Faceting variables must have at least one value Run `rlang::last_error()` to see where the error occurred.問題一解決方法：
samtools index xxx.bam
問題二的原因是實驗提供的bam文件，染色質坐標的名字存在chr前綴，但注釋文件染色質坐標沒有前綴，只是數字。
解決方法（參考：https://github.com/dariober/cnv_facets/issues/9）：
以上解決方案行不通，出現以下報錯：
繼續采用染色質坐標加chr前綴的bam文件，參考https://www.jianshu.com/p/1c66cfb423af，直接查看chipObj 數據，但出現如下報錯：
chipObj Samples: 6 : Nanog-1 Nanog-2 ... RIP-Input-2 RIP-Input-3 Factor Control Replicate Peaks Nanog-1 Nanog RIP-Input-1 1 2359 Nanog-2 Nanog RIP-Input-2 2 1596 Nanog-3 Nanog RIP-Input-3 3 2299 Error in h(simpleError(msg, call)) : 在為't'函數選擇方法時評估'x'參數出了錯: 'names'屬性的長度[9]必需和矢量的長度[7]一樣最終解決：參考大佬：https://github.com/shengqh/ChIPQC
查找網友求助發現很多人遇到同樣報錯——error: “names’ attribute [9] must be the same length as the vector [7]”
shengqh/ChIPQC修復了這個bug——需要卸載原有ChIPQC，重新安裝。
2. ChIPQC: quality metrics report
ChIPQC是一個Bioconductor包，輸入文件包括BAM和peak文件，可以自動計算一些質量評估值，并產生質量報告。
2.1 ChIPQC report
- The QC summary table:
這張表有些列在之前沒有提過，SSD、RiP 和 RiBL 列下的信息是由ENCODE consortium開發和表述的指標，由CHIPQC進行計算。These will allow us to assess the distribution of signal within enriched regions, within/across expected annotations, across the whole genome, and within known artefact regions. 這些指標主要分類如下： - Read characteristics
- Enrichment of reads in peaks
- Peak signal strength
- Peak profiles
NOTE：這里給出的評估指標只是反映數據質量的好壞，符合閾值的并不意味著實驗是成功的，不符合閾值的也不一定意味著失敗。
2.1.1 Metrics: Read characteristics
這部分指標包括：read depth, read length, and duplication rate.需要注意read depth和length，尤其是在樣本之間似乎存在很大差異的情況下。由于我們已經過濾了 BAM 文件，所以duplication rate對我們來說意義不大。
2.1.2 Metrics: Enrichment of reads in peaks
我們通常會探索一些指標，包括 RiP、SSD 和 RiBL，有利于確定是否在峰值中具有很強的reads富集。
- RiP (Reads in Peaks)
This metric represents the percentage of reads that overlap ‘called peaks’（這句暫時不能理解）。它可以被認為是一種“信噪比”衡量文庫中由fragments from binding sites vs. background reads的比值。
RiP (also called FRiP，Fraction of reads in peaks) values will vary depending on the protein of interest:
- 成功富集的典型TF（尖峰/窄峰）的 RiP 約為 5% 或更高
- 數據質量好的Pol2（尖峰/窄峰和分散/寬峰的混合）將顯示 30% 或更高的 RiP
還有一些已知的 RiP < 1% 的良好數據集示例（即 RNAPIII 或結合少數位點的蛋白質）
在示例數據集中，與 Pou5f1 相比，Nanog replicates的 RiP 百分比更高，Pou5f1-rep2 非常低
這表明對 Nanog replicates有更高的富集程度，但仍然需要探索其他指標。
有兩個圖總結了the number of Reads in Peaks。 具有良好富集的 ChIP 樣本將有更高比例的reads overlapping called peaks。
盡管 Nanog 中的 RiP 較高，但與 Pou5f1 相比，Nanog 樣本的箱線圖顯示，replicates之間的分布完全不同，這可能是由于讀長和測序深度的不同導致的。
- SSD
該指標代表了reads覆蓋整個基因組的均勻性。 具體來說，它確定了信號堆積以基因組標準化為讀取總數的標準偏差。
在本示例中，相對于 Nanog replicates，Pou5f1 replicates具有更高的得分，這可能表明 Pou5f1 樣本的富集程度更高，但我們需要仔細查看 ChIPQC 的其余輸出，以確保 Pou5f1 中的高 SSD 不是由于某些未知的人為因素（unknown artifact）。
使用報告中的’Coverage histogram’可以可視化使用 SSD 探索的覆蓋率均勻性。 x 軸表示堿基對位置的讀取堆積高度，y 軸表示有多少位置具有該堆積高度。 這是對數規模。
具有良好富集的 ChIP 樣本應具有合理的“尾部”，或具有更多位置（y 軸上的更高值）的更高測序深度。具有低富集（即Input）的樣本，主要由背景讀數組成，在基因組中的大多數位置（y 軸上的高值）具有低堆積（低 x 軸值）。
在示例的數據集中，與 Pou5f1 相比，Nanog 樣本的尾部非常大（Heavy tail），尤其是rep2。“重尾”是指曲線比指數曲線大（heavier），曲線下體積更大。這表明 Nanog 樣本深度更高的在基因組中的位置更多。然而，Pou5f1 的 SSD 分數更高。 當 SSD 很高但覆蓋率看起來很低時，這可能是由于存在高深度的廣泛區域（large regions of high depth）和基因組區域黑名單的標志（ a flag for blacklisting of genomic regions）。
此處補充一下blacklisting of genomic regions的講解
- RiBL (Reads overlapping in Blacklisted Regions)
具有已知人為高信號的讀數重疊區域的百分比（可能是由于過度的非結構化異常讀數映射）。Lower RiBL percentages are better than higher.
在示例的實驗中，RiBL 百分比看起來是合理的，因為它們并不高。較高的 SSD 值可能是由于更容易碎裂的開放染色質區域或超 ChIPable 區域，這些區域是在許多不相關蛋白質的 ChIP-seq 實驗中富集的高表達位點 - 系統性假陽性。
注意：如果你在peak calling之前過濾掉了blacklisted regions，并且這些過濾后的 BAM 文件用作 ChIPQC 的輸入，則無需評估此指標。
“blacklists regions”是如何編制的？ 這些blacklists 是通過結合使用自動啟發式和手動管理從大量數據中根據經驗得出的。為包括人類、小鼠、蠕蟲和蒼蠅在內的各種物種和基因組版本生成了blacklists。The lists can be downloaded here。對于human，使用了80個開放染色質軌道(DNase和FAIRE數據集)和12個ChIP-seq Input/Control 軌道，共跨越約60個細胞系。這些黑名單適用于基于短讀測序(20-100bp reads)的功能基因組數據。這些并不直接適用于RNA-seq或任何其他轉錄組數據類型。
下圖顯示了blacklisting的效果圖，其中包含在blacklists區域內部（深藍色）或外部（淺藍色）的讀取比例。 如果blacklists區域在call peaking之前已經被過濾掉，則該圖將不會提供有用信息。
2.1.3 Metrics: Peak profiles
- Relative Enrichment of Genomic Intervals (REGI)
使用基因組區域確定的called peaks以及基因組注釋信息，我們可以獲得根據各種基因組特征查看讀取映射的位置。
This is represented as a heatmap showing the enrichment of reads compared to the background levels of the feature. This plot is useful when you expect enrichment of specific genomic regions.
在示例的數據集中，“Promoters500”和“All5UTRs”類別的富集水平最高，這符合我們對 Nanog 和 Pou5f1 應該作為轉錄因子結合的預期。
- Peak Profile
The peak profile plot shows the average peak profiles,以每個峰的峰頂（最高堆積點）為中心。
些圖譜的形狀可能因所研究的標記類型而異——轉錄因子、組蛋白標記或其他 DNA 結合蛋白（如聚合酶）——但類似的標記通常在成功的 ChIPs 中具有獨特的特征。
- Sample similarity
最后，圖中顯示了樣本的相似程度。當我們在差異富集分析中繪制它們時，我們將更仔細地查看它們。
2.1.4 Metrics: Peak signal strength
與峰值信號強度（Peak signal strength）相關的指標是 FragLength 和 RelCC(also called Relative strand cross-correlation coefficient or RSC)。這兩個值都是通過計算strand cross-correlation確定的，這是一個獨立于peak calling的質量評估（quality metric）。
所有 ChIP 樣本的 RelCC 值均大于 1 表明信噪比良好，并且 FragL 值應與您在文庫制備期間的大小選擇步驟中選擇的片段長度大致相同。
- Strand cross-correlation
高質量的 ChIP-seq 實驗將在感興趣的蛋白質結合的位置產生顯著的富集 DNA 序列標簽/讀數聚類； 期望是我們可以觀察到正向和反向鏈上讀取（序列標簽）的雙峰富集（bimodal enrichment of reads）。
How are the Cross-Correlation scores calculated?
以一個小的基因組窗口為例，讓我們來看看下面cross-correlation的細節。 重要的是要注意，通過一次系統地將負鏈移動 k 個堿基對，將cross-correlation度量計算為每個互補堿基（即在負鏈和正鏈上）的覆蓋率之間的 Pearson 線性相關性。 這種轉變一遍又一遍地進行，以獲得基因組給定區域的correlation。
圖 1：鏈位移為零時，兩個向量之間的 Pearson 相關性為 0。
圖 2：鏈位移為 100bp 時，兩個向量之間的 Pearson 相關性為 0.389。
圖 3：鏈位移為 175bp 時，兩個向量之間的 Pearson 相關性為 0.831。
當這個過程完成時，將生成一個值表，將每個堿基對偏移映射到 Pearson 相關值。 為每個染色體的每個峰計算這些 Pearson 相關值，并將值乘以比例因子，然后在所有染色體上求和。
一旦計算出最終的cross-correlation values，就可以根據偏移值（X 軸）繪制cross-correlation values（Y 軸）以生成cross-correlation圖！
cross-correlation圖通常會產生兩個峰：對應于主要片段長度(predominant fragment length)（最高相關值）的富集峰和對應于讀取長度(read length)的峰（“幻影/phantom”峰）。 示例的cross-correlation圖：
在示例中，Nanog 中的最大相關值（較大峰的最高點）高于 Pou5f1，表明信號量更高。 該最大相關性的相應移位值（x 軸）為我們提供了估計的片段長度（estimated fragment length）。
RelCC（或 RSC）值是使用最小和最大cross-correlation values計算的。 要獲得有關如何計算 RSC 和 NSC（另一種cross-correlation values的度量）的更多詳細信息，除了圍繞“幻象峰/phantom peak”現象的討論之外，請查看這些材料。 RSC 值低可能是由于 ChIP 失敗或質量差、讀取序列質量低以及因此存在大量錯誤映射、測序深度淺或這些原因的組合。 此外，可能由于生物學原因（即真正僅結合特定組織類型中的少數位點的因素）具有少量結合位點（< 200）的數據集將輸出低 RSC 分數。
Examples of Cross-Correlation Plots
Strong signal:
與read-length peak相比，高質量的ChIP-seq數據集往往具有更大的fragment-length peak峰值。下面使用來自人類細胞中 CTCF（鋅指轉錄因子）的數據顯示了一個強信號的示例。 使用好的抗體，轉錄因子通常會產生 45,000 - 60,000 個peaks。紅色垂直線顯示真實峰移處的主峰，藍色垂直線處的小凸起代表reads長度。
Weak signal:
下面使用 Pol2 數據演示了一個較弱信號的示例。 在這里，這種抗體不是很有效，呈現廣泛的分散峰。 觀察到兩個峰：一個在真實峰移（true peak shift）處（~185-200 bp），另一個在讀取長度處（read-length peak）。 對于弱信號數據集，讀取長度峰值（read-length peak）將占主導。
No signal:
A failed experiment will resemble a cross-correlation plot using input only, in which we observe little or no peak for fragment length.
2.2 Final takehome from ChIPQC
首先是用于單獨評估每個樣本，以確保我們觀察到的值足夠好，我們可以放心地繼續前進，但還要進一步比較replicates之間和groups之間的異同。
-
Within sample group:
每個樣本組似乎都有一個replicate表現出比另一個更強的信號。 這種差異在基于the cross-correlation plots and coverage plots的 Nanog replicate中更為明顯。 盡管重復之間的信號/富集量存在差異，但看到類似的趨勢。 然而，顯示完全不同的cross-correlation plots的兩個重復表明出現了問題。 -
Between sample groups:
比較一個樣本組與另一個樣本組之間的指標，很難得出一個是否比另一個更好的結論。Pou5f1的SSD和RelCC得分較高，表明富集較好，而cross-correlation plots and coverage plots則表明Nanog樣品中有更多的信號。這兩組之間的差異會在之后differential enrichment and visualization中討論。
2.3 低質量數據的可能來源是什么?
- 免疫沉淀的強度/效率和特異性
ChIP實驗的質量最終取決于抗體的特異性和親和沉淀步驟中達到的富集程度。抗體缺陷主要有兩種類型: - 碎片化/消化
進行超聲處理的方式可能會導致不同的片段大小分布，從而導致樣品特異性染色質配置產生偏差。 因此，如果未與 ChIP 樣品一起進行超聲處理，則不建議使用單個Input樣品作為 ChIP-seq peak calling的對照。
總結
以上是生活随笔為你收集整理的ChIPQC——对ChIP-seq的质量评估的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 视频监控安防平台-国标35114（GB3
- 下一篇： [ 数据集 ] VOC 2012 数据集