ATAC-Seq分析教程：对ATAC-Seq/ChIP-seq的质量评估（一）phantompeakqualtools—科研必备表观遗传学知识

ATAC-Seq剖析教程系列

ATAC-Seq剖析教程：ATAC-seq的布景介绍以及与ChIP-Seq的异同

ATAC-Seq剖析教程：原始数据的质控、比对和过滤

ATAC-Seq剖析教程：用MACS2软件call peaks

ATAC-Seq剖析教程：对ATAC-Seq/ChIP-seq的质量评价（一）phantompeakqualtools

ATAC-Seq剖析教程：对ATAC-Seq/ChIP-seq的质量评价（二）ChIPQC

ATAC-Seq剖析教程：重复样本的处理-IDR

ATAC-Seq剖析教程：用ChIPseeker对peaks进行注释和可视化

ATAC-Seq剖析教程：用网页版工具做功能剖析和motif剖析

ATAC-Seq剖析教程：差异peaks剖析——DiffBind

ATAC-Seq剖析教程：ATAC-Seq、ChIP-Seq、RNA-Seq整合剖析

学习目标：

讨论ChIP-seq数据质量低的来源
理解链穿插相关性（ strand cross-correlation）
运用phantompeakqualtools核算穿插相关性和其他相关的质控衡量值
评价穿插相关图

ChIP-Seq质量评价

在下游剖析前，最好是先对peak calling 后的ChIP-Seq数据进行质量评价。

ATAC-Seq剖析教程：对ATAC-Seq/ChIP-seq的质量评价（一）phantompeakqualtools

链穿插相关（Strand cross-correlation）

链穿插相关是一个有用的评价ChIP-Seq质量的办法，它不依靠于peak calling，而是根据ChIP-Seq试验。假如ChIP-Seq试验成功，DNA富集序列标签（蛋白质相互作用的序列）会在reads的双峰富会集发生明显的集合。
发生reads的双峰富集的原因如下：
在ChIP-Seq试验中，DNA被片段化，蛋白质结合的片段会被免疫沉淀，所以发生了有蛋白质结合的DNA片段（fragments ）。
DNA的正链从5’端开始被测序（如下图赤色reads），DNA负链也从5’结尾被测序发生如下图所示的蓝色reads。

因为从DNA片段的5′结尾测序，使链reads的富集（下图中的蓝色部分）与负链reads的富集（下图赤色部分）有少量的相互抵消区域。咱们需要确认峰位移多少碱基数目能够在两个峰间发生最大的相关性。咱们能够用穿插相关的衡量值（cross-correlation metric）核算发生最大相关的位移。

穿插相关性衡量值
穿插相关衡量是在Watson移动k个碱基后，核算Crick链与Watson链之间的Pearson线性相关。如下面的示意图：
首先在链位移为0时，两个向量之间的Pearson相关值为0.539。

在链位移5bp，两个向量之间的Pearson相关值为0.931。

持续移动这些向量，关于每个链位移核算一个相关值。

最后，咱们将有一个每个碱基对移位与皮尔森相关值的对应表。这是针对每个染色体的每一个峰核算的，然后该值乘以一个缩放因子，再对所有染色体的值相加，就能够绘穿插相关值（y轴）相关于移位值（x轴）生成的穿插相关图。
典型的穿插相关图会发生两个峰：一个富集峰与主要的片段长度（predominant fragment length）相关（高相关性），另一个与read 长度（read length）相关，这个峰也被称为虚幻峰（“phantom” peak）。

质量好的ChIP-Seq数据集倾向发生一个大的片段长度峰（fragment-length peak），下图展示了一个来自于人细胞CTCF（zinc-finger transcription factor）的强信号。假如有好的抗体，转录因子一般发生45,000~60,000个peaks。下图赤色的垂直线表明主峰的真实位移，蓝色的垂直线处有一个小的崎岖表明read lenngth。

下图是weaker signal的一个示例。这里Pol2的抗体不是很有用，有涣散的峰。穿插相关图中有两个峰，一个是真实的峰位移（185-200bp），别的一个在read length。信号弱的数据会集read length的峰会成为主峰。

一个失利的试验发生的穿插相关图类似于input，在fragment length处很少或没有峰，在read length处有信号十分强的。这种现象的原因或许是在结合位点附近fragments没有明显集合。

穿插相关性质量评价衡量值

穿插相关谱图能够核算评价ChIP_Seq试验信噪比的衡量值，并且试验设计保证fragment length精确。低信噪比和不精确的fragment length 标明ChIP-Seq试验或许有问题。
Normalized strand cross-correlation coefficent (NSC)：
NSC是最大穿插相关值除以布景穿插相关的比率(所有或许的链转移的最小穿插相关值)。NSC值越大标明富集效果越好，NSC值低于1.1
标明较弱的富集，小于1表明无富集。NSC值略微低于1.05，有较低的信噪比或很少的峰，这肯能是生物学真实现象，比方有的因子在特定组织类型中只要很少的结合位点；也或许确实是数据质量差。
Relative strand cross-correlation coefficient (RSC):
RSC是片段长度相关值减去布景相关值除以phantom-peak相关值减去布景相关值。RSC的最小值或许是0，表明无信号；富集好的试验RSC值大于1；低于1表明质量低。

phantompeakqualtools

phantompeakqualtools 是一个用于核算ChIP-Seq数据富集和质量衡量值的一个工具包。咱们将运用该包来核算根据链穿插相关峰的主要插入巨细(fragment length)和根据相对phantom peak的数据质量衡量值。phantompeakqualtools是一个R包，依靠samtools。
下载phantompeakqualtools

wget https://storage.googleapis.com/google-code-archive-downloads/v2/code.google.com/phantompeakqualtools/ccQualityControl.v.1.1.tar.gz
tar -xzf ccQualityControl.v.1.1.tar.gz
cd phantompeakqualtools
# 查看README
less README.txt

Linux下装置

R 
install.packages(\"caTools\", lib=\"~/R/library\")

运行phantompeakqualtools

mkdir -p logs qual

for bam in bam_dir/sample1.final.bam bam_dir/sample2.final.bam
do
bam2=`basename $bam .final.bam`
Rscript run_spp_nodups.R -c=$bam -savp -out=qual/${bam2}.qual > logs/${bam2}.Rout
done

参数含义：

-c: 比对过滤后的bam文件的全路径和名字
-savp:保存穿插相关图
-out：会发生数据集重要特征值的输出文件

输出文件解读
输出文件会发生一个tab切割的名为qual的文件,包括的信息如下：

COL1：Filename：比对过滤的bam文件名
COL2：numReads ：有用的测序深度
COL3：estFragLen：逗号分隔的穿插相关峰以相关性递减顺序排列的值
COL4: corr_estFragLen: 逗号分隔的以递减顺序排列穿插相关值
COL5: phantomPeak: Read length/phantom peak链位移
COL6: corr_phantomPeak: phantom peak相关值
COL7: argmin_corr:穿插相关最小的链位移
COL8: min_corr:穿插相关最小值
COL9: Normalized strand cross-correlation coefficient (NSC) = COL4 / COL8
COL10: Relative strand cross-correlation coefficient (RSC) = (COL4 - COL8) / (COL6 - COL8)
COL11: QualityTag: Quality tag based on thresholded RSC (codes: -2:veryLow,-1:Low,0:Medium,1:High,2:veryHigh)
咱们最重视的值是第9列和第11列。
Cross-correlation plots
课程中的示例数据Nanog_rep1的穿插相关图