GATK4这个工具我真的想疯狂吐槽一下，不断更新确实是好事，但是官网上关于同一个tool的各种manual以及Q&A给出的例子都不一致，简直抓狂，记录下自己的实践情况吧，非常辛酸。
2019.10.15：本记录中使用的是GATK4.1.2.0，不完全适用于最新版本的GATK4。
2021.12.03：教程更新为GATK4.2.3.0，增加hg38相关
2022.7.13：GATK官网持续改版，部分链接失效，更新链接

基本流程简单来说，分为三部分：数据预处理——变异检测——结果文件过滤。其中数据预处理部分是一样的，后两部分又分为寻找种系突变和体细胞突变。本文侧重寻找体细胞突变。
建议在代码化的时候设置好各项变量，比如参考文件路径，软件路径，输出文件路径等等，这样在批量出结果的时候就可以脉络清晰，方便查看及引用。

数据预处理

下机数据处理获得dup.bam文件

由于我拿到的数据是dup.bam（即经过比对、排序、去重、索引后），所以数据预处理这部分基本不需要做什么了，但是有以下几点注意事项。

  1. dup.bam文件必须有文件头，此处指的是要包含@RG信息，没有这个信息的话无法使用，可以用picard的AddOrReplaceReadGroups手动加。
$JAVA/java -Xmx2g -jar $PICARD/picard.jar AddOrReplaceReadGroups \
I=$sample.dup.bam O=$sample.ah.bam \
ID=$sample LB=$sample PL=illumina PU=$sample SM=$sample
  1. 看到很多bwa比对的，但我这里用的是bowtie2，如果要在bowtie2步骤把@RG加上，比对获得sam文件，排序去重建索引后，获得的dup.bam文件直接就可以进行下游分析了。
$BOWTIE/bowtie2 -q --phred33 --sensitive --end-to-end \
-x $ref/hg19/Bowtie2Index/human_hg19 \
--rg-id $sample \
--rg "PL:illumina" \
--rg "LB:$sample" \
--rg "PU:$sample" \
--rg "SM:$sample" \
-1 ${sample}_1.clean.fastq.gz -2 ${sample}_2.clean.fastq.gz -S ${sample}.sam
  1. 做变异检测应该都是文件头必须包含@RG，我测试Deepvariant时它也需要带@RG的bam文件。

碱基质量校正（BQSR）

BQSR:Base Quality Score Recalibration.
GATK对有@RG文件头的dup.bam文件进行碱基质量校正。

  1. 所需文件准备，可以从GATK官网Resource Bundle下载需要的参考文件。或者根据自己的文件生成，以hg19为例，需要hg19.fa，hg19.fa.fai，hg19.dict。没有.fa.fai可以用samtools的faidx生成，没有.dict可以用picard的CreateSequenceDictionary生成。
$SAMTOOLS/samtools faidx hg19.fa
$JAVA/java -Xmx2g -jar $PICARD/picard.jar CreateSequenceDictionary R=hg19.fa O=hg19.dict
  1. 碱基质量校正分两步。
  • 第一步，BaseRecalibrator的--known-sites参数需要的文件可以从GATK官网Resource Bundle下载，参考了github上的官方流程以及一些国内教程，几个文件的选择未有定论，官方只用了两个，其他教程有用三个的。具体参看了公众号“碱基旷工”的选择。几个文件的解释有相关介绍：官网解释。
$GATK BaseRecalibrator \
-R hg19.fa \
-I ${sample}.ah.bam -O ${sample}.recal.table \
--known-sites dbsnp_138.hg19.vcf \
--known-sites 1000G_phase1.indels.hg19.sites.vcf \
--known-sites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf
  • 第二步，获得校正后的文件，我这里的后缀为ah.recal.bam。记得samtools建立索引。
$GATK ApplyBQSR \
-bqsr ${sample}.recal.table \
-R hg19.fa \
-I ${sample}.ah.bam -O ${sample}.ah.recal.bam

可选步骤

可以生成一个pdf文件，里面有一些统计和图，其实我没看懂这些图。找到一个资源，有一些关于这些图的讲解，仅供参考GATK使用方法详解，系列文章里有不少关于原理方面的知识值得学习。

$GATK AnalyzeCovariates -bqsr ${sample}.recal.table -plots ${sample}.recal.pdf

变异检测

变异检测分为种系突变和体细胞突变，过程复杂到令人崩溃。

种系突变（HaplotypeCaller）

相对来说寻找种系突变还算简单点。以下方式是计算较慢，需要速度提升的话可以拆成对每个染色体做snp calling，最后合并。这里我考虑到chr*_random和chrUn_*可能有其意义，所以也没有拆成每个染色体寻找突变。加--dbsnp可以在结果vcf文件中注释一些研究过的snp。

