拿到 VCF 文件以后，怎么过滤才能得到比较可靠的变异位点呢？本文主要对 BWA + GATK 流程得到的 VCF 文件进行探讨，且仅涉及 SNPs 位点的过滤、不涉及 INDELs。
VCF 文件的格式
VCF 文件的格式请看这里。
VCF 文件的过滤
VCF 文件的过滤可以考虑以下几个方面：
1) 过滤 INDEL 附近的 SNPs；
2) 提取 SNP 位点；
3) GATK hard filtering；
4) 过滤 DP、GQ；
5) 过滤 maf、missing；
6) 仅保留具有2个等位基因的位点；
7) 过滤掉高杂和度的位点。
这些过滤信息仅供大家参考，每个工作的侧重点不同，所做的过滤也不一样。接下来看一下每一步过滤如何操作：
1 过滤 INDEL 附近的 SNPs
bcftools filter --SnpGap 5 -O v --output filtered_SNPs_near_INDEL.vcf raw.vcf
即使运行了 Indel realignment，并且使用 haplotypecaller，GATK 也不能完美地处理 INDELs，换言之，INDELs 附近的 SNPs 可能是不准确的，可以考虑在分析之前过滤掉。
上述命令中 raw.vcf 为初始的未过滤的 VCF 文件；--SnpGap 代表过滤掉 INDEL 附近 5 bp 以内的 SNPs；-O v 表示输出格式为 vcf；如果文件较大可以考虑设为 -O z，输出压缩的 VCF 文件； --output 用来指明输出文件的名称。
2 提取 SNP 位点
GATK=/your/path/to/GenomeAnalysisTK.jar
ref=/your/path/to/ref.fa
java -jar "$GATK" -T SelectVariants -R "$ref" -V filtered_SNPs_near_INDEL.vcf -selectType SNP -o SNPs_only.vcf
vcftools --vcf filtered_SNPs_near_INDEL.vcf --remove-indels --recode --recode-INFO-all --out SNPs_only
上面用两种方式实现提取 SNP 位点的操作；
如果使用 vcftools，记得加上 --recode-INFO-all 参数，这样才能保留 INFO 的信息，后面还需要根据这列信息进行过滤。如果不加，整列信息将变为 '.'！
3 GATK hard filtering
java -jar "$GATK" -T VariantFiltration -R "$ref" -V SNPs_only.vcf --filterExpression "QD < 2.0" --filterName "QD_snp_filter" --filterExpression "FS > 60.0" --filterName "FS_snp_filter" --filterExpression "SOR > 4.0" --filterName "SOR_snp_filter" --filterExpression "MQ < 40.0" --filterName "MQ_snp_filter" --filterExpression "MQRankSum < -12.5" --filterName "MQRankSum_snp_filter" --filterExpression "ReadPosRankSum < -8.0" --filterName "ReadPosRankSum_snp_filter" -o hard_filtering.snps.vcf
grep -v 'snp_filter' hard_filtering.snps.vcf > hard_filtered.snps.vcf
对于 GATK hard filtering 的参数，请参考官网的介绍；
这步运行后，GATK 并没有直接删掉不符合阈值的位点，仅仅是在 FILTER 那列标注 'PASS' 或标注每个过滤参数自定义的 'filterName'；
在最初尝试过滤参数组合的时候，给每个过滤参数的命名可参考：特有 (FS) + 共有 (_snp_filter) 的方式，以便查看每个过滤参数对数据的影响。
4 过滤 DP、GQ
vcftools --vcf hard_filtered.snps.vcf --minGQ 20 --minDP 5 --recode --out all.filter.DP.GQ
GQ 代表 genotype quality，和碱基质量值的表示一样 ，20 表示错误率为 1%；
DP 是过滤后的覆盖度，一般 2 倍体不能低于 5， 单倍体不能低于 3， 如果测序深度很高，可以考虑用 10；
需要注意的是：这步过滤仅仅是把不符合要求的个体基因型改为 './.'，并没有删掉任何的 SNPs！这也是需要在过滤 maf、missing 之前过滤 DP、GQ 的原因。
5 过滤 maf、missing
vcftools --vcf all.filter.DP.GQ.recode.vcf --max-missing 0.8 --maf 0.05 --max-maf 0.95 --recode --out all.rmMafMissing
注意 vcftools 的 --max-missing 1.0 代表没有 missing，设为 0.8 代表 missing 率为 20%；
6 仅保留具有2个等位基因的位点
vcftools --vcf all.rmMafMissing.recode.vcf --min-alleles 2 --max-alleles 2 --recode --out biallilic
含有多个等位基因的位点也可能是同源基因造成的，过滤后可提高 SNPs 准确性；
很多分析不支持含有 2 个以上的等位基因的位点；
7 过滤掉高杂和度的位点
如果一个位点的杂和度过高，很可能是旁系同源基因或其它因素的影响，需要进行过滤以保证 SNPs 的准确性。可以考虑使用 STACKS 软件中的 populations 命令来过滤，或者自己写一个简单的脚本来过滤。
相关阅读：
BWA + GATK 最佳实践及经验介绍
重测序分析流程及分析经验介绍
