前一章，介绍了文件的准备，然后接下来要做的是跑流程的了。由于这些bam比对的是hg19的参考文件，所以使用GRCh37的内容。GRCh37.bed，exclude.hg37.bed，hg37.ncbiRefSeq.gtf

由于上一章给的都是男性样本，缺少女性样本，所以在额外加三个样本，里面包含一个女性样本。

• SRX27331111
• SRX27331110
• SRX27331108

autobin

autobin会快速估算BAM文件中的读取计数或深度，以估算合理的on-target bin和off-target bin大小。
但是对于WGS比较适用，然而对于WES的话，则需要使用target和antitarget。

cnvkit.py autobin -m wgs SRR31975853.bam SRR31975862.bam SRR31975867.bam -b 50000 -g ~/tmp/data/access.hg19.bed

这里的-b是指，每个 bin（小区间） 应该包含的 平均测序碱基数，默认是10万，可以根据需要进行提前设定。

如果是 靶向扩增测序那么可以这么写：

cnvkit.py autobin -m amplicon SRR31975853.bam SRR31975862.bam SRR31975867.bam -t ~/tmp/data/target.bed 

target和antitarget

对于WES来说，我们需要捕获外显子，所以需要生成特定的target.bed和antitarget.bed

cnvkit.py target ~/tmp/data/hg37.bed --annotate hg37.ncbiRefSeq.gtf --split --short-names -o ~/tmp/data/hg37.target.bed
cnvkit.py access ~/tmp/data/hg37.fa -o ~/tmp/data/access.hg37.bed
cnvkit.py antitarget ~/tmp/data/hg37.target.bed -g ~/tmp/data/access.hg37.bed -o ~/tmp/data/hg37.antitarget.bed

这里需要使用--annotate将用相应的基因名称标记每个区域。也是在ucsc上进行下载就行。这一步是为了提前做准备用的，并不需要多次生成，因此可以相当于用做一个固定的文件。

coverage

上一步生成的hg37.target.bed和hg37.antitarget.bed是我们接下来的依据，对于后续的CNV分析至关重要。

cnvkit.py coverage ~/tmp/cnv/batch/SRR31975862.bam ~/tmp/data/hg37.target.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975862.targetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975862.bam ~/tmp/data/hg37.antitarget.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975862.antitargetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975853.bam ~/tmp/data/hg37.target.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975853.targetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975853.bam ~/tmp/data/hg37.antitarget.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975853.antitargetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975867.bam ~/tmp/data/hg37.target.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975867.targetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975867.bam ~/tmp/data/hg37.antitarget.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975867.antitargetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975851.bam ~/tmp/data/hg37.target.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975851.targetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975851.bam ~/tmp/data/hg37.antitarget.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975851.antitargetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975852.bam ~/tmp/data/hg37.target.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975852.targetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975852.bam ~/tmp/data/hg37.antitarget.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975852.antitargetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975854.bam ~/tmp/data/hg37.target.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975854.targetcoverage.cnn
cnvkit.py coverage ~/tmp/cnv/batch/SRR31975854.bam ~/tmp/data/hg37.antitarget.bed -q 30 -p 4 -o ~/tmp/cnv/batch/cnn/SRR31975854.antitargetcoverage.cnn

这一步是 CNVkit 中计算样本覆盖度的步骤，目的是生成两个文件：target 区域的覆盖度 和 antitarget 区域的覆盖度。这些覆盖度文件将用于后续的 CNV 检测，帮助识别样本中是否存在拷贝数变异（CNV）。

• target coverage 区域 包含了基因、exons 或通过靶向捕获富集的区域，也就是target.bed区域，计算这些区域的覆盖度有助于评估这些区域的拷贝数变异（CNV）。
• antitarget coverage 区域 用作背景对照，以便估计目标区域的覆盖度变化是否可能由于测序深度、GC含量等因素的偏差所引起，而非真实的拷贝数变异。

gender

