关于snp-calling中GATK软件的最佳线程（核心）数

一、前言

对于在服务器上跑测序数据的生信初学者来说，每一个步骤和命令应该请求多少线程总会让人一头雾水。盲目求多不可取，浪费运算资源（和浪费钱~），可是如果线程数太少又不知要跑到猴年马月。本文测试对比了gatk部分命令在不同线程下的运行速度，以供参考。主要是想解决以下三个问题：

1. 请求了32线程，可是运行的命令只调用了8线程
2. 软件占用32线程数，反而运行速度比16线程和8线程更慢
3. 软件在32线程下比16线程下只快了1.5倍，应该该选哪一个？

（节约时间的话直接拉到最后看结果）

二、gatk命令测试

1. gatk HaplotypeCaller

印象里做snp-Calling的时候比较费时间的就是这一步了，可以从官网查阅得知，HaplotypeCaller的默认调用的线程数就是4，所以如果我们提交任务的时候不额外指定，那么不管找服务器要几个线程，它都只调用4个，运行如下命令。

gatk  HaplotypeCaller -R XX -I XX --emit-ref-confidence GVCF --min-base-quality-score 10 -O z --tmp-dir vcf/

60分钟的结果如下：

20:05:27.799 INFO  IntelPairHmm - Available threads: 4
20:05:27.800 INFO  IntelPairHmm - Requested threads: 4
20:05:27.800 INFO  PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
20:05:27.860 INFO  ProgressMeter - Starting traversal
20:05:27.860 INFO  ProgressMeter -        Current Locus  Elapsed Minutes     Regions Processed   Regions/Minute
21:03:52.075 INFO  ProgressMeter -           1:31721093             58.4                237390           4064.6
21:04:02.379 INFO  ProgressMeter -           1:31794779             58.6                237970           4062.6
21:04:12.623 INFO  ProgressMeter -           1:31865291             58.7                238530           4060.4
21:04:22.939 INFO  ProgressMeter -           1:31948694             58.9                239200           4059.9
21:04:32.942 INFO  ProgressMeter -           1:32028143             59.1                239850           4059.4
21:04:43.333 INFO  ProgressMeter -           1:32144224             59.3                240760           4062.9
21:04:53.965 INFO  ProgressMeter -           1:32202234             59.4                241210           4058.4
21:05:04.326 INFO  ProgressMeter -           1:32291309             59.6                241900           4058.2
21:05:14.379 INFO  ProgressMeter -           1:32337518             59.8                242290           4053.3
21:05:24.386 INFO  ProgressMeter -           1:32398011             59.9                242790           4050.4
21:05:34.470 INFO  ProgressMeter -           1:32461380             60.1                243310           4047.7

下面这种情况是提交的任务请求了8核心，但是软件只用了4核心，我们可以看到gatk运行前的信息报告如下Available threads: 8，Requested threads: 4，我们一定要避免这种情况。

15:38:40.151 INFO  IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
15:38:40.152 INFO  IntelPairHmm - Available threads: 8
15:38:40.153 INFO  IntelPairHmm - Requested threads: 4
15:38:40.153 INFO  PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
15:38:40.199 INFO  ProgressMeter - Starting traversal
15:38:40.200 INFO  ProgressMeter -        Current Locus  Elapsed Minutes     Regions Processed   Regions/Minute

我们加一句–native-pair-hmm-threads 8让gatk调用8个线程试试

gatk  HaplotypeCaller --native-pair-hmm-threads 8 -R XX -I XX --emit-ref-confidence GVCF --min-base-quality-score 10 -O z --tmp-dir vcf/

60分钟结果如下：

