简介
基因组测序，最重要的就是检测变异位点，对于家系数据、遗传病研究，更多的是关心 Germline mutation 生殖突变。当然，部分肿瘤研究也会关注 Germline mutation。GATK 对这类变异的检测有一整套流程，主要用到的工具是：HaplotypeCaller 、GenomicsDBImport、GenotypeGVCFs、VariantRecalibrator、 ApplyVQSR 等工具
流程图
Germline mutation 分析，对样本没有太多的要求，肿瘤非配对样本也可以分析。不过方法上有两种，单个样本和多个样本（队列）略有不同。流程图是： 对于多样本或队列样本：
对于单个样本：
在这里，仅介绍多样本或者队列样本的 GATK Germline mutation 流程，基于 BQSR 得到的 BAM 进行分析，主要有以下几步：
HaplotypeCaller
第一步先对每个样本 call 突变，用到了 HaplotypeCaller ，而且是在 GVCF 模式下，代码是：
代码语言：javascript复制${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" HaplotypeCaller
-R ${GENOME}
-ERC GVCF
-I 5.gatk/${id}_bqsr.bam
-O 6.gvcf/${id}.g.vcf.gz
--intervals ${bed}
1>./6.gvcf/${id}_HC.log 2>&1GenomicsDBImport
第二步是对第一步结果的整合，如果有50个样本，就有 50 个 *g.vcf.gz 文件，最后得到一个类似数据库的文件夹。代码是：
${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" GenomicsDBImport
-R ${GENOME}
$(ls 6.gvcf/*g.vcf.gz | awk '{print "-V "$0" "}')
-L ${bed}
--merge-input-intervals TRUE
--genomicsdb-workspace-path ./6.gvcf/gvcfs_db
1>./6.gvcf/GenomicsDBImport.log 2>&1
需要注意的是，对于外显子数据，除了指定 bed 文件，还要加上一个参数 --merge-input-intervals TRUE 。如果不加，对于每一个 bed 文件上的坐标（即bed文件的每一行），程序就会循环一次，并在 ./6.gvcf/gvcfs_db 文件夹中生成一个子文件夹，如果 bed 文件有 20W 行，就会有 20W 个文件夹，每个文件夹大小~100M，这个数据量是非常恐怖的。不仅如此，运行时间也大大增加。而加上参数 --merge-input-intervals TRUE 后，程序会对 bed 文件中的坐标进行整合，同一条染色体会整合到一起运行，并将结果保存到同一个文件夹中。
GenotypeGVCFs
第三步成为联合基因分型，按照官方文档所述：在这一步，收集所有每个样本的 GVCF 结果（已经在上一步整合为数据库）并将它们一起传递给联合基因分型工具 GenotypeGVCF。这会产生一组联合调用的 SNP 和 indel ，准备进行过滤。这种队列范围的分析即使在困难的位点也能灵敏地检测变异，并产生一个平方的基因型矩阵，该矩阵提供有关所考虑的所有样本中所有感兴趣位点的信息，这对于许多下游分析很重要。 此步骤运行速度非常快，可以在将样本添加到队列时随时重新运行，从而解决所谓的 N 1 问题。 代码是：
代码语言：javascript复制${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" GenotypeGVCFs
-R ${GENOME}
-V gendb://6.gvcf/gvcfs_db
-O ./6.gvcf/germline_merge.vcf.gz
1>./6.gvcf/GenotypeGVCFs.log 2>&1VariantRecalibrator
在数据预处理有一步是 BQSR 碱基质量重矫正，而这一步则是变异质量重矫正，两者是两回事，用到的工具也不同，需要区分开。 这一步实际上是基于机器学习的方法，对原始的 vcf 文件进行变异质量重矫正并且进行过滤。不过存在一个的缺点：该算法需要高质量的已知变体集作为训练和真实资源，而对于许多生物来说，这些资源尚不可用。它还需要相当多的数据来了解好与坏变体的概况，因此在仅涉及一个或几个样本的小数据集、靶向测序数据、RNAseq 上使用可能很困难甚至不可能使用，以及非模式生物。对于上述提到的情况，需要改用硬过滤的方法，可以参考：Hard-filtering germline short variants 代码是：
代码语言：javascript复制 ## 首先是对 SNP 位点运行 VariantRecalibrator
${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" VariantRecalibrator
-R ${GENOME}
-V 6.gvcf/germline_merge.vcf.gz
--max-gaussians 4
--resource:hapmap,known=false,training=true,truth=true,prior=15.0 ${hapmap}
--resource:omni,known=false,training=true,truth=false,prior=12.0 ${omni}
--resource:1000G,known=false,training=true,truth=false,prior=10.0 ${phase1}
--resource:dbsnp,known=true,training=false,truth=false,prior=2.0 ${dbsnp}
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR
-mode SNP
-O 6.gvcf/snp.recal
--tranches-file 6.gvcf/snp.tranches
--rscript-file 6.gvcf/snp.plots.R
1>./6.gvcf/SNP_VariantRecalibrator.log 2>&1
## 接着是对 INDEL 位点运行 VariantRecalibrator
${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" VariantRecalibrator
-R ${GENOME}
-V 6.gvcf/germline_merge.vcf.gz
--resource:mills,known=false,training=true,truth=true,prior=12.0 ${indel}
--resource:dbsnp,known=true,training=false,truth=false,prior=2.0 ${dbsnp}
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR
--resource:axiomPoly,known=false,training=true,truth=false,prior=10.0 ${Axiom_Exome}
-mode INDEL
-O 6.gvcf/indel.recal
--tranches-file 6.gvcf/indel.tranches
--rscript-file 6.gvcf/indel.plots.R
1>./6.gvcf/INDEL_VariantRecalibrator.log 2>&1
ApplyVQSR
代码语言：javascript复制 ## 然后就是对 SNP 运行 ApplyVQSR
${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" ApplyVQSR
-V 6.gvcf/germline_merge.vcf.gz
--recal-file 6.gvcf/snp.recal
--tranches-file 6.gvcf/snp.tranches
--truth-sensitivity-filter-level 99.0
--create-output-variant-index true
-mode SNP
-O 6.gvcf/snp.vqsr.vcf.gz
1>./6.gvcf/SNP_ApplyVQSR.log 2>&1
## 对 INDEL 运行 ApplyVQSR，注意，这里输入的时候，采用上一步校正过 SNP 质量后得到的结果
${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" ApplyVQSR
-V 6.gvcf/snp.vqsr.vcf.gz
--recal-file 6.gvcf/indel.recal
--tranches-file 6.gvcf/indel.tranches
--truth-sensitivity-filter-level 99.0
--create-output-variant-index true
-mode INDEL
-O 6.gvcf/snp.indel.vqsr.vcf.gz
1>>./6.gvcf/INDEL_ApplyVQSR.log 2>&1
最后经过过滤的位点在 FILTER 会标记上 PASS 标签。
参考： https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels-
