最近看了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