性别鉴定，在后续的cnv的鉴定中，需要根据性别的不同，采取不同的策略。在cnvkit.py流程中，会单独列出一个sex模块，来进行性别的鉴定。

cnvkit.py sex /root/tmp/cnv/batch/cnn/SRR31975862.targetcoverage.cnn -o ~/tmp/cnv/batch/gender/SRR31975862.gender.txt
cnvkit.py sex /root/tmp/cnv/batch/cnn/SRR31975853.targetcoverage.cnn -o ~/tmp/cnv/batch/gender/SRR31975853.gender.txt
cnvkit.py sex /root/tmp/cnv/batch/cnn/SRR31975867.targetcoverage.cnn -o ~/tmp/cnv/batch/gender/SRR31975867.gender.txt
cnvkit.py sex /root/tmp/cnv/batch/cnn/SRR31975851.targetcoverage.cnn -o ~/tmp/cnv/batch/gender/SRR31975851.gender.txt
cnvkit.py sex /root/tmp/cnv/batch/cnn/SRR31975852.targetcoverage.cnn -o ~/tmp/cnv/batch/gender/SRR31975852.gender.txt
cnvkit.py sex /root/tmp/cnv/batch/cnn/SRR31975854.targetcoverage.cnn -o ~/tmp/cnv/batch/gender/SRR31975854.gender.txt

除了cnvkit.py提供的sex性别检测之外，还可以使用程序判断cnn文件中的SRY基因的深度，SRY基因仅在男性样本中存在深度，女性样本则不存在深度，后续再补充。

reference

对于同一批次的样本，需要建立这一批次的样本的基准线（因为测序仪器，人为操作，样本制备等因素的影响，会导致不同批次的深度，所以要建立改批次的基准线，以降低这些因素的影响，主要是时间上的不同，有玄学在。），这个基准线就是用来判断是否有上调或者是下调等内容。不过，这里有区别，一个是多个样本，一个是单个样本的区别。

多个样本

多个样本时，可以计算平均的基准线。

