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
文件。
接下来就可以进入正题,开始正式的可视化以及分析了。
网友评论