GATK4 —— 获取短变异 （call SNP+indel）

GATK是一款用于基因组数据分析的软件，其强大的处理引擎和高性能计算功能使其能够承担任何规模的项目。

GATK的功能非常强大，这里不详细介绍，大家可以根据自己的要求，从首页进入对应的模块，说明书还是很详细的。

官网：
https://gatk.broadinstitute.org/hc/en-us

GATK 4.2.1.0 命令说明：
https://gatk.broadinstitute.org/hc/en-us/sections/4404575821339-4-2-1-0?page=10#articles

1. 下载安装

1.1下载

下载链接：
https://github.com/broadinstitute/gatk/releases

1.2 安装

解压：

$ unzip gatk-4.2.4.1.zip

查看帮助：

$ cd /your/path/gatk-4.2.4.1
$ ./gatk --help

查看工具列表：

$ ./gatk --list

查看指定工具帮助：

$ ./gatk ToolName --help

2.获取短变异（SNP+indel）

说明书：
https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels-
主要流程如下：

2.1 准备数据

• 这里使用在之前已经比对好并排序去重后的bam文件（.bam），详情参见以下链接：
  测序质量分析 —— FastQC - 简书 (jianshu.com)
  gff格式与gtf格式转换——NBISweden / AGAT - 简书 (jianshu.com)
  序列比对 —— Hisat2 - 简书 (jianshu.com)
  sam文件转换为bam文件——SAMtools - 简书 (jianshu.com)
  SAMtools——bam文件排序 - 简书 (jianshu.com)
  SAMtools——bam文件去重 - 简书 (jianshu.com)

• bam文件索引（.bai）

$ samtools index LPF1_R1_MP.sort.dup.bam

• 参考基因组（.fa）
• 参考基因组索引（.fai）

$ samtools faidx Mparg_v2.0.fa
$ less Mparg_v2.0.fa.fai
Mpar_chr1       41449654        11      100     101
Mpar_chr2       50537509        41864173        100     101
Mpar_chr3       46570116        92907069        100     101
Mpar_chr4       46692298        139942898       100     101
Mpar_chr5       50691827        187102130       100     101
Mpar_chr6       65498417        238300887       100     101
Mpar_chr7       34026340        304454300       100     101
Mpar_chr8       48943667        338820915       100     101

五列分别对应，序列名、序列长度、第一个碱基的偏移量、行的碱基数、行宽。

• 参考基因组索引（.dict）

$ gatk CreateSequenceDictionary -R Mparg_v2.0.fa -O Mparg_v2.0.dict

2.2 HaplotypeCaller

HaplotypeCaller是GATK4中的核心工具，可以利用单倍型区域的局部从头组装同时调用种系snv和小indel。详细原理参见说明书。
https://gatk.broadinstitute.org/hc/en-us/articles/360050814612-HaplotypeCaller#--sample-name

$ ./gatk HaplotypeCaller --help

Required Arguments:

--input,-I <GATKPath>         BAM/SAM/CRAM file containing reads  This argument must be specified at least once.
                              Required.

--output,-O <GATKPath>        File to which variants should be written  Required.

--reference,-R <GATKPath>     Reference sequence file  Required.

示例：

$ gatk --java-options "-Xmx100g -XX:ParallelGCThreads=4" HaplotypeCaller \
    -I ~/LPF1_R1_MP.sort.dup.bam \
    -R ~/Mparg_v2.0.fa \
    -ERC GVCF \
    -O ~/vcf/LPF1_R1_MP.g.vcf.gz

报错：A USER ERROR has occurred: Argument emit-ref-confidence has a bad value: Can only be used in single sample mode currently. Use the --sample-name argument to run on a single sample out of a multi-sample BAM file.

解决办法：
picard——修改BAM文件的Read Group - 简书 (jianshu.com)

更改命令：

$ samtools index LPF1_R1_MP.rg.sort.bam
$ gatk --java-options "-Xmx100g -XX:ParallelGCThreads=4" HaplotypeCaller \
    -I ~/LPF1_R1_MP.rg.sort.bam \
    -R ~/Mparg_v2.0.fa \
    -ERC GVCF \
    -O ~/LPF1_R1_MP.g.vcf.gz

gatk好像最多只支持4个线程，如果我我理解有错误的话，还请大佬指正。

这里注意HaplotypeCaller只能处理单样本文件，当有多样本时，官方建议使用HaplotypeCaller对单bam文件分别进行变异检测，生成GVCF文件，GVCF会记录每一个位点到情况，包括有无突变，VCF只记录突变位点情况，之后在下一步对GVCF文件进行合并。

VCF文件和GVCF文件的格式的详细说明参见下链接：
VCF和GVCF格式说明 - 青萍，你好 - 博客园 (cnblogs.com)

3.GVCF文件合并

官方说明书：https://gatk.broadinstitute.org/hc/en-us/articles/4404604801947-CombineGVCFs

gatk --java-options "-Xmx100g -XX:ParallelGCThreads=4" CombineGVCFs \
    -R ~/ref/Mparg_v2.0.fa \
    -V LPF1_R1_MP.g.vcf.gz \
    -V LPF1_R2_MP.g.vcf.gz \
    -V LPF1_R3_MP.g.vcf.gz \
    -O LPF1_MP.g.vcf.gz

参考资料：
https://www.melbournebioinformatics.org.au/tutorials/tutorials/variant_calling_gatk1/variant_calling_gatk1/

https://gencore.bio.nyu.edu/variant-calling-pipeline-gatk4/