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比对

  • 先建bwa的索引

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

    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

  • 先用FTP链接下GATK的resource database

    主机:ftp.broadinstitute.org
    用户:gsapubftp-anonymous
    
  • 下载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
    
  • 运行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
    
    ~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" ApplyBQSR \
    -R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta \
    -I KPGP-00216_L1.sorted.marked.fixed.bam \
    --bqsr-recal-file recal_data.table \
    -O KPGP-00216_L1.sorted.marked.fixed.BQSR.bam
    

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

    ~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" 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=truetraining=false(文档中写着是true,参数建议中写着的是false。。。我就按照参数上的来了),Q12 (93.69%)
  • 1000G,1000G_phase1.snps.high_confidence.hg38.vcf.gz,truth=false表示VQSR考虑到在1000G数据集中的不仅包含了true variants还有false positives,training=trueQ10 (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=truetraining=trueQ12 (93.69%)
  • dbSNP,dbsnp_146.hg38.vcf.gz,truth=falsetraining=falseknown=trueQ2 (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那块,有些地方还云里雾里的)

参考资料: Variant Quality Score Recalibration (VQSR)
(howto) Recalibrate variant quality scores = run VQSR
GATK流程
(How to) Consolidate GVCFs for joint calling with GenotypeGVCFs
Germline short variant discovery (SNPs + Indels)

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