GATK4.0 和之前的版本相比还是有较大的不同，更加趋于流程化。

软件安装

点击此处查看最新版本

wget https://github.com/broadinstitute/gatk/releases/download/4.1.5.0/gatk-4.1.5.0.zip
unzip gatk-4.1.5.0.zip

GATK 简单说明

## 帮助信息
gat --help

## 列出所有的工具
gatk --list

## 工具的说明，比如以VariantAnnotator 为例
gatk VariantAnnotator --help

GATK分析简要流程

所需数据

ref.fa
reads1.fq
reads2.fq

建立索引

bwa index ref.fa
samtools  faidx ref.fa
gatk CreateSequenceDictionary -R ref.fa -O ref.dict

##
-R Input reference fasta or fasta.gz  Required
-O  输出文件

比对

## bwa 比对
bwa mem -t 4 -R '@RG\tID:id1\tPL:illumina\tSM:test' ref.fa test_1.fq test_2.fq | samtools view -bS - >test.bam

##参数
-R 设置reads group，gatk必须要的信息，其中ID，PL和SM信息是必须要的

## 排序
samtools sort -@ 3 -o test.sorted.bam test.bam
rm test.bam


## 也可以用GATK的SortSam进行排序，可以对SAM，或者BAM直接排序
gatk SortSam -I test.bam -O test.sorted.bam -SO coordinate --CREATE_INDEX true
## c参数
-I： 输入bam或者sam
-O： 输出文件
-SO：排序方式：queryname 或者coordinate
--CREATE_INDEX： 是否建立索引

GATK 要求read group的格式

ID = Read group identifier
　每一个read group 独有的ID，每一对reads 均有一个独特的ID，可以自定义命名；
PL = Platform
　　测序平台；ILLUMINA, SOLID, LS454, HELICOS and PACBIO，不区分大小写；
SM = sample
　　reads属于的样品名；SM要设定正确，因为GATK产生的VCF文件也使用这个名字;
LB = DNA preparation library identifier
　　对一个read group的reads进行重复序列标记时，需要使用LB来区分reads来自那条lane;有时候，同一个库可能在不同的lane上完成测序;为了加以区分，
　　同一个或不同库只要是在不同的lane产生的reads都要单独给一个ID. 一般无特殊说明，成对儿read属于同一库，可自定义，比如：library1

若是忘记添加read group信息还以通过 AddOrReplaceReadGroups 添加

gatk AddOrReplaceReadGroups -I .bam -O .add.bam -LB library1 -PL illumina -PU pl1 -SM name

##参数
-I Input file (BAM or SAM or a GA4GH url)；
-O  Output file (BAM or SAM)；
-LB Read-Group library；
-PL  Read-Group platform (e.g. ILLUMINA, SOLID)；
-PU Read-Group platform unit (eg. run barcode)；
-SM Read-Group sample name

标记重复序列

gatk  MarkDuplicates -I test.sorted.bam -O test.sorted.markdup.bam -M test.sorted.markdup_metrics.txt
##参数
-I 排序后的bam或者sam文件
-M 输出重复矩阵
-O 输出文件

## 建立索引
samtools index test.sorted.markup.bam

检测变异

##两种方法

##（1）多样本一起call，此次只有一个样本，若有多个样本，则继续用 -I 参数添加即可
gatk --java-options -Xmx4G HaplotypeCaller -I test.sorted.markup.bam -O test.vcf -R ref.fa

##（2）单个样本call，然后在合并
## 生成中间文件gvcf
gatk --java-options -Xmx4G HaplotypeCaller -I test.sorted.markup.bam -O test.g.vcf -R ref.fa --emit-ref-confidence GVCF

##通过gvcf检测变异, -V 添加上步得到的gvcf
gatk GenotypeGVCFs -R ref.fa -V test.g.vcf -O test.vcf

##参数
-I BAM/SAM/CRAM file
-O  输出文件
-R 参考基因组
--java-options: 若设置java则需要添加
-Xmx4G：内存为4G，防止内存太大
-V  A VCF file containing variants
-L 第一种方法可单独对染色体分开进行call，而后用GatherVcfs可以合并，可加快速度

$\color{red}{==补充==}$
4.0以后GenotypeGVCFs只能接受single-sample GVCF from HaplotypeCaller or a multi-sample GVCF from CombineGVCFs orGenomicsDBImport，若有多个g.vcf 可以使用上述两种工具进行合并成一个单独的文件即可.

提取SNP，INDEL

## 提取SNP
gatk SelectVariants -V test.vcf -O test.snp.vcf --select-type-to-include SNP

## 提取INDEL
gatk SelectVariants -V test.vcf -O test.indel.vcf --select-type-to-include INDEL

##参数
-O 输出vcf文件
-V 输入vcf文件
--select-type-to-include 选择提取的变异类型{NO_VARIATION, SNP, MNP, INDEL,
                              SYMBOLIC, MIXED}

对vcf文件进行过滤

gatk VariantFiltration -O test.snp.fil.vcf.temp -V test.snp.vcf --filter-expression 'QUAL < 30.0 || QD < 2.0 || FS > 60.0 ||  SOR > 4.0' \
　　　　--filter-name lowQualFilter --cluster-window-size 10  --cluster-size 3 --missing-values-evaluate-as-failing

## 参数
-O 输出filt.vcf文件
-V 输入vcf文件
--filter-expression 过滤条件, VCF INFO 信息
--cluster-window-size 以10个碱基为一个窗口
--cluster-size 10个碱基为窗口，若存在3以上个则过滤
--filter-name 被过滤掉的SNP不会删除，而是给一个标签， 比如 Filter
--missing-values-evaluate-as-failing 当筛选标准比较多的时候，可能有一些位点没有筛选条件当中的一条或几条，例如下面的这个表达式；QUAL < 30.0 || QD < 2.0 || FS > 60.0 || MQ < 40.0 || HaplotypeScore > 13.0 并不一定所有位点都有这些信息，这种情况下GATK运行的时候会报很多WARNING信息，用这个参数可以把这些缺少某些FLAG的位点也给标记成没有通过筛选的。
## QualByDepth(QD): 变异位点可信度除以未过滤的非参考read数
## FisherStrand (FS): Fisher精确检验评估当前变异是strand bias的可能性，这个值在0-60间
# RMSMappingQuality (MQ): 所有样本中比对质量的平方根
# MappingQualityRankSumTest (MQRankSum): 根据REF和ALT的read的比对质量来评估可信度
# ReadPosRankSumTest (ReadPosRankSum) : 通过变异在read的位置来评估变异可信度，通常在read的两端的错误率比较高
# StrandOddsRatio (SOR) : 综合评估strand bias的可能性

筛选PASS的SNP，INDEL

## 根据FILTER那列信息进行筛选
grep PASS test.snp.fil.vcf.temp >  test.snp.fil.vcf

参考

Individual identifier (optional) - The previous column told us to expect to see genotypes here. The genotype is in the form 0|1, where 0 indicates the reference allele and 1 indicates the alternative allele, i.e it is heterozygous. The vertical pipe | indicates that the genotype is phased, and is used to indicate which chromosome the alleles are on. If this is a slash / rather than a vertical pipe, it means we don’t know which chromosome they are on.

GATK4.1 call SNP