用GATK进行二代测序数据 SNP Calling 流程：（四）变异过滤

GATK推荐的最好的过滤方式是用 VQSR功能，它通过机器学习算法来判断SNP的优劣，因此至少需要两个已存在的 SNP 数据集，一个是经过验证的高质量 SNP 数据集作为真集（如 HapMap），还需要一个质量不是特别高，允许存在小部分假阳性的数据集做训练集（如，1000G）。这些数据集在人类研究中很容易找到，但是在植物中比较困难，因此本流程用硬过滤（hard-filtering）的方法进行变异过滤。

提取SNP和INDEL

SNP 和 INDEL 的过滤参数有所不同，因此分开过滤。

#vcf索引
nohup gatk SelectVariants -V raw.vcf.gz -select-type SNP -O snp..vcf.gz &
nohup gatk SelectVariants -V raw.vcf.gz -select-type INDEL -O indel.vcf.gz &

SNP 过滤

查看过滤参数的分布情况

可以先查看一下想用的过滤参数的分布情况

#使用gatk的variantstotable功能提取过滤参数信息
nohup gatk VariantsToTable -V snp.vcf.gz -F CHROM -F POS -F QD -F QUAL -F SOR -F FS -F MQ -F MQRankSum -F ReadPosRankSum -O snps.recode.table &
nohup gatk VariantsToTable -V indel.vcf.gz -F CHROM -F POS -F QD -F QUAL -F SOR -F FS -F MQ -F MQRankSum -F ReadPosRankSum -O indel.recode.table &

使用gatk VariantFiltration 过滤 SNP

版本：The Genome Analysis Toolkit (GATK) v4.1.8.0、Picard Version: 2.22.8

gatk=gatk
# gz格式需要用tabix索引
vcf=snp.vcf.gz
out_prefix=snp.hardfiltereannotated

$gatk VariantFiltration \
    -V $vcf \
    -filter "QD < 2.0" --filter-name "QD2" \
    -filter "QUAL < 30.0" --filter-name "QUAL30" \
    -filter "FS > 60.0" --filter-name "FS60" \
    -filter "MQ < 40.0" --filter-name "MQ40" \
    -filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
    -filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \
    -O ${out_prefix}.vcf

1. 过滤参数使用GATK建议参数。SOR参数在很多文章中没有用到，且与FS参数功能重合，此处也不用了
2. MQRankSum 和ReadPosRankSum的计算需要至少有一个杂合样本，所以有些位点可能会出现没有这两个INFO的情况，如果用复合过滤式会将这些位点过滤掉，因此不用复合式
3. gatk输出文件格式会自动识别后缀，为方便后期awk处理，不再输出.gz压缩格式

VariantFiltration 功能仅会在 VCF 文件中对每个变异位点添加过滤注释，因此还需要跟过滤注释提取出通过过滤标准的变异。

#提取通过过滤的变异（PASS），对于没有注释 MQRankSum 和ReadPosRankSum的位点也保留下来
nohup awk '/^#/||$7=="PASS"||$7=="MQRankSum-12.5;ReadPosRankSum-8"||$7=="MQRankSum-12.5"||$7=="ReadPosRankSum-8"' snp.hardfiltereannotated.vcf > snp.hardfiltered.vcf &

# 剔除非2个等位基因的位点。对于二倍体物种可以仅保留双等位基因位点。
nohup awk '$5 !~ /([[:alpha:]]|*)+,([[:alpha:]]|*)/{print}' snp.hardfiltered.vcf > snp.hardfiltered.biallelic.vcf &
# 用 VCFtools 分别过滤一份完整度80%和50%的vcf
nohup vcftools --vcf snp_hardfiltered_biallelic.vcf --maf 0.05 --max-missing 0.8 --recode --recode-INFO-all --out snp_hardfiltered_biallelic_maf5_ms80 &

nohup vcftools --vcf snp_hardfiltered_biallelic.vcf --maf 0.05 --max-missing 0.5 --recode --recode-INFO-all --out snp_hardfiltered_biallelic_maf5_ms50 &

查看个体杂合、HWE、个体和位点深度

#nohup vcftools --vcf snp_hardfiltered_biallelic_maf5_ms50.vcf.recode.vcf --het --out snp_hardfiltered_biallelic_maf5_ms50 &
#nohup vcftools --vcf snp_hardfiltered_biallelic_maf5_ms50.vcf.recode.vcf --depth --out snp_hardfiltered_biallelic_maf5_ms50 &
#nohup vcftools --vcf snp_hardfiltered_biallelic_maf5_ms50.vcf.recode.vcf --site-depth --out snp_hardfiltered_biallelic_maf5_ms50 &
#nohup vcftools --vcf snp_hardfiltered_355_biallelic_maf5_ms50.vcf.recode.vcf --hardy --out snp_hardfiltered_biallelic_maf5_ms50 &

INDEL过滤

筛选出<10bp，双等位基因的 INDEL

nohup gatk SelectVariants -V indel.utx.vcf.gz --O indel.biallelic.vcf --max-indel-size 9 --restrict-alleles-to BIALLELIC &

使用gatk VariantFiltration 过滤 indel

## 过滤注释
nohup gatk VariantFiltration -V indel.biallelic.vcf -O indel.biallelic.10bp.hardfilterannotated.vcf -filter "QD < 2.0" --filter-name "QD2" -filter "QUAL < 30.0" --filter-name "QUAL30" -filter "FS > 200.0" --filter-name "FS200" -filter "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20"

# 提取通过过滤的变异
nohup awk '/^#/||$7=="PASS"||$7=="ReadPosRankSum-20"' indel.biallelic.10bp.hardfilterannotated.vcf > indel.biallelic.10bp.hardfiltered.vcf &
# 过滤一份完整度 50% 的SNP
nohup vcftools --vcf indel.biallelic.10bp.hardfiltered.vcf --maf 0.05 --max-missing 0.5 --recode --recode-INFO-all --out indel.biallelic.10bp.hardfiltered.maf5ms50 &