ChIPQC——对ChIP-seq的质量评估
第4篇：對(duì)ChIP-seq的質(zhì)量評(píng)估
學(xué)習(xí)目標(biāo)：
文章目錄
- 第4篇：對(duì)ChIP-seq的質(zhì)量評(píng)估
- 一、ChIP-Seq質(zhì)量評(píng)估
- 二、鏈交叉相關(guān)（Strand cross-correlation）
- 1. phantompeakqualtools
- 1.1 軟件安裝
- 1.2 使用phantompeakqualtools計(jì)算交叉相關(guān)性和其他相關(guān)的質(zhì)控度量值
- 三、ChIPQC對(duì)數(shù)據(jù)進(jìn)行質(zhì)量評(píng)估
- 1. 實(shí)際操作
- 1.1 設(shè)置
- 1.2 軟件安裝
- 1.3 Running ChIPQC
- 1.4 運(yùn)行 ChIPQC時(shí)報(bào)錯(cuò)及解決
- 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 低質(zhì)量數(shù)據(jù)的可能來源是什么?
一、ChIP-Seq質(zhì)量評(píng)估
在下游分析前，最好是先對(duì)peak calling 后的ChIP-Seq數(shù)據(jù)進(jìn)行質(zhì)量評(píng)估。
二、鏈交叉相關(guān)（Strand cross-correlation）
1. phantompeakqualtools
1.1 軟件安裝
安裝方法一：https://anaconda.org/bioconda/phantompeakqualtools
安裝方法二：https://github.com/kundajelab/phantompeakqualtools
這里采用方法一進(jìn)行安裝：
conda install -c bioconda phantompeakqualtools1.2 使用phantompeakqualtools計(jì)算交叉相關(guān)性和其他相關(guān)的質(zhì)控度量值
參考：第4篇：對(duì)ATAC-Seq/ChIP-seq的質(zhì)量評(píng)估（一）——phantompeakqualtools
三、ChIPQC對(duì)數(shù)據(jù)進(jìn)行質(zhì)量評(píng)估
學(xué)習(xí)目標(biāo)
參考一：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. 實(shí)際操作
1.1 設(shè)置
- mkdir一個(gè)新文件夾進(jìn)行操作，將需要的文件轉(zhuǎn)移到新文件夾中；需要的文件有：
- BAM files
首先對(duì)比對(duì)過濾后的bam數(shù)據(jù)建索引，然后將bam和index文件轉(zhuǎn)移到新建的工作目錄中 - peak files
將macs2 peak calling 后的narrowPeak 文件移動(dòng)到工作目錄中 - sampleSheet file
sampleSheet file是需要根據(jù)實(shí)驗(yàn)設(shè)計(jì)和數(shù)據(jù)存儲(chǔ)地址等信息創(chuàng)建的一個(gè)csv格式文件（bam,peak文件分別在比對(duì)和call peak的步驟產(chǎn)生）。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") #一開始我是進(jìn)行上述安裝步驟，但后面頻頻出現(xiàn)報(bào)錯(cuò)，后面找到原因是該程序包出現(xiàn)bug，shengqh/ChIPQC修復(fù)該bug，于是改為 BiocManager::install("shengqh/ChIPQC") #或者 library(devtools) install_github("shengqh/ChIPQC")安裝時(shí)最好設(shè)置CRAN使用國內(nèi)鏡像，如清華大學(xué)鏡像：
options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")))1.3 Running ChIPQC
ChIPQC只需要三步就可以完成質(zhì)量評(píng)估和報(bào)告生成。
首先載入包和sampleSheet信息
>library(ChIPQC) >#載入ChIPQC時(shí)需要先退出conda環(huán)境 #Load sample data samples <- read.csv('meta/samplesheet_chr12.csv') View(samples) #創(chuàng)建ChIPQC對(duì)象 #利用sampleSheet的信息讀取每個(gè)樣本的bam和narrowpeak文件，并計(jì)算質(zhì)量評(píng)估值，結(jié)果存在一個(gè)對(duì)象里。 #Create ChIPQC object chipObj <- ChIPQC(samples, annotation="mm10") ## Create ChIPQC report ChIPQCreport(chipObj, reportName="ChIP QC report: Nanog and Pou5f1", reportFolder="ChIPQCreport")1.4 運(yùn)行 ChIPQC時(shí)報(bào)錯(cuò)及解決
我的數(shù)據(jù)分析時(shí)出現(xiàn)以下問題：
#問題一： [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
問題二的原因是實(shí)驗(yàn)提供的bam文件，染色質(zhì)坐標(biāo)的名字存在chr前綴，但注釋文件染色質(zhì)坐標(biāo)沒有前綴，只是數(shù)字。
解決方法（參考：https://github.com/dariober/cnv_facets/issues/9）：
以上解決方案行不通，出現(xiàn)以下報(bào)錯(cuò)：
繼續(xù)采用染色質(zhì)坐標(biāo)加chr前綴的bam文件，參考https://www.jianshu.com/p/1c66cfb423af，直接查看chipObj 數(shù)據(jù)，但出現(xiàn)如下報(bào)錯(cuò)：
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'函數(shù)選擇方法時(shí)評(píng)估'x'參數(shù)出了錯(cuò): 'names'屬性的長度[9]必需和矢量的長度[7]一樣最終解決：參考大佬：https://github.com/shengqh/ChIPQC
查找網(wǎng)友求助發(fā)現(xiàn)很多人遇到同樣報(bào)錯(cuò)——error: “names’ attribute [9] must be the same length as the vector [7]”
shengqh/ChIPQC修復(fù)了這個(gè)bug——需要卸載原有ChIPQC，重新安裝。
2. ChIPQC: quality metrics report
ChIPQC是一個(gè)Bioconductor包，輸入文件包括BAM和peak文件，可以自動(dòng)計(jì)算一些質(zhì)量評(píng)估值，并產(chǎn)生質(zhì)量報(bào)告。
2.1 ChIPQC report
- The QC summary table:
這張表有些列在之前沒有提過，SSD、RiP 和 RiBL 列下的信息是由ENCODE consortium開發(fā)和表述的指標(biāo)，由CHIPQC進(jìn)行計(jì)算。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. 這些指標(biāo)主要分類如下： - Read characteristics
- Enrichment of reads in peaks
- Peak signal strength
- Peak profiles
NOTE：這里給出的評(píng)估指標(biāo)只是反映數(shù)據(jù)質(zhì)量的好壞，符合閾值的并不意味著實(shí)驗(yàn)是成功的，不符合閾值的也不一定意味著失敗。
2.1.1 Metrics: Read characteristics
這部分指標(biāo)包括：read depth, read length, and duplication rate.需要注意read depth和length，尤其是在樣本之間似乎存在很大差異的情況下。由于我們已經(jīng)過濾了 BAM 文件，所以duplication rate對(duì)我們來說意義不大。
2.1.2 Metrics: Enrichment of reads in peaks
我們通常會(huì)探索一些指標(biāo)，包括 RiP、SSD 和 RiBL，有利于確定是否在峰值中具有很強(qiáng)的reads富集。
- RiP (Reads in Peaks)
This metric represents the percentage of reads that overlap ‘called peaks’（這句暫時(shí)不能理解）。它可以被認(rèn)為是一種“信噪比”衡量文庫中由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% 或更高
- 數(shù)據(jù)質(zhì)量好的Pol2（尖峰/窄峰和分散/寬峰的混合）將顯示 30% 或更高的 RiP
還有一些已知的 RiP < 1% 的良好數(shù)據(jù)集示例（即 RNAPIII 或結(jié)合少數(shù)位點(diǎn)的蛋白質(zhì)）
在示例數(shù)據(jù)集中，與 Pou5f1 相比，Nanog replicates的 RiP 百分比更高，Pou5f1-rep2 非常低
這表明對(duì) Nanog replicates有更高的富集程度，但仍然需要探索其他指標(biāo)。
有兩個(gè)圖總結(jié)了the number of Reads in Peaks。 具有良好富集的 ChIP 樣本將有更高比例的reads overlapping called peaks。
盡管 Nanog 中的 RiP 較高，但與 Pou5f1 相比，Nanog 樣本的箱線圖顯示，replicates之間的分布完全不同，這可能是由于讀長和測序深度的不同導(dǎo)致的。
- SSD
該指標(biāo)代表了reads覆蓋整個(gè)基因組的均勻性。 具體來說，它確定了信號(hào)堆積以基因組標(biāo)準(zhǔn)化為讀取總數(shù)的標(biāo)準(zhǔn)偏差。
在本示例中，相對(duì)于 Nanog replicates，Pou5f1 replicates具有更高的得分，這可能表明 Pou5f1 樣本的富集程度更高，但我們需要仔細(xì)查看 ChIPQC 的其余輸出，以確保 Pou5f1 中的高 SSD 不是由于某些未知的人為因素（unknown artifact）。
使用報(bào)告中的’Coverage histogram’可以可視化使用 SSD 探索的覆蓋率均勻性。 x 軸表示堿基對(duì)位置的讀取堆積高度，y 軸表示有多少位置具有該堆積高度。 這是對(duì)數(shù)規(guī)模。
具有良好富集的 ChIP 樣本應(yīng)具有合理的“尾部”，或具有更多位置（y 軸上的更高值）的更高測序深度。具有低富集（即Input）的樣本，主要由背景讀數(shù)組成，在基因組中的大多數(shù)位置（y 軸上的高值）具有低堆積（低 x 軸值）。
在示例的數(shù)據(jù)集中，與 Pou5f1 相比，Nanog 樣本的尾部非常大（Heavy tail），尤其是rep2。“重尾”是指曲線比指數(shù)曲線大（heavier），曲線下體積更大。這表明 Nanog 樣本深度更高的在基因組中的位置更多。然而，Pou5f1 的 SSD 分?jǐn)?shù)更高。 當(dāng) SSD 很高但覆蓋率看起來很低時(shí)，這可能是由于存在高深度的廣泛區(qū)域（large regions of high depth）和基因組區(qū)域黑名單的標(biāo)志（ a flag for blacklisting of genomic regions）。
此處補(bǔ)充一下blacklisting of genomic regions的講解
- RiBL (Reads overlapping in Blacklisted Regions)
具有已知人為高信號(hào)的讀數(shù)重疊區(qū)域的百分比（可能是由于過度的非結(jié)構(gòu)化異常讀數(shù)映射）。Lower RiBL percentages are better than higher.
在示例的實(shí)驗(yàn)中，RiBL 百分比看起來是合理的，因?yàn)樗鼈儾⒉桓摺］^高的 SSD 值可能是由于更容易碎裂的開放染色質(zhì)區(qū)域或超 ChIPable 區(qū)域，這些區(qū)域是在許多不相關(guān)蛋白質(zhì)的 ChIP-seq 實(shí)驗(yàn)中富集的高表達(dá)位點(diǎn) - 系統(tǒng)性假陽性。
注意：如果你在peak calling之前過濾掉了blacklisted regions，并且這些過濾后的 BAM 文件用作 ChIPQC 的輸入，則無需評(píng)估此指標(biāo)。
“blacklists regions”是如何編制的？ 這些blacklists 是通過結(jié)合使用自動(dòng)啟發(fā)式和手動(dòng)管理從大量數(shù)據(jù)中根據(jù)經(jīng)驗(yàn)得出的。為包括人類、小鼠、蠕蟲和蒼蠅在內(nèi)的各種物種和基因組版本生成了blacklists。The lists can be downloaded here。對(duì)于human，使用了80個(gè)開放染色質(zhì)軌道(DNase和FAIRE數(shù)據(jù)集)和12個(gè)ChIP-seq Input/Control 軌道，共跨越約60個(gè)細(xì)胞系。這些黑名單適用于基于短讀測序(20-100bp reads)的功能基因組數(shù)據(jù)。這些并不直接適用于RNA-seq或任何其他轉(zhuǎn)錄組數(shù)據(jù)類型。
下圖顯示了blacklisting的效果圖，其中包含在blacklists區(qū)域內(nèi)部（深藍(lán)色）或外部（淺藍(lán)色）的讀取比例。 如果blacklists區(qū)域在call peaking之前已經(jīng)被過濾掉，則該圖將不會(huì)提供有用信息。
2.1.3 Metrics: Peak profiles
- Relative Enrichment of Genomic Intervals (REGI)
使用基因組區(qū)域確定的called peaks以及基因組注釋信息，我們可以獲得根據(jù)各種基因組特征查看讀取映射的位置。
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.
在示例的數(shù)據(jù)集中，“Promoters500”和“All5UTRs”類別的富集水平最高，這符合我們對(duì) Nanog 和 Pou5f1 應(yīng)該作為轉(zhuǎn)錄因子結(jié)合的預(yù)期。
- Peak Profile
The peak profile plot shows the average peak profiles,以每個(gè)峰的峰頂（最高堆積點(diǎn)）為中心。
些圖譜的形狀可能因所研究的標(biāo)記類型而異——轉(zhuǎn)錄因子、組蛋白標(biāo)記或其他 DNA 結(jié)合蛋白（如聚合酶）——但類似的標(biāo)記通常在成功的 ChIPs 中具有獨(dú)特的特征。
- Sample similarity
最后，圖中顯示了樣本的相似程度。當(dāng)我們?cè)诓町惛患治鲋欣L制它們時(shí)，我們將更仔細(xì)地查看它們。
2.1.4 Metrics: Peak signal strength
與峰值信號(hào)強(qiáng)度（Peak signal strength）相關(guān)的指標(biāo)是 FragLength 和 RelCC(also called Relative strand cross-correlation coefficient or RSC)。這兩個(gè)值都是通過計(jì)算strand cross-correlation確定的，這是一個(gè)獨(dú)立于peak calling的質(zhì)量評(píng)估（quality metric）。
所有 ChIP 樣本的 RelCC 值均大于 1 表明信噪比良好，并且 FragL 值應(yīng)與您在文庫制備期間的大小選擇步驟中選擇的片段長度大致相同。
- Strand cross-correlation
高質(zhì)量的 ChIP-seq 實(shí)驗(yàn)將在感興趣的蛋白質(zhì)結(jié)合的位置產(chǎn)生顯著的富集 DNA 序列標(biāo)簽/讀數(shù)聚類； 期望是我們可以觀察到正向和反向鏈上讀取（序列標(biāo)簽）的雙峰富集（bimodal enrichment of reads）。
How are the Cross-Correlation scores calculated?
以一個(gè)小的基因組窗口為例，讓我們來看看下面cross-correlation的細(xì)節(jié)。 重要的是要注意，通過一次系統(tǒng)地將負(fù)鏈移動(dòng) k 個(gè)堿基對(duì)，將cross-correlation度量計(jì)算為每個(gè)互補(bǔ)堿基（即在負(fù)鏈和正鏈上）的覆蓋率之間的 Pearson 線性相關(guān)性。 這種轉(zhuǎn)變一遍又一遍地進(jìn)行，以獲得基因組給定區(qū)域的correlation。
圖 1：鏈位移為零時(shí)，兩個(gè)向量之間的 Pearson 相關(guān)性為 0。
圖 2：鏈位移為 100bp 時(shí)，兩個(gè)向量之間的 Pearson 相關(guān)性為 0.389。
圖 3：鏈位移為 175bp 時(shí)，兩個(gè)向量之間的 Pearson 相關(guān)性為 0.831。
當(dāng)這個(gè)過程完成時(shí)，將生成一個(gè)值表，將每個(gè)堿基對(duì)偏移映射到 Pearson 相關(guān)值。 為每個(gè)染色體的每個(gè)峰計(jì)算這些 Pearson 相關(guān)值，并將值乘以比例因子，然后在所有染色體上求和。
一旦計(jì)算出最終的cross-correlation values，就可以根據(jù)偏移值（X 軸）繪制cross-correlation values（Y 軸）以生成cross-correlation圖！
cross-correlation圖通常會(huì)產(chǎn)生兩個(gè)峰：對(duì)應(yīng)于主要片段長度(predominant fragment length)（最高相關(guān)值）的富集峰和對(duì)應(yīng)于讀取長度(read length)的峰（“幻影/phantom”峰）。 示例的cross-correlation圖：
在示例中，Nanog 中的最大相關(guān)值（較大峰的最高點(diǎn)）高于 Pou5f1，表明信號(hào)量更高。 該最大相關(guān)性的相應(yīng)移位值（x 軸）為我們提供了估計(jì)的片段長度（estimated fragment length）。
RelCC（或 RSC）值是使用最小和最大cross-correlation values計(jì)算的。 要獲得有關(guān)如何計(jì)算 RSC 和 NSC（另一種cross-correlation values的度量）的更多詳細(xì)信息，除了圍繞“幻象峰/phantom peak”現(xiàn)象的討論之外，請(qǐng)查看這些材料。 RSC 值低可能是由于 ChIP 失敗或質(zhì)量差、讀取序列質(zhì)量低以及因此存在大量錯(cuò)誤映射、測序深度淺或這些原因的組合。 此外，可能由于生物學(xué)原因（即真正僅結(jié)合特定組織類型中的少數(shù)位點(diǎn)的因素）具有少量結(jié)合位點(diǎn)（< 200）的數(shù)據(jù)集將輸出低 RSC 分?jǐn)?shù)。
Examples of Cross-Correlation Plots
Strong signal:
與read-length peak相比，高質(zhì)量的ChIP-seq數(shù)據(jù)集往往具有更大的fragment-length peak峰值。下面使用來自人類細(xì)胞中 CTCF（鋅指轉(zhuǎn)錄因子）的數(shù)據(jù)顯示了一個(gè)強(qiáng)信號(hào)的示例。 使用好的抗體，轉(zhuǎn)錄因子通常會(huì)產(chǎn)生 45,000 - 60,000 個(gè)peaks。紅色垂直線顯示真實(shí)峰移處的主峰，藍(lán)色垂直線處的小凸起代表reads長度。
Weak signal:
下面使用 Pol2 數(shù)據(jù)演示了一個(gè)較弱信號(hào)的示例。 在這里，這種抗體不是很有效，呈現(xiàn)廣泛的分散峰。 觀察到兩個(gè)峰：一個(gè)在真實(shí)峰移（true peak shift）處（~185-200 bp），另一個(gè)在讀取長度處（read-length peak）。 對(duì)于弱信號(hào)數(shù)據(jù)集，讀取長度峰值（read-length peak）將占主導(dǎo)。
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
首先是用于單獨(dú)評(píng)估每個(gè)樣本，以確保我們觀察到的值足夠好，我們可以放心地繼續(xù)前進(jìn)，但還要進(jìn)一步比較replicates之間和groups之間的異同。
-
Within sample group:
每個(gè)樣本組似乎都有一個(gè)replicate表現(xiàn)出比另一個(gè)更強(qiáng)的信號(hào)。 這種差異在基于the cross-correlation plots and coverage plots的 Nanog replicate中更為明顯。 盡管重復(fù)之間的信號(hào)/富集量存在差異，但看到類似的趨勢。 然而，顯示完全不同的cross-correlation plots的兩個(gè)重復(fù)表明出現(xiàn)了問題。 -
Between sample groups:
比較一個(gè)樣本組與另一個(gè)樣本組之間的指標(biāo)，很難得出一個(gè)是否比另一個(gè)更好的結(jié)論。Pou5f1的SSD和RelCC得分較高，表明富集較好，而cross-correlation plots and coverage plots則表明Nanog樣品中有更多的信號(hào)。這兩組之間的差異會(huì)在之后differential enrichment and visualization中討論。
2.3 低質(zhì)量數(shù)據(jù)的可能來源是什么?
- 免疫沉淀的強(qiáng)度/效率和特異性
ChIP實(shí)驗(yàn)的質(zhì)量最終取決于抗體的特異性和親和沉淀步驟中達(dá)到的富集程度。抗體缺陷主要有兩種類型: - 碎片化/消化
進(jìn)行超聲處理的方式可能會(huì)導(dǎo)致不同的片段大小分布，從而導(dǎo)致樣品特異性染色質(zhì)配置產(chǎn)生偏差。 因此，如果未與 ChIP 樣品一起進(jìn)行超聲處理，則不建議使用單個(gè)Input樣品作為 ChIP-seq peak calling的對(duì)照。
總結(jié)
以上是生活随笔為你收集整理的ChIPQC——对ChIP-seq的质量评估的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 视频监控安防平台-国标35114（GB3
- 下一篇： [ 数据集 ] VOC 2012 数据集