目的：提取群体中多态性较好的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
