GATK 4.0 WGS germline call variant

GATK升级4.0版了，作为人类call variant的金标准软件，加上其强大的团队，每次重大更新都会给使用者带来一点新的东西（或者说是改变），我也正好整理下，将GATK基本分析流程过渡到4.0版本

GATK4.0最明显的变化是其命令调用发生了改变，可以看看这个就明白了 https://software.broadinstitute.org/gatk/documentation/quickstart

GATK4.0还集成了picard工具以及增加了SortSam功能，这样Germline short variant discovery流程只要GATK一个软件即可完成sam/bam文件到vcf文件全流程

整体上Germline call variant分析思路没有变，所以我按照 Germline short variant discovery (SNPs + Indels) 教程上所示流程图，将GATK4.0基本用法再过一遍，主要是想再对GATK的VQSR质控理解下（之前一直没好好看其说明）

GATK call variant（snp/indel）

测试数据KPGP，下载地址： ftp://ftp.kobic.re.kr/pub/KPGP/2017_release_candidate_update_only/WGS/KPGP-00216/

参考基因组：GATK的hg38文件：

nohup wget wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.fai &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.dict &

使用fastp对测序数据质控，非常好用的一个质控软件，海普洛斯开发的

~/biosoft/fastp/fastp -i KPGP-00216_L1_R1.fq -I KPGP-00216_L1_R2.fq -o KPGP-00216_L1_R1.clean.fq -O KPGP-00216_L1_R2.clean.fq

可以看看质控报告fastp.html，蛮好的用的

压缩下原始数据，节约下空间，暂时不删（删了也可以）

tar cf - KPGP-00216_L1_R1.fq |pigz >KPGP-00216_L1_R1.fq.tar.gz

bwa比对

1. 先建bwa的索引

   bwa index -a bwtsw -p gatk_hg38 ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta &

2. 比对

   nohup bwa mem -t 4 -R '@RG/tID:KPGP00216/tPL:illumina/tLB:WGS/tSM:KPGP00216' ~/reference/index/bwa/gatk_hg38/gatk_hg38 KPGP-00216_L1_R1.clean.fq KPGP-00216_L1_R2.clean.fq 1>KPGP-00216_L1.sam 2>/dev/null &

将sam格式转化为bam格式，并且排序

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" SortSam -I KPGP-00216_L1.sam -O KPGP-00216_L1.sorted.bam -SO coordinate --CREATE_INDEX true

标记PCR重复序列

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" MarkDuplicates -I KPGP-00216_L1.sorted.bam -O KPGP-00216_L1.sorted.marked.bam -M KPGP-00216_L1.metrics

这个自己理解下哈，GATK里并没有建议一定要做

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" FixMateInformation -I KPGP-00216_L1.sorted.marked.bam -O KPGP-00216_L1.sorted.marked.fixed.bam -SO coordinate --CREATE_INDEX true

对参考基因组fa生成dict文件（其实在下载基因组时就该将这个dict文件一起下载下来）

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk CreateSequenceDictionary -R Homo_sapiens_assembly38.fasta -O Homo_sapiens_assembly38.dict

BQSR，也就是Recalibration Base Quality Score

1. 先用FTP链接下GATK的resource database

   主机：ftp.broadinstitute.org  
   用户：gsapubftp-anonymous

2. 下载BQSR需要的注释位点VCF文件

   wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz
   wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi
   wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
   wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi
   wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
   wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi

3. 运行BQSR，先BaseRecalibrator生成一个recal.data.table中间文件，再运行ApplyBQSR

   ~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" BaseRecalibrator -R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta -I KPGP-00216_L1.sorted.marked.fixed.bam -O recal_data.table --known-sites ~/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz --known-sites ~/annotation/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --known-sites ~/annotation/GATK/resources/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz

接下来则是使用HaplotypeCaller对上述的bam文件进行call variant，就是寻找变异的过程，生成一个VCF文件用于后续位点的质控和注释，最简单的单样本call variant如命令所示，从一个bam到一个vcf文件

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G HaplotypeCaller /
    -R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
    -I KPGP-00216_L1.sorted.marked.fixed.BQSR.bam /
    --dbsnp ~/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz /
    -O KPGP-00216_L1.raw.vcf

按照我的电脑配置，用了21h(组装的台式机，非服务器~)，为了加快速度，可以采用下述分染色体来call variant(也就是一个染色体生成一个VCF文件，22条常染色体 + XY性染色体，总共24个VCF文件)，当然不是一个个运行，使用shell命令行for循环加上 & 提交后台运行，如下所示：

for i in {1..22} X Y;do (~/biosoft/GATK4.0/gatk-4.0.5.1/gatk HaplotypeCaller -R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta -I KPGP-00216_L1.sorted.marked.fixed.BQSR.bam --dbsnp ~/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz -L chr$i -O KPGP-00216_L1.chr$i.raw.vcf &);done

