GATK流程生成多态性的 SNP/InDel

目的：提取群体中多态性较好的SNP和InDel。将GVCF文件进行合并，并提取SNP和InDel，然后使用vcftools获得在群体中多态性较好的SNP和InDel。

第一步：将GVCF文件添加到GATK的GenomicsDBImport中:

#!/usr/bin/bash

gatk=/120t/caix/src/gatk-4.1.3.0/gatk
# 必需的文件，具体格式参考GATK
interval_list=CXB.bed
# Java内存限制
Xmx=200g
Xms=200g
# genomicsdb，必须是未创建的或者空的文件夹
my_database=CXBGenomicsdb
batch_size=10    ## 每批次读入多少样本
map_file=gvcfs.lst  ##- gvcf文件列表   格式：名称\t文件路径
reader_threads=10 # 每个batch的并行线程
num=10    ##-  同时读入的最大intervals数量
#缓存文件
tmp=tmpdir

mkdir ${tmp}

#第一次建库用的参数：
# --genomicsdb-workspace-path $my_database \
#向已存在的库添加新数据的参数：
#--genomicsdb-update-workspace-path $my_database \

$gatk --java-options "-Xmx$Xmx -Xms$Xms" \
GenomicsDBImport \
--genomicsdb-workspace-path $my_database \
--intervals $interval_list \
--batch-size  $batch_size \
--sample-name-map $map_file \
--reader-threads $reader_threads \
--max-num-intervals-to-import-in-parallel $num \
--tmp-dir ${tmp}

第二步：提取群体的vcf文件：

/120t/caix/src/gatk-4.1.3.0/gatk GenotypeGVCFs -R CXB.fasta -V gendb://CXBGenomicsdb -O CXB_output.vcf.gz

第三步：提取SNP和InDel位点：

/120t/caix/src/gatk-4.2.4.1/gatk  SelectVariants   -V   A04_output.vcf.gz  -O  A04.snp.vcf  --select-type-to-include SNP
/120t/caix/src/gatk-4.2.4.1/gatk   SelectVariants    -V   A04_output.vcf.gz   -O   A04.indel.vcf   --select-type-to-include INDEL

第四步：hard-filter 标记低质量的SNP和InDel，并使用vcftools进行过滤

#!/bin/bash
snp=$1
indel=$2
prefix=$3

gatk="/120t/caix/src/gatk-4.2.4.1/gatk"
vcftools="/10t/caix/miniconda2/bin/vcftools"

$gatk   VariantFiltration \
        -V  ${snp}  \
        -filter "QD < 2.0" --filter-name "QD2" \
        -filter "QUAL < 30.0" --filter-name "QUAL30" \
        -filter "SOR > 3.0" --filter-name "SOR3" \
        -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 ${prefix}.snps_ft.vcf.gz


$gatk   VariantFiltration \
        -V  ${indel}   \
        -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" \
        -O ${prefix}.indels_ft.vcf.gz

$vcftools  \
        --minDP 3 \
        --maxDP 500 \
        --minGQ  10 \
        --minQ 30 \
        --min-meanDP 3 \
        --keep-filtered PASS \
        --out   ${prefix}.snps_PASS_miss0.5.maf0.05.vcf  \
        --gzvcf ${prefix}.snps_ft.vcf.gz  \
        --recode --recode-INFO-all \
        --max-missing 0.5 \
        --maf 0.05

$vcftools  \
        --minDP 3 \
        --maxDP 500 \
        --minGQ  10 \
        --minQ 30 \
        --min-meanDP 3 \
        --keep-filtered PASS \
        --out   ${prefix}.indels_PASS_miss0.5.maf0.05.vcf  \
        --gzvcf ${prefix}.indels_ft.vcf.gz   \
        --recode --recode-INFO-all \
        --max-missing 0.5 \
        --maf 0.05