GATK4全基因组数据分析最佳实践

本文转载自微信公众号解螺旋的矿工，作者为黄树嘉，已获得授权。黄树嘉写了WGS系列的文章，堪称教科书级别的生物信息学习材料。我们将陆续转载给大家。大家也可以关注公众号解螺旋的矿工。

这是我根据之前的WGS系列和GATK4实践文章进行重新梳理之后确定下来的分析流程，这是一个WGS的最佳实践，它基于GATK4和我的实际经验，稍作修改即可应用到实际的项目中。我以这篇文章为标记，终结当前WGS系列数据分析的主体流程问题。

首先，要提醒大家的是，我在这里所呈现的仅仅只是一个shell流程，它并不符合正规的IT设计规范，你也需要对其参数做适当的修改，比如修改相关软件路径和数据路径，放在这里主要有以下两个目的：

• 第一、我想清楚地告知大家一个完整的WGS数据分析流程里面到底都有啥，执行顺序是怎么样的，先给出一个直接可用的版本供参考，大家在有需要的时候，就不用从零开始了，有据可依，少走一些弯路，可以结合实际项目和流程中的注释信息稍作修改即可使用，提高效率；

• 第二、打下一个最基本的WGS流程基准，为以后在这个基础之上跌代修改提供依据。

另外，GATK4虽然已经正式发布了一段时间，用法也确实和GATK3有些不同，但是差距很小，主要是参数形式和调用形式上的改变（我之前的两篇GATK4实践文章对此已有过演示）。它主要是整合了多个现有工具、新增一部分新功能（包括Spark功能），但基本的核心内容并没有改变。

之前我听说有些同学觉得改换了GATK4之后，担心原来基于GATK3的WGS分析方法就不适用了。这其实有些多虑了，WGS数据分析和解读是我们的目的，GATK是帮助我们达成该目的的一个重要工具，但如果有更好/更合适的工具，我们随时准备替换，甚至重写。所以GATK也好，BWA也罢，对于我们而言都只是“术”，重要的是，我们要知道该如何对数据进行分析和解读，这是根本之“道”。这也是我在写作过程中努力坚持的一个宗旨和原则。

另外，我在以下流程的注释中留下了很多重要的信息，以及一些步骤的使用条件，例如某些步骤（比如HaplotypeCaller）可以有多种不同的实现路径，这个可以根据实际的需要进行选择。

好了，现在进入正篇。

我之前已经用四五万字写了一系列的WGS文章，对其中的很多原理和原则都做了详细的解释，这里就不再赘述了，如有需要可以在本文下方的推荐阅读中查阅。另外，在这个流程中我做了一些默认设置，比如，参考序列和GATK4所需的bundle数据都是选择了目前最新的hg38。

那么，接下来，就直接上可用的流程代码吧，我设计为两种模式：

第一：单样本模式

你可以把下面这一段代码，写到一个shell脚本中（注意代码中的“\”是接下一行的符号，它后面不能有任何空格或者其它字符），我们这里把文件名定为，wgs_single.sh。这个流程假设你只有一个样本，这个样本只有一对Illumina测序的PE fastq数据文件。流程有6个参数，你可以在代码中清楚地看到，分别是：

• read1的路径

• read2的路径

• Read Group ID

• 测序文库编号

• 样本编号

• 输出目录路径

用起来也比较简单，直接在命令行中执行，并接上以上对应的参数即可，比如：

$ sh wgs_single.sh read.1.fq.gz read.2.fq.gz Test_RG Test_lib Test_sample Test_outdir

以下是完整的流程代码：

第二：多样本模式

以上单样本的流程比较简单直接，如果你没有做成一键式产品的需要，那么基本上是可以用上面的流程一步步完成你的WGS分析的。但随着目前测序行业的发展，大规模人群的测序会变得越来越普遍，并且基于群体的数据分析要远优于单样本，意义也更加深远，因此这种多样本的模式将成为常态。

多样本的流程和单样本相比会有些不同，首先是不适合把执行过程都封装在同一个shell脚本中。在这里我提供了一种实现方式，可供参考，我分为两步：

• 第一步，单独为每个样本生成后续分析所需的中间文件——gVCF文件。这一步中包含了对原始fastq数据的质控、比对、排序、标记重复序列、BQSR和HaplotypeCaller gVCF等过程。这些过程全部都适合在单样本维度下独立完成。值得注意的是，与单样本模式不同，该模式中每个样本的gVCF应该成为这类流程的标配，在后续的步骤中我们可以通过gVCF很方便地完成群体的Joint Calling；

• 第二步，依据第一步完成的gVCF对这个群体进行Joint Calling，从而得到这个群体的变异结果和每个人准确的基因型（Genotype），最后使用VQSR完成变异的质控。这两个步骤其实还包含了许多细节，具体可见我在流程中的注释。

以下是第一步的代码，参数和单样本模式一样。这个代码我们可以把它存放在一个名为wgs_fastq_to_gvcf.sh的文件中：

接着，这是第二步的代码，我们把它存放在一个名为wgs_gvcf_to_vcf.sh的文件中。需要强调的是，它的输入和输出参数需要与第一步保持一致，流程中我对此做了额外的注释：

使用方式可以参考上文“单样本模式”的例子，直接在命令行中完成即可，不再赘述。

变异注释

变异的注释这一个步骤，我把它单独拎出来，原因是它完全可以独立于以上的流程。任何SNP和Indel数据（VCF格式），只要有需要你都可以随时完成这个注释。由于比较简单，流程的细节我在这里就不再多说了，最难的其实只是如何安装好VEP和它需要的相关数据集（cachedir目录下的数据）。

## 使用VEP完成变异的注释
VEP=/your_path_to/ensembl-vep/vep
time $VEP --fasta $reference/Homo_sapiens_assembly38.fasta \
 --vcf --merged --fork 10 --hgvs --force_overwrite --everything \
   --offline --dir_cache /your_path_to/ensembl-vep/cachedir \
   -i $outdir/gatk/${sample}.HC.VQSR.vcf.gz \
   -o $outdir/gatk/${sample}.HC.VQSR.VEP.vcf.gz

小结

好了，这篇文章就到此结束了，祝你数据分析顺利~

※  ※  ※

---探索生命科学的技术---

生信圈教科书式必读公众号

推荐阅读

这是我的知识星球：解螺旋技术交流圈，是一个与读者朋友们的私人朋友圈，欢迎你的加入。我有9年前沿而完整的生物信息学、NGS领域的工作经历，在该领域发有多篇Nature级别的科学文章。

这是知识星球上第一个真正与基因组学和生物信息学强相关的圈子。我旨在营造一个高质量的组学知识圈和人脉圈，通过提问、彼此分享、交流经验、心得等，彼此更好地学习生信知识，提升基因组数据分析和解读的能力。

在这里你可以结识到全国优秀的基因组学和生物信息学专家，同时可以分享你的经验、见解和思考，有问题也可以向我提问和圈里的星友们提问。