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
  • 测序文库编号
  • 样本编号
  • 输出目录路径

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

1
$ 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目录下的数据)。

1
2
3
4
5
6
7
## 使用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级别的科学文章。

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

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

知识星球邀请链接:「解螺旋技术交流圈」