1. 参照序列和注释数据
CNVkit要求有两组基础数据才能进行分析
- Fasta格式的参考基因序列
- 基因注释数据库。需要Bed或者RefFlat格式。
以上的数据都要求是无压缩的独立文件。
fasta格式的序列数据应该都有，但是Bed和RefFlat就不一定了。
Bed格式大家都熟悉，但是RefFlat就不一定了，RefFlat格式是这样的。从左到右是染色体名，基因起点位置，终点位置，基因名字。
chr1 1508981 1509154 SSU72
chr1 2407978 2408183 PLCH2
chr1 2409866 2410095 PLCH2
如果像我这样用的是全基因组WGS数据，既没有BED也没有RefFlat，只有一个gff文件。就需要gff文件转RefFlat格式的工具。估计有的同学手里的事gtf文件，至于gtf和gff文件有什么区别，以后有机会在细说，在这里先来解决格式这个问题。
gff文件的话需要两步走。
- gff转换成gtf
需要用到Cufflinks里的函数gffread
gffread in.gff -T -o out.gtf
- gtf转换成refFlat
需要用到一个叫gtfToGenePred的包，Linux系统不知道为什么卡了半天，mac一下就搞定了。
conda install -c bioconda/label/cf201901 ucsc-gtftogenepred
装完以后用就是了。原理就是把复杂的gft文件精简化了一下。
gtfToGenePred -genePredExt -geneNameAsName2 genes.gtf refFlat.tmp.txt
paste <(cut -f 12 refFlat.tmp.txt) <(cut -f 1-10 refFlat.tmp.txt) > refFlat.txt
rm refFlat.tmp.txt
gzip refFlat.txt
这样就得到了一个refFlat.txt格式的注释文件。
2. BAM对比数据
用最常规的BWA对短序列进行mapping, 比对就可以得到BAM文件。
注意: 不要对bam文件进行去重处理，这样可能导致CNVkit无法检测出拷贝数。其实在GATK4等管道里是否应该进行去重就存在着一定的争议，看SNPs没影响，但是看CNV的话肯定不能去重。
3. 建立对比数据(WGS全基因组的情况)
接下来就是使用batch来计算出各个拷贝数的相对值。
cnvkit.py batch Sample1.bam Sample2.bam -n Control1.bam Control2.bam \
-m wgs -f hg19.fasta --annotate refFlat.txt
这样就会根据每个样本生成包含拷贝数具体信息的.cnn文件。
接下来就可以进入正题，开始正式的可视化以及分析了。