GATK best practice每个步骤耗时如何?

上次我们介绍了完整的 GATK best practice(请点击) 在我的基因组重测续数据分析流程,详细讲解了每个步骤的代码,输入输出文件,准备文件,以及耗时。 但是对同一个样本的多个lane的数据合并的问题,缺失了一个重要步骤,而且有热心的读者咨询整个流程的耗时问题,所以特出此番外篇作为补充。

再次介绍一下我的这次基因组重测续数据共5条lane,都是单独的PE150的fastq文件。

如果仅对L1这条lane数据进行处理

首先是BWA的比对耗时如下:

[main] Real time: 15870.794 sec; CPU: 77463.156 secpicard.sam.SortSam done. Elapsed time: 45.60 minutes.picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 64.20 minutes.picard.sam.FixMateInformation done. Elapsed time: 58.05 minutes.

然后是GATK对bam文件的一些预处理,步骤是:

RealignerTargetCreator --> IndelRealigner --> BaseRecalibrator --> PrintReads 

后面我会讲到这些步骤是否是必须的。

耗时如下:

INFO  15:50:24,097 ProgressMeter - Total runtime 1165.34 secs, 19.42 min, 0.32 hours INFO  17:21:00,917 ProgressMeter - Total runtime 4265.44 secs, 71.09 min, 1.18 hours INFO  19:58:23,969 ProgressMeter - Total runtime 9436.69 secs, 157.28 min, 2.62 hours INFO  23:41:00,540 ProgressMeter - Total runtime 13349.77 secs, 222.50 min, 3.71 hours 

可以看到对L1这条lane数据来说,整个流程耗时才不到10个小时,还算是可接受范围内的。

接下来用HaplotypeCaller找变异

这个步骤我对realign的bam和recal的bam分别处理了,耗时如下:

INFO  20:40:49,063 ProgressMeter - Total runtime 39243.88 secs, 654.06 min, 10.90 hours INFO  08:53:17,633 ProgressMeter - Total runtime 43939.69 secs, 732.33 min, 12.21 hours 

对bam文件找变异的时间取决于reads数量的多少以及这些reads覆盖的基因组区域大小,虽然一条lane的数据量很小,但它仍然是全基因组测序,如果是全外显子测序这个耗时是不一样的。

整个L1这条lane数据(120M的reads)处理后的文件大小如下所示:

3.0M Jun  7 01:14 L1.bam.bai 8.8G Jun  7 02:33 L1_marked.bam9.0G Jun  7 03:31 L1_marked_fixed.bam3.3M Jun  7 03:36 L1_marked_fixed.bam.bai2.6K Jun  7 02:33 L1.metrics8.2M Jun  5 17:21 L1_realigned.bai9.0G Jun  5 17:21 L1_realigned.bam8.2M Jun  5 23:41 L1_recal.bai 27G Jun  5 23:41 L1_recal.bam8.1M Jun  5 23:53 L1_recal.bam.bai 39M Jun  5 15:50 L1_target.intervals320K Jun  5 19:58 L1_temp.table

因为我的这次基因组重测续数据共5条lane,这5条lane的数据其实是可以并行完成这几个步骤的,最长耗时约12小时。 每个数据处理我都分配了 5个线程, 40G的内存。

merge后不进行AddOrReplaceReadGroups处理

意味着要把5条lane的数据当做是不同的样本,这样后续处理这5个lane的bam矫正以及找变异都是独立进行的,虽然最后生成的vcf文件只有一个,但是每条lane都有独立的基因型。 realign和recal耗时如下:

INFO  17:54:59,396 ProgressMeter - Total runtime 5194.10 secs, 86.57 min, 1.44 hours INFO  02:04:10,907 ProgressMeter - Total runtime 22558.84 secs, 375.98 min, 6.27 hours INFO  18:48:45,762 ProgressMeter - Total runtime 60267.18 secs, 1004.45 min, 16.74 hours INFO  21:30:54,519 ProgressMeter - Total runtime 96119.87 secs, 1602.00 min, 26.70 hours 

同样还是对realign的bam和recal的bam分别用HaplotypeCaller找变异

INFO  19:36:31,960 ProgressMeter - Total runtime 201031.47 secs, 3350.52 min, 55.84 hours INFO  22:59:21,693 ProgressMeter - Total runtime 184944.78 secs, 3082.41 min, 51.37 hours 

可以看到这个时候的耗时相比仅针对一条lane已经是非常恐怖了,近100个小时。不仅仅是因为这个时候reads数量增多,而且是因为5条lane独立样本处理。

merge后需要AddOrReplaceReadGroups处理

这个才是正确的做法,因为不同的lane出来的数据都是我本人的全基因组重测续数据,后续处理应该是当做一个样本的,所有需要AddOrReplaceReadGroups处理,代码是:

### AddOrReplaceReadGroups ###java -Djava.io.tmpdir=$TMPDIR    -Xmx40g -jar $PICARD AddOrReplaceReadGroups \INPUT=${sample}.bam OUTPUT=${sample}_tmp.bam   RGID=jmzeng  RGLB=lib_all  RGPL=illumina RGPU=x10  RGSM=jmzengmv ${sample}_tmp.bam ${sample}.bamsamtools index ${sample}.bam 

这个代码需要结合前面的GATK best practice一起理解。

realign和recal耗时如下:

