deeptools工具系列——ChIP数据质控富集程度

前言

1. 原理简单介绍
2. 工具介绍
   2.1 参数介绍
   2.2 代码模板
3. 其他说明

前言

当我们拿到多个ChIP-seq数据时，我们会考虑，这么多数据，哪些数据是更好的呢？这些数据中有哪些是可以用于下游分析的呢？

这其实就是关于ChIP-seq数据质控的问题。在2012年ENCODE在Genome Research上发表过一篇文章介绍用于做质控的参数，具体可以参考：

ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia

里面的参数有些依赖于call peak后才能计算，例如FRiP值的计算。这些需要call peak的计算都会受到不同软件call peak结果的影响。
我找到了deeptools工具中的一个质控软件——plotFingerprint

该工具无需call peak即可给出一个确定ChIP-seq数据是否有富集的初始判断。

1. 原理简单介绍

其思想也相对简单：

按照bin的大小，将基因组进行分割后，统计每个bin中有read的数目。

将所有bin按照read数目多少进行排序，由小到大。

以x轴为累积bin的百分比，y轴为累积reads的百分比进行绘制结果。

因为x和y轴都是累积百分比，所以：

对于Input control来说，如果它的信号是均匀随机分布在整个基因组上，那么它绘制出来的结果图应该就是一条对角线。

对于实验组IP数据来说，如果它的信号是不均匀随机分布在整个基因组上，在某些区域有很强的富集程度，那么在x轴前部分的bin中累积read相对较低，而最后一小部分bin中的累积read则会出现陡升。绘制出来的结果图应该是一条类似指数分布曲线的样子。

结果图举例：

image.png

可以看到：

input数据的富集效果最差，而红色的H3K4me3的富集效果最好。

并且我们甚至可以根据曲线下面积进行定量计算。而这个曲线下面积在 plotFingerprint 工具里甚至已经帮我们计算好了！

2. 工具介绍

2.1 参数介绍

usage: plotFingerprint -b treatment.bam control.bam -plot fingerprint.png

必要参数:
  --bamfiles/-b                                指定bam文件
  --plotFile/-plot/-o                        指定输出文件名，可以根据文件后缀自动判断文件格式，后缀可以是{png, eps, pdf, svg}
  --outRawCounts                            指定文件名保存每个bin中read的数目
  --outFileSortedRegions                    如果这里给了参数，那么排序后热图的y轴坐标值就会被保存。默认为None


Read处理可选参数:
  --extendReads/-e                    对于SE数据，根据实验过程中打断的片段大小则直接填写该数值，需要手动指定；对于PE数据可以不指定具体数据，根据BAM文件信息自行计算。PE数据中没有配对的数据以及配对长度超过4*fragment的reads而言，则当作SE数据处理。
  --ignoreDuplicates                起点和终点相同的read当作一条进行处理
  --minMappingQuality                至少达到最低mapping质量得分的read才被考虑
  --centerReads                        相对于片段长度，reads处于中心位置。对于PE数据会自动计算，对于SE数据则仅当设置了--extendReads才会有效
  --samFlagInclude                    根据bam文件Flag进行过滤
  --samFlagExclude                    根据bam文件Flag进行过滤
  --minFragmentLength                主要用于ATACseq
  --maxFragmentLength                主要用于ATACseq

可选参数:
  --labels/-l                        指定样本标签，如果不指定则以文件名进行命名，多个label之间以空格间隔
  --smartLabels                        自动移除路径和扩展名后使用文件名作为文件标签
  --binSize/-bs INT                    手动指定bin的大小，默认为10000
  --numberOfSamples/-n                抽样数目，计算中并不是用全部的bin进行计算，而是从中随机抽样进行计算
  --plotFileFormat                    输出文件格式
  --plotTitle/-T                    输出文件title
  --skipZeros                        是否忽略0，默认不忽略。如果设置则忽略
  --outQualityMetrics                指定文件名输出绘制结果图中的具体统计结果值。具体统计结果值含义参考：https://deeptools.readthedocs.io/en/latest/content/feature/plotFingerprint_QC_metrics.html
  --JSDsample                        指定统计分析所需要的control。具体统计结果值含义参考：https://deeptools.readthedocs.io/en/latest/content/feature/plotFingerprint_QC_metrics.html
  --region/-r CHR:START:END            指定基因组上有限区域，只分析这部分数据，可以节省时间，多用于测试代码
  --blackListFileName/-bl            指定blackList的bed文件路径  
  --numberOfProcessors/-p            指定线程

2.2 代码模板

$ plotFingerprint \
 -b testFiles/*bam \
--labels H3K27me3 H3K4me1 H3K4me3 H3K9me3 input \
--minMappingQuality 30 --skipZeros \
--region 19 --numberOfSamples 50000 \
-T "Fingerprints of different samples"  \
--plotFile fingerprints.png \
--outRawCounts fingerprints.tab \
--outQualityMetrics fingerprints_qc_metrics.tab

$ head fingerprints.tab
  #plotFingerprint --outRawCounts
  'H3K27me3'  'H3K4me1'       'H3K4me3'       'H3K9me3'       'input'
      1       0       0       0       0
      0       0       0       0       1
      0       1       0       0       0
      12      0       0       3       3
      3       0       1       1       0
      6       4       0       1       0
      1       0       0       0       0
      4       1       1       1       0
      1       0       0       0       0

3. 其他说明

image.png

第一张图可以很清晰给我们解释这个图的原理，前97%的bin中只包含了55%的reads，而剩下3%的bin中却包含了剩下约45%的reads，说明了该ChIP-seq的富集程度高。

第二张图则告诉我们，在最开始的一些bin中没有reads。我们可以根据图中前面没有reads的bin所占百分比有一个大概了解。

第二张图告诉我们，有些broad peak的ChIP-seq数据，如果input和IP差异不大，我们也不能说该IP是失败的，还是要具体情况具体分析。