$GATK HaplotypeCaller \
-R hg19.fa \
--dbsnp dbsnp_138.hg19.vcf \
-I ${sample}.ah.recal.bam \
-O ${sample}.hc.vcf \
-ERC GVCF \
-G StandardAnnotation \
-G AS_StandardAnnotation \
-bamout ${sample}.hc.bam

体细胞突变（Mutect2）

做肿瘤研究这一块更侧重做体细胞突变。Mutect2主要是根据对正常-肿瘤样本进行位点比较寻找突变，如果没有正常样本，那么软件能使用，但假阳性会很高。

直接对每个样本寻找体细胞突变

即不进行对照，也很简单快捷。这里要记得加--germline-resource。后面会提到该文件的获取途径。

$GATK Mutect2 \
-R hg19.fa \
--germline-resource af-only-gnomad.raw.sites.hg19.vcf.gz \
-I ${sample}.ah.recal.bam \
-O ${sample}.mu.vcf

正常-肿瘤样本对照寻找snp

非常复杂，需要构建PoN（Panel Of Normals），下载germline-resource。另外还有去除污染和链方向偏差矫正，这两块我没有做。

  1. Panel of Normals (PON)：官方建议使用健康的年轻的至少40例样本构建PoN，且测序实验技术上应尽可能与肿瘤样本相似，构建PoN的目的一是排除种系突变，二是排除测序技术误差。来自同一组织的更佳，不足40例样本也可以构建，聊胜于无，比没有强，我用了8个样本构建。官方也提供了Public GATK Panel of Normals，可以从其谷歌云存储里下载。
    PS：但是这里有一个问题，根据PoN的作用，我觉得构建PoN的健康样本的测序时间及测序条件应该与要寻找突变的肿瘤样本相近，这点不利于长期实践和标准流程构建，如果在生产过程中，要加入寻找突变的检测，长期下来，可能需设定PoN更新周期来更新PoN。
  • 第一步，Mutect2的Tumor-only mode对每个样本寻找体细胞突变。构建PoN的第一步不能加--germline-resource参数。此处由于后续工具有beta版的，导致这里必须加--max-mnp-distance 0参数，否则下一步报错或无法获得结果文件。（此处官网不同页面下的两个例子不一致，第一次我没加这个参数得到了血泪教训，重跑了好久）
    在GATK4.1.4.0的CreateSomaticPanelOfNormals部分有详细介绍，以下为引用原文：

Note that as of May, 2019 -max-mnp-distance must be set to zero to avoid a bug in GenomicsDBImport.

$GATK Mutect2 \
-R hg19.fa \
--max-mnp-distance 0 \
-I ${sample}.ah.recal.bam \
-O ${sample}.re.vcf.gz
  • 第二步，把第一步获得的结果合并。
merge_vcfs=""
while read sample
do
	  merge_vcfs=${merge_vcfs}" -V ${sample}.re.vcf.gz "
done<samplelist_pon

$GATK GenomicsDBImport \
-R hg19.fa \
-L hg19.interval_list \
--genomicsdb-workspace-path pon_db \
${merge_vcfs}

这里又有坑了GenomicsDBImport：在GATK4.1.2.0中该工具的-L是必须参数，我是针对全基因组的，实在不知道该传什么，就下载了hg19.bed，用picard的BedToIntervalList获得了hg19.interval_list。（picard功能好多，简直是救我狗命，好崩溃呀）
补充：对于外显子数据，还要加上一个参数--merge-input-intervals true，同一条染色体会整合到一起运行，并将结果保存到同一个文件夹中。另外如果在/tmp文件夹存储过小的服务器，还需指定临时文件目录，否则会因临时文件目录被写满而终止程序，使用--tmp-dir /data/group/tmp。

$JAVA/java -Xmx10g -jar $PICARD/picard.jar BedToIntervalList \
I=hg19.bed \
O=hg19.interval_list \
SD=human.dict
  • 第三步，将自己构建的pon_db和germline-resource相结合获得PoN。这个工具目前是beta版。
$GATK CreateSomaticPanelOfNormals \
-R hg19.fa \
--germline-resource af-only-gnomad.raw.sites.hg19.vcf.gz \
-V gendb://pon_db \
-O pon_db.vcf.gz

