提到检测CNV的软件，当下最热门的应该是CNVkit。关于其基本使用如下：
cnvkit.py batch *Tumor.bam -n *Normal.bam -t my_baits.bed -f hg19.fasta \
--annotate refFlat.txt --access data/access-5kb-mappable.hg19.bed \
--output-reference my_flat_reference.cnn -d example3/
batch串起来后产生文件挺多，这篇文章介绍了其一步一步运行的作用，对理解整个流程很有帮助。归根结底，cnv的检测都是基于log2 ratio的计算，cnvkit官网说是做了一些矫正。具体怎么一步步算出来，我也不知道。（找遍全网也没有啊，除非去扒源代码）
log radio= log2( A/B)，在cnv检测里，A是突变样本的指定区域的深度，B是正常样本的指定区域的深度。
网上有一个关于log radio很浅显的解释。举个例子，1号染色体某个区域，深度由100变化到10000 ，深度由100 变化到1如果用倍数关系描述其变化程度，分别为100倍和0.01倍。由少变多时，这个数字接近无穷大；由多变少时，这个数字只会无限接近于0。而用log2 radio表示，则是log2(10000/100)和log2(1/100)，即6.6438与-6.6438。
cnvkit的batch运行文件为xxx.call.cns,格式如
chromosome start end gene log2 cn depth p_ttest probes weight
chr1 227917 76342690 TNFRSF14,RPL22,KIF1B,MTOR,MTHFR,TNFRSF8,SPEN,EPHA2,SZRD1,SDHB,CDA,CDC42,C1QA,C1QC,C1QB,EPHB2,CD52,ARID1A,LCK,HDAC1,PSMB2,MYCL,MPL,PTCH2,MUTYH,RAD54L,CYP4B1,CYP4A11,CDKN2C,JUN,JAK1,MSH4 -0.0104248 2 723.668 0.0132858 912 852.077
chr1 76343828 76365396 MSH4 1.863 8 6013.96 2.90539e-06 8 7.56257
....
chr20 57428315 57486252 GNAS 1.88257 8 5905.01 1.25501e-14 21 19.6902
....
cnv用bed格式存储可以方便后续分析
cnvkit.py export bed xxx.call.cns -o xxx.bed --label-genes
--label-genes选项将会展示基因名，但是只能在输入一个文件时使用，因为默认是显示文件前缀。
过滤掉cn为2的，即可得到所需变异。
chr1 76343828 76365396 MSH4 8
chr16 60500 397030 AXIN1 7
chr20 57428315 57486252 GNAS 8
到这里，还需要注释出外显子区域，因此，需要下载gtf文件
下载后，提取每个基因外显子的坐标文件
zcat hg19.refGene.gtf.gz|awk '($3=="exon"){print$0} ' |awk -v OFS='\t' '{print $1,$4,$5,$12,$10"_"$14}' |sed 's/"//g'|sed 's/;//g' >exon.bed
每列为pos/start/end/transcriptid/gene_exonpos
chr1 11869 12227 NR_148357 LOC102725121_1
chr1 12613 12721 NR_148357 LOC102725121_2
chr1 13221 14362 NR_148357 LOC102725121_3
chr1 11874 12227 NR_046018 DDX11L1_1
chr1 12613 12721 NR_046018 DDX11L1_2
chr1 13221 14409 NR_046018 DDX11L1_3
使用bedtools取交集
bedtools -a xxx.bed -b exon.bed -wa -wb >ann.bed
此时会出现两个问题:
1.转录本不同，对每个基因外显子区域的划分不同
2.在chr16:60500-397030区域，远不止cnvkit已经给出注释的AXIN1这一个基因
问题一中，转录本NR，NM的含义可以看ncbi的说明，一般会选择区域最长的或者是在HGNC上查（或许室间质评也应该规定汇报用的参考转录本）。
问题二，查询bed可知,目标区域实际不包括AXIN1前的许多基因，这些基因cnvkit划分的antitarget部分，但是cnv的划分segement时，划在了一起。显然此区域只应该回报AXIN1。
image.png
选择转录本后，用bedtools 分组整理下结果
bedtools -a xxx.bed -b exon.bed -wa -wb|bedtools groupby -i - -g 1-5 -c 10 -o collapse
大功告成，除了和标准结果由些许不同外.......
(cnvkit误我啊，我要改用exomedepth去)
网友评论