cnvkit.py reference ~/tmp/cnv/batch/cnn/*coverage.cnn -f ~/tmp/data/hg37.fa -o ~/tmp/cnv/batch/ref/Reference.cnn

单个样本

单个样本时，由于基准线只有他自身，因此要用GRCh37.fa参考基因组建立基准线。

cnvkit.py reference -o FlatReference.cnn -f hg37.fa -t targets.bed -a antitargets.bed

fix

修正 和 标准化 拷贝数数据将使用来自参考基因组的内容对原始的覆盖度数据进行修正，以去除由于 区域覆盖度 和 GC含量 等因素引起的偏差，从而得到更准确的拷贝数变异结果。

cnvkit.py fix /root/tmp/cnv/batch/cnn/SRR31975862.targetcoverage.cnn /root/tmp/cnv/batch/cnn/SRR31975862.antitargetcoverage.cnn ~/tmp/cnv/batch/ref/Reference.cnn -o ~/tmp/cnv/batch/cnr/SRR31975862.cnr
cnvkit.py fix /root/tmp/cnv/batch/cnn/SRR31975853.targetcoverage.cnn /root/tmp/cnv/batch/cnn/SRR31975853.antitargetcoverage.cnn ~/tmp/cnv/batch/ref/Reference.cnn -o ~/tmp/cnv/batch/cnr/SRR31975853.cnr
cnvkit.py fix /root/tmp/cnv/batch/cnn/SRR31975867.targetcoverage.cnn /root/tmp/cnv/batch/cnn/SRR31975867.antitargetcoverage.cnn ~/tmp/cnv/batch/ref/Reference.cnn -o ~/tmp/cnv/batch/cnr/SRR31975867.cnr
cnvkit.py fix /root/tmp/cnv/batch/cnn/SRR31975851.targetcoverage.cnn /root/tmp/cnv/batch/cnn/SRR31975851.antitargetcoverage.cnn ~/tmp/cnv/batch/ref/Reference.cnn -o ~/tmp/cnv/batch/cnr/SRR31975851.cnr
cnvkit.py fix /root/tmp/cnv/batch/cnn/SRR31975852.targetcoverage.cnn /root/tmp/cnv/batch/cnn/SRR31975852.antitargetcoverage.cnn ~/tmp/cnv/batch/ref/Reference.cnn -o ~/tmp/cnv/batch/cnr/SRR31975852.cnr
cnvkit.py fix /root/tmp/cnv/batch/cnn/SRR31975854.targetcoverage.cnn /root/tmp/cnv/batch/cnn/SRR31975854.antitargetcoverage.cnn ~/tmp/cnv/batch/ref/Reference.cnn -o ~/tmp/cnv/batch/cnr/SRR31975854.cnr

由于方便文件管理，我在这一步之前，手动调整了一些文件路径。

segment

为了更加精确的识别CNV，从拷贝数比率表，推断 置信区间，为call变异做准备。

cnvkit.py segment /root/tmp/cnv/batch/cnr/SRR31975862.cnr -p 4 -m cbs --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975862.cns.tmp
cnvkit.py segment /root/tmp/cnv/batch/cnr/SRR31975853.cnr -p 4 -m cbs --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975853.cns.tmp
cnvkit.py segment /root/tmp/cnv/batch/cnr/SRR31975867.cnr -p 4 -m cbs --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975867.cns.tmp
cnvkit.py segment /root/tmp/cnv/batch/cnr/SRR31975851.cnr -p 4 -m cbs --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975851.cns.tmp
cnvkit.py segment /root/tmp/cnv/batch/cnr/SRR31975852.cnr -p 4 -m cbs --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975852.cns.tmp
cnvkit.py segment /root/tmp/cnv/batch/cnr/SRR31975854.cnr -p 4 -m cbs --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975854.cns.tmp

cnvkit.py segmetrics /root/tmp/cnv/batch/cnr/SRR31975862.cnr -s /root/tmp/cnv/batch/cns/SRR31975862.cns.tmp --ci --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975862.cns
cnvkit.py segmetrics /root/tmp/cnv/batch/cnr/SRR31975853.cnr -s /root/tmp/cnv/batch/cns/SRR31975853.cns.tmp --ci --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975853.cns
cnvkit.py segmetrics /root/tmp/cnv/batch/cnr/SRR31975867.cnr -s /root/tmp/cnv/batch/cns/SRR31975867.cns.tmp --ci --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975867.cns
cnvkit.py segmetrics /root/tmp/cnv/batch/cnr/SRR31975851.cnr -s /root/tmp/cnv/batch/cns/SRR31975851.cns.tmp --ci --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975851.cns
cnvkit.py segmetrics /root/tmp/cnv/batch/cnr/SRR31975852.cnr -s /root/tmp/cnv/batch/cns/SRR31975852.cns.tmp --ci --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975852.cns
cnvkit.py segmetrics /root/tmp/cnv/batch/cnr/SRR31975854.cnr -s /root/tmp/cnv/batch/cns/SRR31975854.cns.tmp --ci --drop-low-coverage -o /root/tmp/cnv/batch/cns/SRR31975854.cns

这里还有一些额外的参数，例如-t和-m

• --threshold/-t：指定分段方法的显著性阈值，较低的阈值（如 -t 0.01）表示更高的灵敏度，更多的变动会被认为是显著的，但也可能引入一些 假阳性。
• -m：分割方法。
  • cbs：圆形二进制分段（CBS），常用与外显子测序（WES），默认使用这个。
  • hmm：马尔可夫模型，适用于大多数样本，速度比 cbs 快，当数据质量较好时，能够在较高精度上处理较复杂的分段。
  • hmm-tumor：马尔可夫模型——肿瘤版，能够识别较细的拷贝数变化，
  • hmm-germline：马尔可夫模型——常染色体版，用于常染色体样本。
• --vcf 如果在这一步的时候，添加一个snp的vcf文件，也就是碱基变异的文件的话，cnvkit.py会在vcf的基础上进一步验证cnv变异的真实性。
• --drop-low-coverage：如果是肿瘤样本，则可以添加这个参数。过滤低质量数据。
• --ci：计算置信区间。

call

call变异

cnvkit.py call /root/tmp/cnv/batch/cns/SRR31975862.cns -t=-1.1,-0.4,0.3,0.7 --filter cn --male-reference -o /root/tmp/cnv/batch/cnv/SRR31975862.call.precise.cns
cnvkit.py call /root/tmp/cnv/batch/cns/SRR31975853.cns -t=-1.1,-0.4,0.3,0.7 --filter cn --male-reference -o /root/tmp/cnv/batch/cnv/SRR31975853.call.precise.cns
cnvkit.py call /root/tmp/cnv/batch/cns/SRR31975867.cns -t=-1.1,-0.4,0.3,0.7 --filter cn --male-reference -o /root/tmp/cnv/batch/cnv/SRR31975867.call.precise.cns
cnvkit.py call /root/tmp/cnv/batch/cns/SRR31975851.cns -t=-1.1,-0.4,0.3,0.7 --filter cn --male-reference -o /root/tmp/cnv/batch/cnv/SRR31975851.call.precise.cns
cnvkit.py call /root/tmp/cnv/batch/cns/SRR31975852.cns -t=-1.1,-0.4,0.3,0.7 --filter cn --male-reference -o /root/tmp/cnv/batch/cnv/SRR31975852.call.precise.cns
cnvkit.py call /root/tmp/cnv/batch/cns/SRR31975854.cns -t=-1.1,-0.4,0.3,0.7 --filter cn --male-reference -o /root/tmp/cnv/batch/cnv/SRR31975854.call.precise.cns

• -t=-1.1,-0.4,0.3,0.7指定拷贝数（copy number）分类的硬阈值，用逗号分隔的数值。
  • -1.1 表示拷贝数小于 1.1 的区域会被标记为 缺失（Deletion）。
  • -0.4 表示拷贝数小于 0.4 的区域会被视为 低拷贝数。
  • 0.3 表示拷贝数大于 0.3 的区域会被视为 正常（Diploid）。
  • 0.7 表示拷贝数大于 0.7 的区域会被视为 扩增（Amplification）。

导出vcf

在后期的注释中，不是使用cns文件，而是用到vcf文件，所以要用cnvkit.py导出vcf文件

cnvkit.py export vcf /root/tmp/cnv/batch/cnv/SRR31975862.call.precise.cns --male-reference -i SRR31975862 -o /root/tmp/cnv/batch/vcf/SRR31975862.call.precise.cns.vcf
cnvkit.py export vcf /root/tmp/cnv/batch/cnv/SRR31975853.call.precise.cns --male-reference -i SRR31975853 -o /root/tmp/cnv/batch/vcf/SRR31975853.call.precise.cns.vcf
cnvkit.py export vcf /root/tmp/cnv/batch/cnv/SRR31975867.call.precise.cns --male-reference -i SRR31975867 -o /root/tmp/cnv/batch/vcf/SRR31975867.call.precise.cns.vcf
cnvkit.py export vcf /root/tmp/cnv/batch/cnv/SRR31975851.call.precise.cns --male-reference -i SRR31975851 -o /root/tmp/cnv/batch/vcf/SRR31975851.call.precise.cns.vcf
cnvkit.py export vcf /root/tmp/cnv/batch/cnv/SRR31975852.call.precise.cns --male-reference -i SRR31975852 -o /root/tmp/cnv/batch/vcf/SRR31975852.call.precise.cns.vcf
cnvkit.py export vcf /root/tmp/cnv/batch/cnv/SRR31975854.call.precise.cns --male-reference -i SRR31975854 -o /root/tmp/cnv/batch/vcf/SRR31975854.call.precise.cns.vcf