然后再用GatherVcfs将多个染色体的vcf文件合并

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk GatherVcfs $(for i in {1..22} X Y;do echo "-I KPGP-00216_L1.chr${i}.raw.vcf";done) -O KPGP-00216_L1_gather.raw.vcf

如果你有多个样本（或者大样本的一个队列），那么首选用g.vcf，其是GTAK HaplotypeCaller生成第一个中间文件，可以用于后续的GenotypeGVCFs进行joint genotyping；g.vcf简单的说是储存了所有位点的变异信息（不做任何可信度等的过滤）的vcf文件；其使用方式就上述HaplotypeCaller工具上加上 -ERC GVCF ，然后输出是XXX.g.vcf

有一点需要注意的是：

按照以前的思路，将多个样本的g.vcf文件一起通过的GenotypeGVCFs进行joint genotyping；但是！GATK 4.0的GenotypeGVCFs只支持a single single-sample GVCF，a single multi-sample GVCF created by CombineGVCFs 以及a GenomicsDB workspace created by GenomicsDBImport；所以之前的方法已经失效了，你在用GenotypeGVCFs前需要将多个样本的g.vcf文件用CombineGVCFs方式或者GenomicsDBImport方式合并成一个文件，前者（比较传统）是一个总的g.vcf文件，后者是一个GenomicsDB（XX.db）

这里对这个先不做详细说明啦，后续笔记会做尝试

另外对于多样本是分开cal还是joint call，可以参考下： 生信笔记：call snp是应该一起call还是分开call？

GTAK的VQSR质控

接下来则的对获得的VCF文件做质控，对于人物种的WGS测序，首选VQSR（Variant Quality Score Recalibration），其原理可以看 Variant Quality Score Recalibration (VQSR)-High-level overview ，个人理解简单的说：

1. GATK认为VQSR比根据各种annotations进行hard-filtering过滤要好，减少了人为阈值的局限性，避免了一刀切的弊端，从而在sensitivity和specificity之间达到一定的平衡
2. VQSR根据机器学习算法从highly validated变异位点数据集（每个位点的annotation profile，一般使用5-8个annotation）中获取到good variants/bad variants
3. 根据上述的位点从我们自己数据集中挑选出一个变异子集（probably true positives）来建模训练，获得一个可识别good variants的模型；bad variants的模型也是如此获得
4. 然后根据上述获得的模型，对自己数据集的每个变异位点进行一个总的打分
5. 最后根据设定的sensitivity阈值对变异位点进行过滤

更简单的说：

利用自身数据和已知变异位点集的overlap，通过GMM模型构建一个分类器来对变异数据进行打分，从而评估每个位点的可行度

如果官网文档原理看的不是很清楚，可以参考这个（个人觉得讲的很好） GATK4.0和全基因组数据分析实践（下） 中的VQSR原理部分

所以可看出，VQSR对你的数据集的变异位点的数目有一定的要求；如果你的变异位点数目过少（如样本数较少的WES或者小panel测序），那么是无法做VQSR的；GATK对于样本不达标的WES，建议用1000 Genomes Project中的数据代替（将1000G的bam文件和WES样本的bam生成g.vcf文件，再一起做joint call）

GTAK的VQSR主要分两步骤进行：

• VariantRecalibrator 对应我上述理解原理的1-4步骤，对每个变异位点打分注释VQSLOD value，从而生成一个recalibration文件以及一个xxx.plots.R.pdf(Gaussian mixture model plots)
• ApplyRecalibration 对应5步骤，根据recalibration文件生成recalibrated VCF文件，并且根据过滤参数进行过滤（标记PASS）

所以当你tranche sensitivity（–truth-sensitivity-filter-level）设置为99.9%时，则表示将VQSLOD值高于整体99.9%的变异位点标记为PASS，表示通过过滤阈值，剩下的位点则被认为是假阳性的了；这个参数可以根据你的具体需求来设定，看你是需要more specific or more sensitive（tranche从90到100，specific随之降低，而sensitive随之升高）；或者生成多个tranche的结果，从中挑选你满意的阈值（可以看结果文件SNP中的xxx.tranches，而indel是没这个文件的；人全基因组Ti/Tv（转换和颠换的比例）的值一般在2.0-2.1…）

GATK的VQSR用的annotations指标有以下几种（VQSR的文档中还提到了InbreedingCoeff，这个是在大于10个样本中才使用的一个指标；如果样本过少或者是closely related samples(such as a family)的话，建议剔除；而且QD不适合用于外显子测序的数据的VQSR）：

• Coverage（DP）
• QualByDepth（QD）
• FisherStrand (FS)
• StrandOddsRatio (SOR)
• RMSMappingQuality (MQ)
• MappingQualityRankSumTest (MQRankSum)
• ReadPosRankSumTest (ReadPosRankSum)

按照官方教程，SNP的VQSR过滤，选用的resource datasets为：

