基于bcftools + vcftools的vcf文件快速过滤流程

目录

群体变异筛选的目的

        群体变异VCF文件即包含一个群体内所有变异信息的文件，通常来自于bwa + gatk等流程获得的gvcf文件Combine而得。同时，由gatk等软件Combine而得的vcf文件往往包含很多低质量的变异，此处的低质量有两个层面：1.在变异位点层面，已有测序的reads在深度和QUAL值等方面不足以证明该位点的存在；2.在群体层面，如果一个变异位点的MAF很低或缺失频率很高，保留该位点可能对后续群体分析造成影响。

变异筛选条件以及软件选择

        在做变异筛选条件之前，推荐详细了解VCF的文件格式

        我第一次接触的群体分析时，组内老师给了我一套现成的分析流程:bwa + gatk来mapping reads和calling变异信息，后续gvcf文件的combine以及筛选均通过gatk的模块进行，参数如下：

gatk --java-options -Xmx20G VariantFiltration -R <ref.genome> -V SNP.vcf --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 30.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 " --filter-name "SF" --genotype-filter-expression "GQ < 30 || DP < 15 || DP > 100" --genotype-filter-name "GF" -O SNP.filter.vcf

gatk --java-options -Xmx20G VariantFiltration -R <ref.genome> -V INDEL.vcf --filter-expression "QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0" --filter-name "my_indel_filter" --genotype-filter-expression "GQ < 30" --genotype-filter-name "GF" -O INDEL.filter.vcf

        熟悉gatk的同学应该都能体会在大容量群体中，gatk从gvcf combine、genotype到variantfiltration这几步，每一步都需要消耗大量的时间,尽管sentieon加速软件很大程度上解决了combine这一主要耗时步骤的问题。但作为一个每周都需要参加组会汇报的苦逼研究生，手握几百份群体数据，每次调试参数运行下来消耗数天时间，这样下去导师那一关肯定不好过了。

        因此，我在这个过程中探索了其它加速的方式。

        首先，将基因组变异按照染色体拆分是大家普遍的选择，vcftools可以很好的做到将gvcf文件按照染色体id拆分(均可通过conda命令安装):

vcftools --gzvcf in.vcf.gz --chr <chr_id> --recode --recode-INFO-all --stdout | bgzip -f > chr_id.vcf.gz

tabix -p vcf chr_id.vcf.gz ##为vcf文件构建索引

        群体变异Combine这个步骤，我是通过sentieon的减速实现的，对于染色体长度60-90Mb的基因组，约900份平均5-10×重测序数据，我每条染色体使用十个线程用了6.5-11.3个小时跑完。

        bcftools将vcf文件拆分为SNP.vcf 和 InDel.vcf:

bcftools view <all.vcf.gz> --threads <n> -v snps -O z -o SNP.vcf.gz
bcftools view <all.vcf.gz> --threads <n> -v indels -O z -o InDel.vcf.gz

        其实不难发现，bcftools筛选的SNP和InDel变异存在一些问题，当位点同时具有SNP和InDel变异的信息时，软件会将两个结果都保留下来，即SNP.vcf和InDel.vcf文件中有着相同的位点信息。如果该问题影响后续分析可以自己写脚本处理一下，或者按照流程稍后去除。

        bcftools分别对SNP和InDel的vcf文件进行筛选:

bcftools filter --threads <n> -s Filter -e \'QD < 1.5 || DP < 5 || FS > 60.0 || MQ < 30.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0\' -s VeryLowQual -e \'QUAL < 30.0\' -s LowQual -e \'QUAL >= 30.0 && QUAL < 50.0\' -O z -o SNP_filter.vcf.gz SNP.vcf.gz
bcftools filter --threads <n> -s Filter -e \'QD < 1.5 || DP < 5 || FS > 200.0 || ReadPosRankSum < -20.0\' -s VeryLowQual -e \'QUAL < 30.0\' -s LowQual -e \'QUAL >= 30.0 && QUAL < 50.0\' --threads 10 -O z -o InDel_filter.vcf.gz InDel.vcf.gz

        其中一些参数和上面的gatk VariantFiltration有所区别，类似的参数都可以根据形式调整，通过筛选条件的变异第七列会有“PASS”的标记，或在满足条件情况下赋予"VeryLowQual" 和 “LowQual”的标记(均可在参数中自行定义)。随后就是提取"PASS"类型的变异:

bcftools filter --threads <n> -i \'FILTER=="PASS"\' --no-version -O z -o SNP_filter.vcf.gz SNP_filted.vcf.gz

        到此，群体变异信息的合并以及初步筛选结束。下一步针对群体中MAF和缺失比例等信息的筛选参数则需要根据不同物种，测序深度以及分析目的来做调整:

vcftools --gzvcf SNP_filted.vcf.gz --remove-indels --maf 0.05 --max-missing 0.5 --recode --stdout |bgzip -f > SNP_maf0.05_miss0.5.vcf.gz
vcftools --gzvcf InDel_filted.vcf.gz --keep-only-indels --maf 0.05 --max-missing 0.5 --recode --stdout |bgzip -f > Indel_maf0.05_miss0.5.vcf.gz
##保留MAF>=0.05，以及最大缺失比例为0.5的变异位点
##此处--max-missing存在表述上的错误，如--max-missing 0.1表示至少在10%的材料中存在的变异，而非允许最多10%的缺失
##此处 --remove-indels 和 --keep-only-indels 的参数可以解决前面提到的位点重复问题

        最终，得到了该群体过滤后的所有变异信息，可以用作群体结构变异，建树等后续分析

PS：此流程为个人分析所得经验总结，欢迎大家交流指正。