队列研究GATK VQSR 的高速执行方式

    变异质控过滤，是下游分析非常重要的一个环节 。通常人类重测序和人类队列分析我们使用GATK VQSR模块来完成这个事情。对人少数样品的重测序，VQSR很容易完成，对于大规模的队列分析就艰难许多。 常常会碰到内存不足，处理缓慢等很多情况。BinBash 科学计算平台带您一起快速解决大规模队列分析的VQSR 过滤。

先简单回顾下VQSR的命令行，其中参数意义，本文不多做讨论。

gatk VariantRecalibrator \
   -R $reference \
   -V $input \
   --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.hg38.vcf \
   --resource:omini,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.hg38.vcf \
   --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg38.vcf \
   --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_138.hg38.vcf \
   -an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
   -mode SNP \
   -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
   --tranches-file snps.tranches \
   -O snps.recal.gz 

gatk ApplyVQSR \
   -R $reference \
   -V $input \
   --ts_filter_level 99.0 \
   --tranches-file snps.tranches \
   --recal-file snps.recal.gz \
   -mode SNP \
   -O snps.VQSR.vcf.gz

gatk VariantRecalibrator \
   -R $reference  \
   -V snps.VQSR.vcf.gz \
   --resource:mills,known=true,training=true,truth=true,prior=12.0 Mills_and_1000G_gold_standard.indels.hg38.vcf \
   --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_146.hg38.vcf \
   -an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
   -mode INDEL \
   --tranches-file indels.tranches \
   -O indels.recal.gz 

gatk ApplyVQSR \
   -R $reference \
   -V snps.VQSR.vcf.gz \
   --ts_filter_level 99.0 \
   --tranches-file indels.tranches \
   --recal-file indels.recal.gz \
   -mode INDEL \
   -O VQSR.vcf.gz 

    大规模的队列分析，往往产生非常大的VCF文件，比如BinBash 平台完成的2000例人的VCF文件大约5TB，100列的人大约650GB，3400例的豆数据大约8TB，2900例棉花的大约5TB， 1000例猴子的大约4TB......。 大规模队列研究产生的如此大规模的数据，直接使用GATK 做VQSR 会非常缓慢。因此，根据VQSR所必须使用的信息，对原始VCF文件进行处理，会带来巨大的时间节省和计算资源节省效益。

    按照VQSR 算法，仅仅使用了VCF的前9列信息，但是为了保证GATK对VCF文件格式要求，因此，对大队列分析数据，提取出前10列数据，形成一个新的包含所有突变位点的小的VCF文件，以供完成VariantRecalibrator  步骤足以。一行简单的脚本即可得到前10列信息：

zcat $input | awk -v OFS="\t" '{if($0~"##"){print $0}else{print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10}}' |bgzip -c > main.vcf.gz 
bcftools index -t main.vcf.gz 

    然后，使用main.vcf.gz  文件作为VariantRecalibrator  输入文件，然后使用原始的vcf 文件作为ApplyVQSR 的输入文件即可很快的完成VQSR 步骤了。

                    

公众号                                          微信