使用cnvkit来对大批量wes样本找cnv
cnvkit被设计来处理同一个批次的多个肿瘤配对样本测序情况,首先对所有的normal数据进行bin处理拿到背景值,然后就这个背景值来处理所有的tumor测序数据计算拷贝数变异情况。
该软件使用比较复杂,建议读一读官网教程。所有的命令都被包装到一个python脚本里面,使用该脚本调用一系列字命令,如下:
- Copy number calling pipeline
每个命令都有自己的特殊功能,需要仔细阅读。
流程图:
流程代码如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | cnvkit.py access baits.bed --fasta hg19.fa -o access.hg19.bed cnvkit.py autobin *.bam -t baits.bed -g access.hg19.bed [--annotate refFlat.txt --short-names] # For each sample... cnvkit.py coverage Sample.bam baits.target.bed -o Sample.targetcoverage.cnn cnvkit.py coverage Sample.bam baits.antitarget.bed -o Sample.antitargetcoverage.cnn # With all normal samples... cnvkit.py reference *Normal.{,anti}targetcoverage.cnn --fasta hg19.fa [--male-reference] -o my_reference.cnn # For each tumor sample... cnvkit.py fix Sample.targetcoverage.cnn Sample.antitargetcoverage.cnn my_reference.cnn -o Sample.cnr cnvkit.py segment Sample.cnr -o Sample.cns # Optionally, with --scatter and --diagram cnvkit.py scatter Sample.cnr -s Sample.cns -o Sample-scatter.pdf cnvkit.py diagram Sample.cnr -s Sample.cns [--male-reference] -o Sample-diagram.pdf |
可以看到软件提供的命令基本上的都用到了,coverage—>fix—>segment—>segment
其实上面这么一大串的命令已经被包装成了一句命令,就是:
1 | cnvkit.py batch *.bam -r my_reference.cnn -p 8 |
这个一句话命令与上面的多行代码是等效的,默认的segment算法是 circular binary segmentation algorithm (CBS),也可以用-m
切换使用其它算法,比如: faster HaarSeg (haar
) or Fused Lasso (flasso
)
上面得到的只是segment的结果,还可以call一下:
1 2 3 4 | cnvkit.py call Sample.cns -o Sample.call.cns cnvkit.py call Sample.cns -y -m threshold -t=-1.1,-0.4,0.3,0.7 -o Sample.call.cns cnvkit.py call Sample.cns -y -m clonal --purity 0.65 -o Sample.call.cns cnvkit.py call Sample.cns -y -v Sample.vcf -m clonal --purity 0.7 -o Sample.call.cns |
这个时候需要考虑到已有的vcf变异文件,或者计算好的tumor纯度,或者倍性等等。把segment计算得到的log2 ratio值还原成 0,1,2,3,4这样的拷贝数。
但是,事实上上面的代码一般来说是不能直接使用的,因为我们的测序数据通常是WES数据,需要加上很多参数。
实践运行cnvkit
上面流程很复杂,命令也很多,但是不知道也没关系,用起来其实就一个batch命令即可,当然这个batch命令本身参数也不少,而且被设计用来处理不同的数据情况。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | # From baits and tumor/normal BAMs ## 同一批次的所有样本N/T测序数据的bam文件一起运行 cnvkit.py batch *Tumor.bam --normal *Normal.bam \ --targets my_baits.bed --annotate refFlat.txt \ --fasta hg19.fasta --access data/access-5kb-mappable.hg19.bed \ --output-reference my_reference.cnn --output-dir results/ \ --diagram --scatter ## 如果新增加了肿瘤测序数据,就运行下面的 # Reusing a reference for additional samples cnvkit.py batch *Tumor.bam -r Reference.cnn -d results/ # Reusing targets and antitargets to build a new reference, but no analysis cnvkit.py batch -n *Normal.bam --output-reference new_reference.cnn \ -t my_targets.bed -a my_antitargets.bed --male-reference \ -f hg19.fasta -g data/access-5kb-mappable.hg19.bed |
值得注意的就是,如果是全基因组测序数据,用batch --method wgs
,如果是捕获基因组测序,包括全外显子,就用 batch --method amplicon
,然后一定要提供捕获区域的bed文件,一般是外显子加上其侧翼上下游的50bp长度。
人类外显子长度平均是200bp,所以默认的bin是267bp,这样可以把比较长的exon给拆分开来。
还有就是 access
参数需要的文件,
输入输出文件
其中 coverage 命令会对每一个 normal 样本都计算 *.targetcoverage.cnn
and *.antitargetcoverage.cnn
files , 说明是: target and antitarget coverage tables (.cnn)
这些文件需要合并起来:
1 | cnvkit.py reference *coverage.cnn -f ucsc.hg19.fa -o Reference.cnn |
然后再校正区域测序深度及GC含量,之后变成 copy number ratios (.cnr) 文件。
1 | cnvkit.py fix Sample.targetcoverage.cnn Sample.antitargetcoverage.cnn Reference.cnn -o Sample.cnr |
最后把 copy number ratios (.cnr) 文件用segment算法跑一下即可,输出cns后缀文件的 segment信息。
对于normal样本只需要输出cnn即可,合并成Reference.cnn
,然后把一个个tumor样品,根据这个Reference.cnn
来计算 cnr,进而计算 cns 。
最后对于cns可以进行call找到真正的拷贝数。
1 2 3 4 5 6 7 8 9 10 11 | $ headNPC_merge_marked_fixed.cns |cut -f 1-3,5,8 chromosomestart end log2weight chr112098 1701806 -0.183469 84.0794 chr11702902 1752401 -0.962192 4.87216 chr11752901 12777601-0.220165 370.756 chr11277810112920301-1.1168810.7699 chr11292030727407686-0.275998 558.214 chr127408186125184087 -0.04474042418.53 chr1143185087 248945922 -0.04226292967.61 chr210500 5692985 0.15197885.3751 chr25692985 90402011-0.03291651874.56 |
上面的segment结果还可以call一下,如果有需要的话。
可视化结果:
很明显可以看到有拷贝数变异的区域了。
相关推荐:
- NCBI在线BLAST使用方法与结果详解 2938
- 神经网络术语:Epoch、Batch Size和迭代 527
- Consed的安装与使用教程 465
- 陈连福的NGS生物信息学培训教材V2.1 277
- WGCNA分析使用教程 272
最新创建圈子
-
原料药研发及国内外注册申报
2019-01-25 10:41圈主:caolianhui 帖子:33 -
制药工程交流
2019-01-25 10:40圈主:polysciences 帖子:30 -
健康管理
2019-01-25 10:40圈主:neuromics 帖子:20 -
发酵技术
2019-01-25 10:39圈主:fitzgerald 帖子:17 -
医学肿瘤学临床试验
2019-01-25 10:39圈主:bma 帖子:58
谢谢提醒。问题已经修复。