CNVkit使用手册

全基因组CNV分析2. 快速上手

2021-09-07  本文已影响0人  Jason数据分析生信教室

1. 参照序列和注释数据

CNVkit要求有两组基础数据才能进行分析

  1. Fasta格式的参考基因序列
  2. 基因注释数据库。需要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文件的话需要两步走。

  1. gff转换成gtf
    需要用到Cufflinks里的函数gffread
gffread in.gff -T -o out.gtf
  1. 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文件。

接下来就可以进入正题,开始正式的可视化以及分析了。

上一篇下一篇

猜你喜欢

热点阅读