高通量测序数据处理学习记录(七):使用ChIPQC包检查ChIP-seq的质量

参考文章

六六
这篇文章的作者看来是详细的观看了ChIPQC manual并且参考了已有文档进行了学习和实践。
Manual见Manual

讨论

问:为什么有了参考文章还要发文呢?
答:因为参考文章只是重现了作者的vignette而已,本人在实战过程中就遇到了无数的问题,花费了一天的时间才完成整个分析(包括重装了R +-+),还在Biostar&Bioconductor forum上发了两篇post:
Question: ChIPQC BiocParallel errors
Question: ChIPQC ERROR with "subscript out of bounds"

分析重现 → 原版pdf

1. 介绍

ChIPQC是一个用来分析ChIP-seq数据的Bioconductor包,它整合了可视化以及统计分析,用以处理比对完成的bam文件(需要peak文件)。
这个包主要有两个函数:

  • ChIPQCsample()函数负责对单个ChIP-seq数据进行质控生成报告
  • ChIPQC()函数可以对一组数据进行质控并生成报告

ChIPQC需要准备一个信息列表,提供给的ChIPQC()函数构建ChIPQCexperiment object

2. 示例

2.1 准备experiment sample sheet

ChIPQC包自带了准备好的sample信息表,但是未准备bam文件及peak文件,因为bam文件太大了,我们可以先看一下sample文件:

> library(ChIPQC)
> samples = read.csv(file.path(system.file("extdata", package="ChIPQC"), "example_QCexperiment.csv"))
> samples
  SampleID  Tissue Factor Replicate            bamReads                           Peaks
1   CTCF_1    A549   CTCF         1 reads/SRR568129.bam peaks/SRR568129_chr22_peaks.bed
2   CTCF_2    A549   CTCF         2 reads/SRR568130.bam peaks/SRR568130_chr22_peaks.bed
3   cMYC_1    A549   cMYC         1 reads/SRR568131.bam peaks/SRR568131_chr22_peaks.bed
4   cMYC_2    A549   cMYC         2 reads/SRR568132.bam peaks/SRR568132_chr22_peaks.bed
5   E2F1_1 HeLa-S3   E2F1         1 reads/SRR502355.bam peaks/SRR502355_chr22_peaks.bed
6   E2F1_2 HeLa-S3   E2F1         2 reads/SRR502356.bam peaks/SRR502356_chr22_peaks.bed

上图包含了所有必须的信息,bamReads和Peaks的参数是相对路径,当然你也可以指定直接路径,可以看出来这就是一个data.frame,在R里面都可以直接自己构建。

2.2 构建一个ChIPQCexperiment object

exampleExp = ChIPQC(samples, annotation = "hg19", blacklist = NULL, chromosome = "chr18")

这个里面需要注意的是,annotation是指你使用的参考基因组,其中一下物种的参考基因组的指定已经内置:

“hg19”: Human, version 19
“hg18”: Human, version 18
“mm10”: Mouse, version 10
“mm9” : Mouse, version 19
“rn4” : Rat, version 4
“ce6” : C. Elgans, version 6
“dm3” : D. Melanogaster, version 3

其中blacklist位置是需要不加以分析的区域,不理解的请去Google一下就可以,每个物种的blacklist是可以从UCSC上下载得到的。
另外,需要注意的是,chromosome参数也很重要,当你不去指定的时候,只会默认读取第一条序列进行统计。
作者的包里面内置了构建好的ChIPQCexperiment object,使用data(example_QCexperiment)调用:

> data(example_QCexperiment)
> exampleExp
Samples: 6 : CTCF_1 CTCF_2 ... E2F1_1 E2F1_2 
        Tissue Factor Replicate Peaks
CTCF_1    A549   CTCF         1  1118
CTCF_2    A549   CTCF         2   648
cMYC_1    A549   cMYC         1   253
cMYC_2    A549   cMYC         2   159
E2F1_1 HeLa-S3   E2F1         1   325
E2F1_2 HeLa-S3   E2F1         2   249
        Reads Map% Filt%   Dup% ReadL FragL RelCC  SSD  RiP% RiBL%
CTCF_1 341055  100  26.3 16.600    28   131 2.740 2.53 31.20  1.33
CTCF_2 303856  100  28.4  7.320    28   131 2.690 1.43 12.80  0.00
cMYC_1 287462  100  31.0 13.600    28    97 0.347 1.47  6.59  1.78
cMYC_2 317537  100  29.9  4.540    28   129 0.386 1.09  2.79  0.00
E2F1_1 223580  100  31.7  1.010    28   101 0.701 1.29  7.80  2.00
E2F1_2 194919  100  31.8  0.663    28   107 0.303 1.40  5.36  2.79

2.3 质量控制指标摘要

我们来依次看一下exampleExp里面的信息(解读的很粗略,寻求更深理解的可以自行阅读manual):

  • Reads 可以知道bam文件里面有多少reads
  • Map% 可以知道bam文件的比对率如何,如上我们可以知道,bam文件是过滤过的,丢弃了为比对上的reads
  • Filt% 指的是根据比对质量,过滤的比例
  • Dup% 很明显指的是重复率,也是要去除的
  • ReadL 读长
  • FragL 根据比对情况估算出来的fragment长度
  • RelCC 一般大于1,越大越好,大于1可以认为富集的不错
  • SSD 同样是富集得分,对高信号值比较敏感,越高越好
  • Rip% 在peaks里面的reads占所有reads的比率,一般来说大于5%
  • RiBL% 在blacklist region里面的reads占所有reads的比率,可以看出背景的多少,太高了就有问题

2.4 可视化数据

最后的操作很简单:

> ChIPQCreport(exampleExp)
Saving 7 x 7 in image
Saving 7 x 7 in image
Saving 7 x 7 in image
Using Sample as id variables
Saving 7 x 7 in image
Using Sample as id variables
Saving 7 x 7 in image
Saving 7 x 7 in image

异常无脑的操作,在当前工作目录就会出现一个ChIPQCreport文件夹,点开后会有图片以及一个报告html,打开查看即可。


以下为高通量测序数据处理系列快速通道:

高通量测序数据处理学习记录(零):NGS分析如何选择合适的参考基因组和注释文件

高通量测序数据处理学习记录(一):比对软件STAR的使用

高通量测序数据处理学习记录(二):Read Counts的提取

高通量测序数据处理学习记录(三):Pathway Analysis及GSEA

高通量测序数据处理学习记录(四):DeepTools学习笔记

高通量测序数据处理学习记录(五):上传二代测序数据到GEO

高通量测序数据处理学习记录(六):什么是测序深度和测序覆盖度?

高通量测序数据处理学习记录(七):使用ChIPQC包检查ChIP-seq的质量

推荐阅读更多精彩内容