CreateSomaticPanelOfNormals是beta版的。这里的坑是用hg19的话，GATK官网Resource Bundle里面没有af-only-gnomad.raw.sites.hg19.vcf.gz文件，在讨论区看到回答guanyu--germline-resource参数，回答里有一个资源hg19资源，未经验证。还提供了另一种方法，官网下载b37的，转为hg19的，这个方法可通过picard的LiftoverVcf转换，需要一个b37tohg19.chain文件。（又是picard，好厉害）
备注：我下载了这个资源，用官方的b37也生成了一个自己的，两个文件大小不一致，统计其中记录的突变数据行数是一样的，在CreateSomaticPanelOfNormals这步获得的pon_db.vcf.gz结果也是一样，diff了一下，差别仅在文件头里里记录的执行语句及任务开始时间不同。所以我认为这个可以用。自己生成的是可以复现溯源。

$JAVA/java -Xmx10g -jar $PICARD/picard.jar LiftoverVcf \
I=af-only-gnomad.raw.sites.b37.vcf.gz \
O=af-only-gnomad.raw.sites.hg19.vcf.gz \
CHAIN=b37tohg19.chain \
REJECT=rejected_variants.vcf \
R=hg19.fa
  1. 对每对正常-肿瘤样本寻找体细胞突变。正常-肿瘤需要是来自同一个体的，我这里用的是TI、TU样本。这里-normal参数即正常样本bam文件头@RG里SM的值。
$GATK Mutect2 \
-R hg19.fa \
-I ${sample_tu}.ah.recal.bam \
-I ${sample_ti}.ah.recal.bam \
-normal ${sample_ti} \
--germline-resource af-only-gnomad.raw.sites.hg19.vcf.gz \
--panel-of-normals pon_db.vcf.gz \
-O ${sample_tu}_${sample_ti}.vcf

结果文件过滤

种系突变结果过滤（VQSR/hard-filtering）

过滤这里是种系突变过滤比较复杂。这里只讲snp的，indel的同理。

  1. 做VQSR过滤是建立了GMM模型进行过滤。要求样本量大，且有足够多的已知变异参考信息（比如人类）。分两步。
  • 第一步。
$GATK VariantRecalibrator \
-R hg19.fa \
-V ${sample}.hc.vcf \
--resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.hg19.sites.vcf \
--resource:omini,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.hg19.sites.vcf \
--resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg19.sites.vcf \
--resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_138.hg19.vcf \
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
-mode SNP \
--rscript-file ${sample}.plots.R \
--tranches-file ${sample}.hc.snp.tranches \
-O ${sample}.hc.snp.recal
  • 第二步。
$GATK ApplyVQSR \
-R hg19.fa \
-V ${sample}.hc.vcf \
--tranches-file ${sample}.hc.snp.tranches \
--recal-file ${sample}.hc.snp.recal \
-mode SNP \
-O ${sample}.hc.VQSR.vcf
  1. 但是不是所有数据都符合做VQSR的要求，不能做的就需要硬过滤（待验证），我没有用硬过滤，这里还要学习如何定义各个指标的阈值，也很复杂。
$GATK SelectVariants \
-select-type SNP \
-V ${sample}.hc.vcf \
-O ${sample}.hc.snp.vcf

$GATK VariantFiltration \
-V ${sample}.hc.snp.vcf \
--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name "Filter" \
-O ${sample}.hc.snp.filter.vcf

体细胞突变过滤

感觉这里可以用先苦后甜来形容了，一行命令搞定。

$GATK FilterMutectCalls \
-R hg19.fa \
-V ${sample}.mu.vcf \
-O ${sample}.mu.filter.vcf

参考文件下载

hg38的germline-resource下载地址，下载的文件叫af-only-gnomad.hg38.vcf.gz，对应索引也要下载。
hg38的interval_list下载地址，下载的文件叫whole_exome_illumina_coding_v1.Homo_sapiens_assembly38.targets.interval_list，论坛上关于该文件的说明，引用其中一句"In general you probably want to be using at the targets file, that's the region of interest"。
GATK官网Resource Bundle，持续改版，目前存放在谷歌云存储里，只包含部分文件，如有其他需求，可在GATK官方论坛里检索。

结语

获得的突变结果文件是vcf格式的，过滤前后的FILTER位有差异，过滤前的文件FILTER位是一个点号，过滤后的通过过滤的位点后为“PASS”，其他的还有各种标记，还需进一步研究。
这个GATK4找突变，真的很难很复杂。还在摸索中。
关于种系突变VQSR里硬过滤参数阈值设定、vcf文件格式说明之类的，相关知识和内容实在太多，在一篇文章中讲不完，再说我也不怎么会。

参考资料：
（1）深入了解snp-calling流程
（2）Mutect2-肿瘤不同转移时期的免疫微环境异质性研究
（3）初试GATK4-mutect2来call somatic mutation
（4）【原创】GATK使用方法详解（包含bwa使用）第三部分
（5）GATK官方网站
（6）GATK 的 Germline mutation 流程