14:44:05.714 INFO  IntelPairHmm - Available threads: 8
14:44:05.714 INFO  IntelPairHmm - Requested threads: 8 #从这里我们可以看出确实用上了8个线程
14:44:05.714 INFO  PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
14:44:05.762 INFO  ProgressMeter - Starting traversal
14:44:05.762 INFO  ProgressMeter -        Current Locus  Elapsed Minutes     Regions Processed   Regions/Minute
15:42:30.129 INFO  ProgressMeter -           1:29557912             58.4                221020           3784.2
15:42:40.504 INFO  ProgressMeter -           1:29682598             58.6                221980           3789.4
15:42:50.932 INFO  ProgressMeter -           1:29786339             58.8                222710           3790.6
15:43:00.934 INFO  ProgressMeter -           1:29933997             58.9                223770           3797.9
15:43:10.987 INFO  ProgressMeter -           1:30044884             59.1                224570           3800.7
15:43:21.112 INFO  ProgressMeter -           1:30166428             59.3                225530           3806.0
15:43:31.294 INFO  ProgressMeter -           1:30245351             59.4                226170           3805.9
15:43:41.962 INFO  ProgressMeter -           1:30338703             59.6                226900           3806.8
15:43:52.111 INFO  ProgressMeter -           1:30388809             59.8                227350           3803.6
15:44:02.181 INFO  ProgressMeter -           1:30443667             59.9                227800           3800.4
15:44:12.181 INFO  ProgressMeter -           1:30538540             60.1                228460           3800.9

可以明显看出来，60分钟跑下来，4线程的运行速度一直比8线程略快，最后前者243310个regions，后者228460个regions，4线程下运行速度提升了6.5％左右。

所以在这里提升线程数对速度也没有帮助，建议小伙伴们提交任务的时候不要请求太多线程喔~。尽量是低线程多任务运行程序

2.MarkDuplicates

标记重复这一步，手册上说如果不额外指定，那么就会调用所有可用核心。（额外指定可以用–executor-cores参数）。值得一提的是，服务器上通常会关闭核心的超线程功能，所以一个核心在这里就是对应一个线程。

这里我查询了一下，gatk该模块是调用的picard，建议跟samtools sort软件保持一致，用4核心就好啦~。

3.MergeVcfs

接下来就不一一举例了，mergevcf我只测试了8,16,32,推荐的线程数是8，增加线程提升较小，印象中使用8核心到16核心速度可以提升1.2-1.4倍，但是32完全无提升。这里就是我们开篇说的第三种情况了，增加了大量核数却没有提升多少速度。(感兴趣的自行测试)

ps：这个最佳线程数的坑我一年前就开了想写，但是由于各种原因没有完成，现在有空闲时间就拿起来写了。但是，部分当年的测试结果已经找不到了，也没有十几个T的空间和资源给我跑calling了（hhh）。

三、官方资料

经过我对耗时最长的部分GATK calling命令测试，得出的结论是跑GATK的全流程中，4-8核心最佳，有的命令“Requested threads: 8”，甚至跑得比4核心更慢了（点名HaplotypeCaller）。

因为不想继续测试，就开始找GATK的官方解释，许多命令没有明确讲多少线程合适，但是找到了官方团队对于多线程的一个评论回复

“Instead of the 24 cores in local mode, we would recommend multiple executors with 8 cores each. ”
GATK4.0开始放弃了自己实现多线程任务，选用了现成的SPARK系统（放弃重复造轮子）。大概意思就是GATK调用SPARK进行多线程的时候超过8核心的话效率会比较低，和我们的测试结果一致。

四、总结

下面贴一下我个人的推荐吧，欢迎指出问题，如果有人做更详尽的测试也可以贴评论区告诉我。

命令	核心数目
bwa	4-8
gatk SortSam（samtools sort）	4
gatk MarkDuplicates	4-8
gatk HaplotypeCaller	4
gatk GenotypeGVCFs	4-8
gatk MergeVcfs	4-8
所有压缩和解压任务、索引建立的任务	1

选用合适的核心数，有时可以避免几倍的资源浪费，甚至还可以节约时间，有人看的话我再去开坑把其他常用软件都测试一下（逃）

五、参考

https://gatk.broadinstitute.org/hc/en-us/community/posts/4445958815003-MarkDuplicatesSpark-speed-estimates

https://github.com/broadinstitute/gatk#running-gatk4-spark-tools-on-a-spark-cluster