• HapMap，hapmap_3.3.hg38.vcf.gz， truth=true 表示VQSR将这个数据集中的变异位点作为真位点true sites， training=true 表示VQSR将true sites用于训练recalibration model，并赋予这些变异位点prior likelihood值为 Q15 (96.84%)
• Omni，1000G_omni2.5.hg38.vcf.gz， truth=true ， training=false （文档中写着是true，参数建议中写着的是false。。。我就按照参数上的来了）， Q12 (93.69%)
• 1000G，1000G_phase1.snps.high_confidence.hg38.vcf.gz， truth=false 表示VQSR考虑到在1000G数据集中的不仅包含了true variants还有false positives， training=true ， Q10 (90%)
• dbSNP，dbsnp_146.hg38.vcf.gz， truth=false 表示VQSR未将dbSNP数据集中的位点作为可信数据集， training=false 表示不用于训练数据集， known=true 表示stratify output metrics such as Ti/Tv ratio by whether variants are present in dbsnp or not， Q2 (36.90%)

INDEL的VQSR过滤，选用的resource datasets为：

• Mills，Mills_and_1000G_gold_standard.indels.hg38.vcf.gz， truth=true ， training=true ， Q12 (93.69%)
• dbSNP，dbsnp_146.hg38.vcf.gz， truth=false ， training=false ， known=true ， Q2 (36.90%)

其中几个resource datasets已经在BQSR中下载了，剩下的也是按照其方法下载即可

对于GATK的VQSR，snp和indel过滤是分开的，两者参数整体上差不多，但有点略微的区别

SNP的VQSR的过滤，一般的参数可以如下：

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk VariantRecalibrator /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-V KPGP-00216_L1_gather.raw.vcf /
--resource hapmap,known=false,training=true,truth=true,prior=15.0:/home/anlan/annotation/GATK/resources/bundle/hg38/hapmap_3.3.hg38.vcf.gz /
--resource omni,known=false,training=true,truth=false,prior=12.0:/home/anlan/annotation/GATK/resources/bundle/hg38/1000G_omni2.5.hg38.vcf.gz /
--resource 1000G,known=false,training=true,truth=false,prior=10.0:/home/anlan/annotation/GATK/resources/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz /
--resource dbsnp,known=true,training=false,truth=false,prior=2.0:/home/anlan/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz /
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP /
-mode SNP /
-O KPGP-00216_L1.snp.recal /
--tranches-file KPGP-00216_L1.snp.tranches /
--rscript-file KPGP-00216_L1.snp.plots.R

这里还有个参数要提一下， -tranche 默认是输出[100,99.9,99.0,90.0]4个tranche阈值的统计结果，如果想看其他阈值的结果，需要自行加上；结果就是看KPGP-00216_L1.snp.tranches（还有图形展示的KPGP-00216_L1.snp.tranches.pdf），而KPGP-00216_L1.snp.recal文件则是用于ApplyVQSR的

Applying recalibration/filtering to SNPs， --truth-sensitivity-filter-level 可以根据自身需求来设定

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk ApplyVQSR /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-V KPGP-00216_L1_gather.raw.vcf /
-O KPGP-00216_L1.snp.VQSR.vcf /
--truth-sensitivity-filter-level 99.5 /
--tranches-file KPGP-00216_L1.snp.tranches /
--recal-file KPGP-00216_L1.snp.recal /
-mode SNP

INDEL的VQSR过滤，这里 -an 参数跟SNP有点不同，没有MQ，并且加了 --max-gaussians 4 用于设定Gaussians（clusters of variants that have similar properties）的数目，即减少聚类的组数，从而使得每个组的变异位点数目达到要求

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk VariantRecalibrator /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-V KPGP-00216_L1_gather.raw.vcf /
--max-gaussians 4 /
--resource mills,known=false,training=true,truth=true,prior=12.0:/home/anlan/annotation/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz /
--resource dbsnp,known=true,training=false,truth=false,prior=2.0:/home/anlan/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz /
-an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum /
-mode INDEL /
-O KPGP-00216_L1.indel.recal /
--tranches-file KPGP-00216_L1.indel.tranches /
--rscript-file KPGP-00216_L1.indel.plots.R

Applying recalibration/filtering to INDEL

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk ApplyVQSR /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-V KPGP-00216_L1_gather.raw.vcf /
-O KPGP-00216_L1.indel.VQSR.vcf /
--truth-sensitivity-filter-level 99.0 /
--tranches-file KPGP-00216_L1.indel.tranches /
--recal-file KPGP-00216_L1.indel.recal /
-mode INDEL

如果VQSR不能使用的话，那么只能用硬过滤（hard-filtering）了，这个注释指标每个人都可能用的不一样，可以先从GATK官方建议的指标入手，然后再看看文献中应用的指标，最后再结合自身需要来设定

以上是个人对GATK的一些理解，有些理解的不够好的地方有请联系我纠正下（尤其的VQSR那块，有些地方还云里雾里的）

本文出自于 http://www.bioinfo-scrounger.com 转载请注明出处