使用原有的BWA+GATK进行SNP calling存在严重的时长慢问题;无论是比对,还是gvcf calling,还是joint calling 都是严重的限速步骤。
可以利用V100/A100等计算卡对这些步骤就行优化;使用的是nvidia的整合包parabricks。
安装自己看官网, GPU没有去借一块或者买一块
https://docs.nvidia.com/clara/parabricks/latest/index.html
对单个样本$sample,进行gvcf calling:
pbrun deepvariant_germline --ref $reference_genome --in-fq 1.fq.gz 2.fq.gz '@RG\tID:$sample\tPL:illumina\tLB:$sample\tSM:$sample\tPU:$sample' --gvcf --out-variants $sample.g.vcf.gz --out-bam $sample.markdup.bam
reference_genome 是参考基因组。
注意下,deepvariant的下游需要用GLnexus进行joint calling,而不是GATK(GVCF缺少GATK需要的信息、)
glnexus_cli -m 120 --config DeepVariant *.vcf.gz > $sample.bcf
然后对BCF进行过滤就行,注意该BCF文件的过滤需要调整参数,推荐参数:
先用bcftools筛选SNPS, 过滤基因型和次等位基因频率
bcftools view --types snps -m 2 -M 2 -q 0.05:minor $sample.bcf > filter1.vcf
再用VCFtools进一步过滤质量:
vcftools --remove-indels --recode --recode-INFO-all --vcf filter1.vcf--maf 0.05 --minDP 5 --maxDP 50 --minGQ 20 --hwe 0.001 --max-missing 0.8 --stdout > filter2.vcf
bgzip filter2.vcf
tabix filter2.vcf.gz
最后过滤一下half call的
bcftools filter -e 'GT="./0" || GT="./1" || GT="./2"' filter2.vcf.gz -o filter3.vcf
这样我们就过滤完了,可以用于下游的各种分析