INFO  15:52:21,739 ProgressMeter - Total runtime 3573.62 secs, 59.56 min, 0.99 hours INFO  22:46:39,615 ProgressMeter - Total runtime 24853.46 secs, 414.22 min, 6.90 hours INFO  16:06:14,340 ProgressMeter - Total runtime 62366.41 secs, 1039.44 min, 17.32 hours INFO  18:18:08,004 ProgressMeter - Total runtime 94304.46 secs, 1571.74 min, 26.20 hours 

这个耗时跟上面把lane当做是独立样本的没有什么区别,因为耗时取决于reads数据量。

接下来找变异,我不仅仅是对realign和recal,还加入了最原始的bam,全部耗时如下:

INFO  02:29:32,680 ProgressMeter - Total runtime 149229.64 secs, 2487.16 min, 41.45 hours INFO  02:15:39,234 ProgressMeter - Total runtime 148379.91 secs, 2473.00 min, 41.22 hours INFO  21:20:51,032 ProgressMeter - Total runtime 130663.06 secs, 2177.72 min, 36.30 hours

可以看到这些耗时都显著的小于把lane当做独立样本。总流程耗时仍然超过了80个小时,这也就是为什么时隔10天我才推出番外~~~

全部流程输出的文件大小如下:

 59G Jun 14 14:17 jmzeng.bam8.4M Jun 14 14:52 jmzeng.bam.bai1.1G Jun 21 02:29 jmzeng_merge_raw.snps.indels.vcf 12M Jun 21 02:29 jmzeng_merge_raw.snps.indels.vcf.idx8.4M Jun 14 22:46 jmzeng_realigned.bai 59G Jun 14 22:46 jmzeng_realigned.bam1.1G Jun 21 02:15 jmzeng_realigned_raw.snps.indels.vcf 12M Jun 21 02:15 jmzeng_realigned_raw.snps.indels.vcf.idx8.5M Jun 16 18:18 jmzeng_recal.bai161G Jun 16 18:18 jmzeng_recal.bam8.5M Jun 16 19:39 jmzeng_recal.bam.bai1.1G Jun 20 21:20 jmzeng_recal_raw.snps.indels.vcf 12M Jun 20 21:20 jmzeng_recal_raw.snps.indels.vcf.idx 47M Jun 14 15:52 jmzeng_target.intervals323K Jun 15 16:06 jmzeng_temp.table

值得注意的是,recal之后的bam文件是原始bam的3倍大小,特别耗费存储空间。

接下来我会讲解realign和recal步骤的必要性,毕竟这两个步骤也的确太耗时了,尤其是recal,不仅仅耗时还特别占硬盘存储。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-06-28

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏Spark学习技巧

Kafka源码系列之使用要点总结及重要错误解决

1,创建一个topic bin/kafka-topics.sh --create --zookeeper localhost:2181 --replicatio...

1936
来自专栏葡萄城控件技术团队

移动APP的自动化测试

 开发移动应用,最耗时耗力的就是手动测试APP的每个功能点或修复bug。有人就会提议App的业务逻辑可以使用nUnit或xUnit测试单元来辅助完成。那用户界面...

1978
来自专栏大魏分享(微信公众号:david-share)

API的计量与限速 | 将一个Web API纳入API管理 |API Management学习第二篇

在本文中,我们将针对:API Management学习第一篇中编写的Restful API,进行纳管。

945
来自专栏非著名程序员

有关使用Universal-Image-Loader的遇到的问题和使用小技巧

? 今天我们来分析一下使用Universal-Image-Loader异步加载图片时遇到的一些问题和解决办法。今天咱们的公众号不分享高大上的原理分析和源码分析...

1808
来自专栏大魏分享(微信公众号:david-share)

后容器时代技术制高点:API管理平台3Scale的架构设计与部署

1832
来自专栏生信技能树

把vcf文件转换为maf格式,肿瘤外显子上游分析教程到此为止

可能还有一些教程我漏掉了,毕竟这些年发布了近万篇教程了,大家直接我去我博客,生信菜鸟团就可以搜索,去我们的论坛,生信技能树里面也可以搜到。

1002
来自专栏杨建荣的学习笔记

SQL审核工具SQL Advisor简单体验

现在的很多大公司,都喜欢招丰富经验的人,从公司的角度来说,能把当前的事务性工作解决了,在这个基础上能够把你的理解和知识沉淀下来,那是极好的,说通俗一些,算是吸...

722
来自专栏NetCore

微信公众平台快速开发框架 For Core 2.0 beta –JCSoft.WX.Core 5.2.0 beta发布

写在前面 最近比较忙,都没有好好维护博客,今天拿个半成品来交代吧。 记不清上次关于微信公众号快速开发框架(简称JCWX)的更新是什么时候了,自从更新到支持.Ne...

2228
来自专栏北京马哥教育

linux存储系统流程简介

存储系统是linux系统非常重要,也是非常基础的知识点。整个存储系统涉及到知识点也非常的多。 本文主要通过磁盘简介->分区管理->文件系统管理->文件存储结构...

3215
来自专栏葡萄城控件技术团队

Visual Studio 2015速递(3)——ASP.NET 新特性

系列文章 Visual Studio 2015速递(1)——C#6.0新特性怎么用 Visual Studio 2015速递(2)——提升效率和质量(VS20...

1816

扫描关注云+社区