使用CNVkit进行CNV分析

CNVkit是一款CNV预测软件，适用于全外显子，目的区域靶向测序等数据的CNV检测，官网如下

https://cnvkit.readthedocs.io/en/stable/

被很多高分文章引用， 示意如下

文章发表在PLos上，链接如下

https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004873

和Excavator类似，CNVKit也将基因组划分为in-target和off-target两部分，流程图示意如下

将in-target和off-target区域划分成小的bin区间，统计bin区间内的测序深度，综合考虑GC含量，目的区域的大小和分布密度， 重复元件等因素，对原始测序深度进行校正，然后计算相对对照样本的log2 ratio, 通过segmentation算法来划分segment, 支持cbs, haar, flasso多种segmentation算法。

CNVkit采用python进行开发，使用方便，集成了可视化功能，可以直观的展示分析结果，支持导出多种格式的结果文件，可以很好的与下游软件相结合。采用了模块化的开发思想，按照功能拆分成了独立的模块，示意如下

每个功能模块对应一个子命令，为了方便调用，将pipeline的完整功能封装在batch子命令中，只通过这一个命令来运行整个pipeline。整个数据分析的过程可以分为以下几个部分

1. 输入区域文件

对于靶向测序而言，肯定需要输入目的区域的文件，格式为bed， target子命令用于处理目的区域的bed文件，可以添加对应的基因注释等信息，用法如下

cnvkit.py target \
my_baits.bed \
--annotate refFlat.txt \
-o my_targets.bed

除了目的区域in-target，还需要计算off-set区域，也称之为antitarget, in-target和off-target区域加起来就是基因组上所有可覆盖的区域。二测测序并不能达到100%的覆盖度，基因组上的高重复区域，端粒，着丝粒等区域就无法覆盖，所以cnvkit通过access子命令来计算基因组上可以覆盖到的区域，命令如下

cnvkit.py access \
hg19.fa \
-x excludes.bed \
-o access.hg19.bed

计算出可覆盖的区域之后，减去in-target区域， 就可以得到off-target区域，通过antitarget子命令来实现，代码如下

cnvkit.py antitarget \
my_targets.bed \
-g access.hg19.bed \
-o my_antitargets.bed

2. 计算样本的测序深度

coverage和autobin两个子命令都可以用来计算测序深度，以coverage为例，用法如下

cnvkit.py coverage \
Sample.bam \
my_targets.bed \
-o Sample.targetcoverage.cnn

cnvkit.py coverage \
Sample.bam \
my_antitargets.bed \
-o Sample.antitargetcoverage.cnn

分别统计target和antitarget区域的测序深度信息，输出结果后缀为cnn,是cnvkit中定义的一种格式，专门用来存储测序深度信息。

3. 构建正常基因组的测序分布模型

通过reference子命令来构建正常基因组的测序分布模型，采用对照样本的测序深度，校正GC含量等系统误差。有多个对照样本时，可以将所有的对照样本合并来创建，用法如下

cnvkit.py reference \
*coverage.cnn \
-f hg19.fa \
-o Reference.cnn

当没有对照样本时，软件可以模拟出一个正常的测序深度分布模型，用法如下

cnvkit.py reference \
-o FlatReference.cnn \
-f hg19.fa \
-t my_targets.bed \
-a my_antitargets.bed

4. 计算实验样本相对正常对照的log2 ratio

通过fix子命令计算log2 ratio，用法如下

cnvkit.py fix \
Sample.targetcoverage.cnn \
Sample.antitargetcoverage.cnn \
Reference.cnn \
-o Sample.cnr

输出结果后缀为cnr, 是cnvkit中定义的一种格式，专门用来存储log2ratio的信息。

5. 划分segment, 计算拷贝数

通过segment子命令进行segment的划分，用法如下

cnvkit.py segment \
Sample.cnr \
-o Sample.cns

输出结果后缀为cns, 是cnvkit中定义的一种格式，和SEG格式类似，用来存储CNV分析的结果。接下来还可以通过call子命令，计算每个segment区域的绝对拷贝数，用法如下

cnvkit.py call \
Sample.cns \
-o Sample.call.cns

关于各种文件格式的具体解释可以参考以下链接

https://cnvkit.readthedocs.io/en/stable/fileformats.html

6. 结果可视化

提供了以下3种可视化的子命令

1. diagram
2. scatter
3. heatmap

diagram用于展示单个样本的CNV在染色体上的分布，用法如下

cnvkit.py diagram \
-s Sample.cns \
Sample.cnr

可视化结果如下

scatter子命令展示单个样本染色体区域上log2 ratio值的分布，用法如下

cnvkit.py scatter \
-s Sample.cns \
Sample.cnr

可视化结果示意如下

heatmap子命令展示多个样本的CNV分布情况，用法如下

cnvkit.py heatmap *.cns

可视化结果示意如下

分布运行步骤繁琐，上面的代码可以通过batch一个子命令搞定，用法如下

cnvkit.py batch \
*Tumor.bam \
--normal *Normal.bam \
--targets my_baits.bed \
--annotate refFlat.txt \
--fasta hg19.fasta \
--access access.hg19.bed \
--output-reference my_reference.cnn \
--output-dir results/ \
--diagram --scatter

以上就是cnvkit的基本用法，更多细节请参考官方文档。