第4篇：对ChIP-seq的质量评估

学习目标：

探讨ChIP-seq数据质量低的来源
主要熟悉ChIPQC对ChIP-seq进行Peak calling之后的质量评估
次要理解链交叉相关性（ strand cross-correlation）
次要使用phantompeakqualtools计算交叉相关性和其他相关的质控度量值

文章目录

第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 phantompeakqualtools

1.2 使用phantompeakqualtools计算交叉相关性和其他相关的质控度量值

参考：第4篇：对ATAC-Seq/ChIP-seq的质量评估（一）——phantompeakqualtools

三、ChIPQC对数据进行质量评估

学习目标

讨论ChIP-seq数据质量评估的其他方法
用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）：
#使用samtools view -H 和shell删除bam文件中的chr前缀
samtools view -H RIP-Input-1_markdup.bam | sed 's/chr//' | samtools reheader - RIP-Input-1_markdup.bam > test.CHR.bam

以上解决方案行不通，出现以下报错：

继续采用染色质坐标加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，重新安装。
library(ChIPQC)
#查看该包的安装位置
find.package("ChIPQC")
#删除包
remove.packages("ChIPQC")
#安装新包
install.packages("devtools")
library(devtools)
install_github("shengqh/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之间的分布完全不同，这可能是由于读长和测序深度的不同导致的。

该指标代表了reads覆盖整个基因组的均匀性。具体来说，它确定了信号堆积以基因组标准化为读取总数的标准偏差。

对于经过 IP 处理的样本，预计读数丰富或覆盖率高的区域以及覆盖率较低的其他区域。而对于对照样本，预计整个基因组的覆盖率差异较小。 “良好”或富集的样本通常具有显着的读取堆积区域（read pile-up）（覆盖范围的较大差异），因此更高的 SSD 能表示更好的富集。
除了真正的 ChIP 富集外，SSD 分数对人为高信号区域（regions of artificially high signal）也很敏感。因此，高 SSD 可能是 ChIP 富集或黑名单区域（blacklisted regions）中某些人为高信号的结果。

在本示例中，相对于 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.

列入Blacklists的区域通常是唯一可映射的，因此简单的可映射性过滤器不会删除它们。这些区域通常存在于特定类型的重复序列中，例如着丝粒、端粒和卫星重复序列。
RiBL 分数可作为 ChIP 或Input背景信号水平的指南。这些区域约占基因组的 0.5%，但占总信号的比例很高（> 10%）。
来自blacklisted regions的信号已被证明会混淆 peak callers and fragment length estimation。我们应该跟踪和过滤映射到这些区域的reads。

在示例的实验中，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的质量评估相关推荐

NanoPlot：三代纳米孔测序数据质量评估
简介二代测序最常用的质量评估软件是FastQC,多样本时可进一步结合MultiQC.此外速度超快的fastp也特别推荐,而且包括质量评估.质量控制等功能,可以说是国产软件之光,详见下方详细教程: 数 ...
人脸质量评估网络推荐
腾讯优图开源的: https://github.com/Tencent/TFace/tree/quality https://github.com/deepcam-cn/FaceQuality 效果图 ...
MNA-CNN: 如何在美学质量评估中储存照片原始信息
传统基于深度学习的图片美学评估的要求输入是fixed-size,待评估的照片经常会经过裁剪.缩放和充填等变换,造成图片质量的下降. 那如何可以实现对原始图片尺寸进行美学质量评估呢? 示例:由上图,可以 ...
字节跳动VQScore算法拿下ICME 2021“压缩UGC视频质量评估”比赛第一名
在ICME 2021国际视频质量评估算法竞赛中,字节跳动拿下第一. 在多媒体领域世界顶级学术会议ICME 2021的「压缩UGC视频质量评估」比赛中,字节跳动旗下火山引擎多媒体实验室组成的" ...
无参考质量评估在视频增强的进展与应用
无参考质量评估在许多无法取得参考信息的实际系统中应用广泛且十分重要,本文整理自腾讯音视频实验室的高孟平在LiveVideoStackCon 2019上海大会中的分享,详细介绍了腾讯音视频实验室团队如何 ...
teablue数据分析_基于大数据分析的茶叶质量评估
引言茶叶是茶加工企业生存和发展的基本保障,是茶产品能够顺利发展的重要资源.茶叶的质量和产量,通常会遭受各种大气条件以及生态条件和生产措施的影响.为了保证茶叶质量,需要监测茶叶生产的整个周期. ...
企业如何进行数据质量评估
一般来说,当企业有了全新的业务需求.重大的技术变更,又或者从一个新的数据来源获取了全新的数据,并期望将它应用在一个具体的业务中的时候,我们都需要进行比较完整的数据质量分析.数据质量评估步骤如下: 1. ...
AAAI2023 | 基于课程学习的机器翻译质量评估去噪预训练
每天给你送来NLP技术干货! 来自:南大NLP 点击这里进群->加入NLP交流群 01 研究动机机器翻译质量评估(Quality Estimation, QE)任务指在没有参考译文的情况下,仅 ...
讯飞机器翻译质量评估挑战赛Baseline(PaddlePaddle)
赛题简介比赛地址:点击直达举办方:科大讯飞股份有限公司任务类型:质量评估(QE).自然语言回归赛事背景机器翻译质量评估(QE)指在没有人工翻译参考下对机器翻译系统译文进行自动打分.一方面,Q ...
腾讯多媒体实验室：基于三维卷积神经网络的全参考视频质量评估算法
腾讯有多个视频业务线,点播视频有腾讯视频.企鹅影视,短视频有微视.K歌,直播类有Now直播.企鹅电竞,实时传输类有QQ和微信的音视频通话.无线投屏和腾讯会议等. 用户对不同的产品有不同程度的期待:比如 ...

ChIPQC——对ChIP-seq的质量评估