GATK Best Practices:通过GATK4 docker运行processing-for-variant-discovery-gatk4.wdl

Run GATK Best Practices for data pre-processing by Cromwell/WDL

与GATK4正式发布的还有WDL(workflow description langaue,https://software.broadinstitute.org/wdl/),WDL将工作流程分为了workflow, task, call, command 和 output。与以往GATK提供Best practice的PPT介绍不同,现在Broad提供的是Best practice(https://software.broadinstitute.org/gatk/best-practices/)的WDL文件。WDL文件运行通过cromwell运行,并且有json格式的参数输入文件指定WDL文件中流程所需要的参数。比如

sudo java -jar cromwell.jar run workflow.wdl --inputs workflow.inputs.json

我们只需要修改json文件中的参数就可以运行gatk4 Best Practices,而不需要自己去搭建流程,简化了工作量,也遵循了Broad提供的推荐设置和流程。本文只介绍突变检测前的序列比对和recalibrate这部分的GATK best practices,该流程生成了用于variant calling的bam文件。

1,文件准备

WDL文件和json文件
Broad在github上提供了进行突变检测call variant之前的数据处理data proceesing流程,见https://github.com/gatk-workflows/gatk4-data-processing
从github上,我们需要下载两个文件
processing-for-variant-discovery-gatk4.wdl (用于data pre-processing 的 pipeline)
processing-for-variant-discovery-gatk4.hg38.wgs.inputs.json(指定WDL的参数文件)

ubam文件:
要求ubam文件中要有RG tag,经过排序sort之后,该文件可以通过picard将fastq文件转换得到

GATK resoure bundle,从中下载GATK需要的dbsnp文件,known site等文件
ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/
https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0/

2, 软件准备

下载cromwell,用于运行WDL
官网:http://cromwell.readthedocs.io/en/develop/tutorials/FiveMinuteIntro/

下载GATK docker镜像,根据json文件中的第5部分DOCKERS,还需要下载python和genomes in the cloud 镜像,pull命令中加入registry.docker-cn.com,可以指定docker从国内源中下载,避免下载速度较慢或不能下载。

sudo docker pull registry.docker-cn.com/broadinstitute/gatk:latest
sudo docker pull registry.docker-cn.com/library/python:2.7
sudo docker pull registry.docker-cn.com/broadinstitute/genomes-in-the-cloud:2.3.1-1512499786

3,修改json文件

json文件中的gs://路径下的文件是google storage cloud上面的文件,可以换成本地电脑的路径。

第一部分 SAMPLE NAME AND UNMAPPED BAMS:
该部分指定样本名和分析需要的起始ubam文件,该部分会经常修改

sample_name 样本名
ref_name 基因组版本
flowcell_unmapped_bams_list ubam的list文件路径,该路径中每一行为一个ubam的绝对路径
unmapped_bam_suffix ubam文件后缀

第二部分 REFERENCE FILES:
修改成参考基因组的fastq文件、dict文件和用bwa index之后生成的alt、ann、bwt文件等路径。该部分只需要修改一次,后续不用在修改。

第三部分 KNOWN SITES RESOURCES:
从GATK resource bundle中下载的文件,按照要求改成本地路径即可。

第四部分 MISC PARAMETERS
不用修改。因为我们下载的docker镜像中含有预先编译好的gatk,bwa,picard等程序

第五部分 DOCKERS
如果不确定,可以通过sudo docker images查看
REPOSITORY TAG IMAGE ID CREATED VIRTUAL SIZE
registry.docker-cn.com/library/python 2.7 b1b87432bf0e 42 hours ago 678.2 MB
registry.docker-cn.com/broadinstitute/gatk latest 5e9cb2ded24b 2 weeks ago 4.882 GB
registry.docker-cn.com/broadinstitute/genomes-in-the-cloud 2.3.1-1512499786 989127e46444 3 months ago 1.879 GB
在json文件中将对应的docker镜像改成本地有的,切记tag要写对。

第六部分 PATHS
ls
不用修改,broadinstitute/genomes-in-the-cloud的iamge镜像生成容器中的路径,
如果使用sudo docker run -it 进入broadinstitute/genomes-in-the-cloud的话,可以看到/usr/gitc/路径中的其他文件
同样,/gatk/gatk是broadinstitute/gatk的image镜像生成的容器中的路径,所以不用修改。

第七到十部分也不用修改。

因为cromwell运行WDL的时候会自动调用docker,所以我们不需要关注docker怎么使用,甚至不用关系流程,因为都是broad提供,我们只需提供json文件即可。如果想在docker上跑gatk4,可以参考https://gatkforums.broadinstitute.org/gatk/discussion/10172/how-to-run-the-gatk4-docker-locally-and-take-a-look-inside

4,运行GATK4 Best practices中的data preprocessing流程

sudo java -jar cromwell.jar run processing-for-variant-discovery-gatk4.wdl --inputs processing-for-variant-discovery-gatk4.hg38.wgs.inputs.json

刚开始时会生成一个workflow的id,中间文件都会在该文件夹内。

如果中间运行出错,可以查看workflow文件夹中的stderr文件。cromwell的文件夹结构如下

运行结束后,会收到
[2018-03-22 19:04:55,04] [info] SingleWorkflowRunnerActor workflow finished with status ‘Succeeded’.的提示。

就这样,我在不了解docker的情况下,在没有关注GATK Best Practices具体流程的情况下,对ubam进行比对,mark duplicate,recalibrate等操作生成了用于variant calling的bam文件。具体的Best practices流程,需要通过研究processing-for-variant-discovery-gatk4.wdl文件。后续在分析其他样本时,只需要将json文件中的第一部分修改即可。

下一步就是通过WDL进行germline variant calling的工作

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

GATK Best Practices:通过GATK4 docker运行processing-for-variant-discovery-gatk4.wdl》上有13条评论

  1. Karen

    你好,我在用cromwell想在firecloud上面直接运行gatk的wdl,json里面的内容是不是只需要修改input部分就可以呢?自己瞎研究了半天没太弄明白,所以想来请教一下您

    回复
    1. zzx 文章作者

      应该是修改input部分就行了,因为我是本地跑的,把resource部分也改了。我和同事本地跑的很麻烦,改了很多东西才行。如果用GATK提供的wdl文件,建议在云上进行。

      回复
  2. xiaoyuan

    你你好,感谢您的关于GATK数据处理过程经验的分享,但是我还有几个问题想要请教您。一、如果做其他物种的SNPs计算,需要将GATK提供最佳实践里的wdl文件里的这段代码里的hg38名称需要修改吗?二、如果有多个个体,应用最佳实践可以批量处理吗?

    回复
    1. zzx 文章作者

      其他物种的不好做variant recalibrate,因为没有用于训练的数据集。我现在不建议wdl,感觉不适合本地跑,要改的东西也多。建议用命令行来做。

      回复
        1. zzx 文章作者

          以前GATK老版本的时候,官网会有best practice的教程,和这个上面的PPT差不多,告诉你怎么做。https://www.broadinstitute.org/partnerships/education/broade/best-practices-variant-calling-gatk-1
          大概分为比对->indel realignment(GATK4移除了)->base recalibrate->variant recalibrate->call。如果是多个样本的话,可以先单独生成gvcf,然后进行joint call

          然后你根据PPT的内容,可以去gatk论坛上搜,论坛上有教程,比如variant recalibrate的教程是下面的这个
          https://gatkforums.broadinstitute.org/gatk/discussion/2805/howto-recalibrate-variant-quality-scores-run-vqsr

          另外就是有些公众号他们也会发布命令行教程,但大部分都是以人重测序为主的。

          回复
  3. Angry_Panda

    生物生,gatk新手,请教一个问题关于你提到的:flowcell_unmapped_bams_list是ubam的list文件路径,该路径中每一行为一个ubam的绝对路径。
    请问是理解为,做一个txt文件,里面每行列出各个ubam文件的绝对路径吗
    请问可以简化如下吗,把所有以上ubam文件放在一个文件夹,直接写那个文件夹的绝对路径。

    回复
    1. zzx 文章作者

      一个txt文件,里面每一行是一个ubam的路径,不能简化。最好不用docker,用docker太麻烦了。可以参考这个,https://www.broadinstitute.org/partnerships/education/broade/best-practices-variant-calling-gatk-1,或者百度下GATK的教程,真的解决不了的话,我也可以提供我们的流程。

      回复
      1. Angry_Panda

        十分感谢你的回复.我又试了几天,仍然是一直在报错. 加之计算机知识比较薄弱, 弄的我十分没有信心了.我知道要流程是一个挺不好的行为,但是现在我有点无从下手了. 不知道能否分享你用的本地跑的流程给我看看? 哪怕只是Best practice的一个小部分,但是我可以从中学习到底哪一部分是应该怎么修改,对应什么.
        不论能否分享给我看看,都先跟你道一声谢谢, 已经从你的帖子学习到了一些东西. 只是自己这块太差了吧,哎.

        回复
  4. shuang

    你好,谢谢你的文章,十分实用。请问你有不用docker 的版本吗,我现在公司要求我做一个不用docker的出来。我自己是小白,实在是没理解那些docker thing是指代的什么。有很多 docker thing in .json and .wdl file, such as “##_COMMENT5”: “DOCKERS” in json file and runtime part in wdl file. 我不知道怎么修改他们。
    十分感谢你的回复和帮助
    br,shuang

    回复
  5. xiaoyun

    您好, 谢谢分享, 请问你文中写道如下
    ubam文件:
    要求ubam文件中要有RG tag,经过排序sort之后,该文件可以通过picard将fastq文件转换得到.

    我准备重复你的实验, 我打开json file 里面的txt file从中找到了24个ubam file谷歌云地址,分别下载下来, 是可以直接作为input跑程序了吗? 还是说我需要额外的操作. 你上文那句话,我没有看懂 感觉弄混了.

    回复
    1. zzx 文章作者

      如果自己测序的话,得到的是fastq文件,需要用picard转成ubam进行下一步分析,你要是有ubam的话,就可以直接进行。

      回复

发表评论

电子邮件地址不会被公开。 必填项已用*标注

This site uses Akismet to reduce spam. Learn how your comment data is processed.