最近看了GATK的网站，以前经常用GATK call somatic mutation， 但是用的版本太老了，另外服务器最近换了盘，很多代码路径都变了，网上关于GATK4的教程比较少，不想做伸手党，于是尝试自己写代码试一试。

1 Getting start with GATK

GATK是Genome AnalysisToolkit的缩写，它搜集了一系列分析高通量测序数据的工具，主要目的是在于发现变异，这个工具可以单独的用，但是也可以和其他的workflow合在一起。

contents：

1. 安装GATK

GATK 工具具有相当简单的软件需求： 一个Unix-style OS 以及java 1.8 版本以上，其他的工具需要额外的R和python的依赖。 可以在这个网址下载： https://github.com/broadinstitute/gatk/releases，
最新版本的GATK可以按照如下命令下载：

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

如果需要其他版本，可以自行在在这个网址下载。 一旦你下载了解压这个文件夹，里面一般会有这四个文件：

gatk
gatk-package-[version]-local.jar
gatk-package-[version]-spark.jar
README.md

这个时候你可能会问，为什么会有两个jar， 就像名字里一样，gatk-package-[version]-spark.jar 是在Spark cluster上专门跑Spark tools 设计的。
而 gatk-package-[version]-local.jar是为了跑其他的设计的。那么这是否意味着你必须每次跑程序时指定你要跑哪一个，大可不必，在文件夹中有gatk这个参数，这是一个可执行的脚本将会根据你的命令行选择适当的jar，如果你想的话，当然可以用一个specific的jar，但是用gatk会简单一些。 它会根据你的情况帮你设置一些你需要的参数。

2. 是否安装成功

为了测试你是否成功的安装了GATK，在你的终端跑下面的命令行，

./gatk -- help 

3.跑GATK和picard 的命令行

 gatk [--java-options "jvm args like -Xmx4G go here"] ToolName [GATK args go here]

比如说：

gatk --java-options "-Xmx8G" HaplotypeCaller -R reference.fasta -I input.bam -O output.vcf

4.GATK Best Practices

GATK的输入文件需要做一些预处理。

（1）数据预处理

在所有call 点突变的步骤前必须的一个阶段，数据的预处理。 他包括处理原始的测序数据（FASTQ或者uBAM格式）去产生可以分析的BAM文件。 这些步骤包含了比对到reference 基因组上，以及一些数据清理的操作，矫正技术偏差，将数据更加合适分析。

输入格式

希望输入的read data的格式是unmapped BAM (uBAM) 格式。 转换工具可以将数据从FASTQ转换到uBAM。

（2）主要步骤

首先，我们把这个测序reads reads 比对到ref 基因组上去产生一个SAM/BAM格式的文件，这个SAM/BAM格式的数据根据坐标顺序排列。下一步，我们标记了重复序列为了减轻由于数据产出步骤而带来的误差（PCR扩增）。 最后我们矫正了碱基的质量得分，因为call 突变的算法严重依赖每条测序read中每个碱基的质量得分。

1 比对

在这里需要用到的工具是BWA

$RG_string="@RG ID:H0164.2  PL:illumina PU:H0164ALXX140820.2    LB:Solexa-272222    PI:0    DT:2014-08-20T00:00:00-0400 SM:NA12878  CN:BI"
$RG_LB="sample_name"
bwa mem -t 16 -M -R "$RG_string"  fq.1.gz fq.2.gz -o $RG_LB.sam"

2 mark duplicates (标记重复)

gatk --java-options "-Xmx20G" MarkDuplicatesSpark
-I $sample.bam -O ${sample}_marked.bam \
-M $sample.metrics \
1>log.mark 2>&1

3 FixMateInformation

gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" FixMateInformation \
    -I ${sample}_marked.bam \
    -O ${sample}_marked_fixed.bam \
    -SO coordinate \
    1>${sample}_log.fix 2>&1 

4 GATK碱基质量重矫正（BQSR）：

下载所需要的参考基因组

wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/*

此时会下载如下文件

1000G_phase1.snps.high_confidence.hg38.vcf
dbsnp_146.hg38.vcf
Mills_and_1000G_gold_standard.indels.hg38.vcf

然后对碱基进行矫正

gatk --java-options "-Xmx20G" BaseRecalibrator \
   -I ${sample}_marked_fixed.bam \
   -R reference.fasta \
   --known-sites 1000G_phase1.snps.high_confidence.hg38.vcf\
   --known-sites dbsnp_146.hg38.vcf\
   --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf\
   -O recal_data.table1

gatk  --java-options "-Xmx20G" ApplyBQSR \
   -R reference.fasta \
   -I ${sample}_marked_fixed.bam \
   --bqsr-recal-file  recal_data.table1 \
   -O ${sample}_marked_fixed_BQSR.bam

gatk --java-options "-Xmx20G" BaseRecalibrator \
   -I ${sample}_marked_fixed_BQSR.bam \
   -R reference.fasta \
   --known-sites 1000G_phase1.snps.high_confidence.hg38.vcf\
   --known-sites dbsnp_146.hg38.vcf\
   --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf\
   -O recal_data.table2

gatk --java-options "-Xmx20G" AnalyzeCovariates \
     -before recal_data.table1 \
     -after recal_data.table2 \
     -plots AnalyzeCovariates.pdf

5. 用mutect2 call点突变

目前mutect2 v4.1.0.0是支持多个样本的综合分析的。 里面有三种模式（1）tumor-normal mode。 （2） tumor-only mode。 （3） mitochondrial mode。

a 在v4.1，不再需要用-tumor 这个参数去特异的设置tumor 样本。 仅仅是需要用-normal 这个参数去设置正常样本，如果你有正常样本的话。

b 在4.0.4.0 开始，GATK提倡使用默认参数--af-of-alleles-not-in-resource， 工具中不同的模型的参数设置不同。 tumor-only calling 设置默认的是5e-8. tumor-normal calling 设置的是1e-6。 以及mitochondrial mode设置的是4e-3。 在以前的版本里，默认这是的是0.001， 这个是人类的平均异质性。 关于其他的物种，把--af-of-alleles-not-in-resource 改成1。

在这里，假设你已经有了normal.bam和tumor.bam, 参考文件是ref.fasta, 目标区域文件是intervals.interval.list .
(0) 构建normal panel
如果你有多个normal.bam做对照，在开始之前，可以利用这些normal bam构建normal panel作为第二步的原始候选变异检测输入集。

gatk Mutect2 -R ref.fasta  -I normal1.bam -max-mnp-distance 0 -O normal1.vcf.gz

gatk GenomicsDBImport \
-R reference.fasta \
-L intervals.interval_list \
--genomicsdb-workspace-path pon_db \
-V normal1.vcf.gz

gatk CreateSomaticPanelOfNormals \
-R reference.fasta \
-V gendb://pon_db \
-O pon.vcf.gz

（I）Tumor with matched normal
如果有一个匹配的normal， mutect2是被设计仅仅是用来call somatic mutation的。这个工具可以自动的跳过那些已知的存在于germline中的突变。 这样可以避免一些计算资源的浪费。 如果一个突变处于germline的边缘。mutect2会把这个变量放到callset数据集中，被FilterMutectCalls用来做进一步的过滤，或者review。

 gatk --java-options "-Xmx20G" Mutect2 \
     -R reference.fa \
     -I tumor.bam \
     -I normal.bam \
     -normal normal_sample_name \
     --germline-resource af-only-gnomad.vcf.gz \
     --panel-of-normals pon.vcf.gz \
     --f1r2-tar-gz f1r2.tar.gz \
     -O somatic.vcf.gz

这个命令产生的输出文件还需要用FiterMutectCalls过滤。
同时在mutect2的v4.1还支持来自于同一个人的多个肿瘤和正常样本。 但是必须要设置这个-I和-normal参数。

gatk --java-options "-Xmx20G" Mutect2 \
     -R reference.fa \
     -I tumor1.bam \
     -I tumor2.bam \
     -I normal1.bam \
     -I normal2.bam \
     -normal normal1_sample_name \
     -normal normal2_sample_name \
     --germline-resource af-only-gnomad.vcf.gz \
     --panel-of-normals pon.vcf.gz \
     -O somatic.vcf.gz

估计配对样本的交叉污染估计

gatk Pileup \
-R reference.fasta \
-L intervals.interval_list \
-I normal.bam \
-O normal-pileups.table

gatk Pileup \   
-R reference.fasta \   
-L intervals.interval_list \   
-I tumor.bam \   
-O tumor-pileups.table

gatk CalculateContamination \   
-I tumor-pileups.table \   
-matched normal-pileups.table \   
-O contamination.table

测序偏好矫正以及过滤

gatk LearnReadOrientationModel \
-I f1r2.tar.gz \
-O read-orientation-model.tar.gz

gatk FilterMutectCalls \
-R reference.fasta \
-V somatic.vcf.gz \
--contamination-table contamination.table \
--ob-priors read-orientation-model.tar.gz \
-O filtered.vcf.gz\

最后得到的filtered.vcf.gz就是过滤好的结果啦，vep注释起来，发现GATK4.1也提供了一个注释的工具Funcotator，有兴趣也可以尝试一下。
（ii）Tumor-only mode
这个模型是针对一个样本。

 gatk --java-options "-Xmx20G"  Mutect2 \
   -R reference.fa \
   -I sample.bam \
   -O single_sample.vcf.gz

(iii) Mitochondrial mode

 gatk Mutect2 \
  -R reference.fa \
  -L chrM \
  --mitochondria-mode \
  -I mitochondria.bam \
  -O mitochondria.vcf.gz