前言
- 原理简单介绍
- 工具介绍
2.1 参数介绍
2.2 代码模板 - 其他说明
前言
当我们拿到多个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是失败的，还是要具体情况具体分析。
网友评论