参考文章
六六
这篇文章的作者看来是详细的观看了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分析如何选择合适的参考基因组和注释文件
高通量测序数据处理学习记录（二）：Read Counts的提取
高通量测序数据处理学习记录（三）：Pathway Analysis及GSEA
高通量测序数据处理学习记录（四）：DeepTools学习笔记
