癌症分析【第4课】GATK Germline变异检测

第四课：GATK Germline 变异检测

当你想要用GATK做分析，当然要看工具文档了。来看看官网上的Tool Documentation是这样滴：

内容很多，就Short Variant Discovery而言就有7个工具。要全部看完才能用？估计就从入门到放弃了……其实，作为一个生信自学者，要先确定几个问题：

• 我的目的是什么？ 需要做什么分析？
• 需要用哪些工具？先用哪个，再用哪个？
• 每一个工具具体的用法、参数设置是什么？

这些问题的答案，都在GATK官网里有。这一课，我们来充分利用一下GATK官网提供的资源吧！

Q：我想用GATK做Germline变异检测，怎么做呢？

A：GATK官网，参考Best Practices提供的流程，用自己的数据体验一遍吧！上一课生成的经过质量值校正的bam文件（bqsr.bam）在这里就可以用上了，为了节省时间，这一课也继续使用chr6.fasta文件作为reference参考基因组文件。

Germline short variant discovery (SNPs + Indels) 
https://software.broadinstitute.org/gatk/best-practices/workflow?id=11145

目前我们需要解决的问题是，对于单个样本的变异检测，即单个样本的Germline SNPs + Indels的检测，在GATK官网Best Practices提供了上述流程。由图可见，主要分为三个步骤：

1. Variant Discovery 变异的发现
   HaplotypeCaller
2. 变异检测Callset Refinement 变异的后处理
   CNNScoreVariants 用模型给变异打分
   FilterVariantTranches
3. 根据打分注释过滤变异Functional Annotation 变异的功能注释 
   Funcotator 变异功能注释

接下来，我们来动手试一下吧！

Step1 Variant Discovery

工具：HaplotypeCaller

GATK4 HaplotypeCaller在可能存在变异的区域，采用局部de-novo组装单倍型的方法，同时找出SNP和InDel变异：

gatk --java-options "-Xmx4g" HaplotypeCaller -R chr6.fasta -I bqsr.bam -O sample.vcf.gz \

注：由于我们只有一个样本，因此只用默认参数，输出VCF格式的变异结果。若你要分析不止一个样本，Best Practices建议使用GVCF模式（参数 -ERC GVCF）。GVCF中的G代表了Genomic，实际上是VCF格式的一种，只是包含了比常规的VCF格式更多的信息。GVCF格式中记录了所有的位点，无论这些位点上有没有发现变异。在多个样本的分析中，采用GenomicsDBImport和GenotypeGVCFs将多个样本的变异文件进行合并分析，可提高变异检测的敏感性。

What is a GVCF and how is it different from a 'regular' VCF?
https://software.broadinstitute.org/gatk/documentation/article.php?id=4017

思考题：如果你有多个样本，可以尝试根据Best Practices的步骤，进行群体（多个样本）的变异检测分析（含单个样本的GVCF生成、多个样本的GVCF合并及联合分析、VQSR过滤及功能注释等步骤）。

Step2 Callset Refinement

工具：CNNScoreVariants和FilterVariantTranches

完成HaplotypeCaller变异检测后，需要对找到的变异集（callset）进行过滤。CNNScoreVariant用已经训练好的CNN模型对变异集进行打分，在用FilterVariantTranches根据打分对变异进行过滤。

CNNScoreVariant有已经预训练好的1D Model和2D Model，本次演示采用1D Model（如果你想试一下2D Model，则还需要输入bam文件，运行需要比较长的时间，别忘了给你的电脑插上电源哦）。

1D Model with pre-trained architecture

gatk CNNScoreVariants -V sample.vcf.gz -R chr6.fasta -O annotated.vcf

注：如果你使用本机安装的GATK，可能使用CNNScoreVariants会出现报错，这是因为该工具的使用需要预先配置运行环境。此处建议使用Docker版本的GATK，因此Docker版本的GATK已经做好了配置，可直接使用。

Apply tranche filters based on the scores in the info field with key CNN_1D

gatk FilterVariantTranches \
-V annotated.vcf \
--resource dbsnp_138.hg19.vcf.gz \
--resource 1000G_phase1.indels.hg19.sites.vcf.gz \
--info-key CNN_1D \
--snp-tranche 99.9 --snp-tranche 99.0 --snp-tranche 95.0 \
--indel-tranche 99.9 --indel-tranche 99.0 --indel-tranche 95.0 \
-O filtered.vcf

变异集过滤完成

思考题：还记得从哪里下载变异数据集吗？提示：GATK官网资源Download, Resource Bundle。

Step3 Functional Annotation

工具：Funcotator
对过滤后的变异集进行功能注释，需先下载注释数据资源（Data Source），Germline的注释数据资源并不大（约14M），可以在Docker容器中，在/gatk目录下使用以下命令下载：

Download Data Source (germline)

gatk FuncotatorDataSourceDownloader \
--germline \
--validate-integrity \
--extract-after-download

下载的注释资源存在Docker容器的/gatk/funcotator_dataSources.v1.6.20190124g/目录下，在下一步功能注释中作为参数输入：

Function Annotation

gatk Funcotator \
-R chr6.fasta \
-V filtered.vcf \
-O filtered.maf \
--output-file-format MAF \
--data-sources-path /gatk/funcotator_dataSources.v1.6.20190124g/ \
--ref-version hg19

功能注释结果文件

温馨提示

本次课程中每一步产生的结果文件均存放在如下地址，可下载且永久有效：  ftp://ftp.cngb.org/pub/course/Hands-on-Lab/Lesson4/

本次小课堂就到这里，下次再见~