善学者,假人之长以补其短

本文总结了生物信息人员在使用BLASTX的一些经验,期待与其他同行交流BLASTX的优化经验。

0 前言

BLAST作为生物信息最重要的局部比对工具,在序列物种注释和基因注释起重要作用。其中BLASTX(将query序列翻译成氨基酸序列和蛋白数据库进行比对)常常作为对核酸序列的基因注释工具。BLASTX先要将DNA序列按照6个读码框翻译成氨基酸序列,从而BLASTX的比对时间和需要计算资源都是比较大的,如何减少运行时间和降低运算成本,成为生信工作人员关心的重要问题之一。另外,本文主要讨论BLASTX工具的优化,其他BLASTX的替代工具,比如DIAMOND等,不在此文的讨论范围。

1 使用ncbi-blast+工具替换blastall

NCBI在1989年引入BLAST软件包,现在叫legacy_blast,主要是通过blastall 去调用 blastx,blastall最后一个版本更新在2016年。2009年,NCBI发布了检索速度更快,输出格式更为多样化的ncbi-blast+软件包。BLAST+的特点为:擅长处理长查询序列,将长序列切成短序列,即减少了内存消耗又更好的利用CPU计算架构,提高了比对效率,详情见Command Line Applications User Manual。BLAST+软件包将BLASTN,BLASTP和BLASTX单独拆分出来并分别进行了优化和处理,因而计算速度比blastall有了很大的提升。具体测试结果可以看ncbi-blast+的user_manual

2 使用blastx-fast模式

ncbi-blast+软件包对blastx-fast的比对模式介绍比较少,没有详细说明fast模式用了什么参数或者算法,本人的猜想可能是做了贪婪算法,增加或者减少一些参数阈值,比如增大word-size参数值等。blast2go是作为核酸比对蛋白序列使用广泛的工具包,在进行序列比对时,采用的BLASTX参数是: -task blastx-fast。本文作者也是用一些数据进行了测试。ncbi-blast+版本: ncbi-blast+ 2.4.0。query数据:MEGA-HIT软件组装双端测序数据,得到的contigs核酸序列和截取的其中10,000条contigs的部分核酸序列(final.contigs.fa和final.contigs.10000.fa),数据详情见下表1

表1 query数据(DNA)
sample_name num_seqs sum_len(bp) min_len(bp) avg_len(bp) max_len(bp)
final.contigs.fa 142,983 188,296,603 200 1,316.9 390,130
final.contigs.100000.fa 10,000 14,499,134 200 1,449.9 388,898

database:从NR库截取100,000条蛋白序列(NR_100000.fa),数据库详情见下表2

表2 database 数据(Protein)
sample_name num_seqs sum_len(bp) min_len(bp) avg_len(bp) max_len(bp)
NR_100000 100,000 32,073,119 6 320.7 34,350

计算所用机器类型:阿里云的弹性计算-普通实例(8CPU和16GB内存)。

为了比较blastx普通模式和blastx-fast模式,运行参数为表3

表3 运行参数
任务类型 参数
blast -num_threads 8 -max_target_seqs 5 -evalue 0.00001 -outfmt 5 -task blastx
blastx-fast -num_threads 8 -max_target_seqs 5 -evalue 0.00001 -outfmt 5 -task blastx-fast

数据的运行时间,如下表4

表4 blastx和blastx-fast 运行时间
software seq_name db_name run_time(min)
blast final.contigs.10000.fa NR_100000 277
blast-fast final.contigs.10000.fa NR_100000 90
blast final.contigs.fa NR_100000 3,621
blastx-fast final.contigs.fa NR_100000 1,105

从运行时间来看,blastx-fast 大约为BLASTX的 1/3。本文采用的是基因组组装结果序列,contigs的长度最长有390多kb,BLAST+软件包对长的查询文件优化提升很多。根据多次的经验来看,BLASTX-fast的运行时间大约是BLASTX的1/4~1/3。一致性比较如下表5

表5 结果一致性比较表
seq_name blastx_query_num blastx-fast_query_num 交集数 (交集数/blastx_query)% (交集数/blastx-fast_query) blastx_特有 blastx-fast_特有
final.contigs.10000.fa 5,836 5,541 5,541 94.945 100.000 295 0
final.contigs.fa 81,262 77,175 77,154 94.944 99.92 295 21

交集/blastx-fast_query是为了比较公共集合占总query数的百分比。从运行结果来看,发现 blastx-fast模式的结果,基本上属于BLASTX结果的子集。在14万条查询序列也能保持很高的一致性。本文推测,blastx-fast比blastx的结果要少一些,漏掉的比例在 5%左右,漏掉的结果可能是假阳性或者真阳性结果。

blast2go在其官网一些案例使用的是blastx-fast模式。官方帮助文档,参考第10页也是使用blastx-fast。

3 充分利用cpu并行任务

对于通过云计算申请的机子,计费按使用时间算,尽量最大限度地使用计算资源,一方面减少运行时间,另一方面减少计算花费。当一条query在和database比对的时候,可能只占了部分内存和部分CPU,如果能在申请的实例,在保证多个任务不会带来太大的IO和线程切换问题,可以尝试尽量多并行几个任务。观察单个任务的CPU和内存的使用状态,发现虽然BLASTX设置的CPU为8个,但是实际使用的却只有200%,空闲600%的CPU。根据上述的观察结果,设计了测试,如表6所示

表6 不同数据切分和cpu的组合试验
切份 seq_name db_name threads(n) *Mem(MB) *time(min)
*1 final.contigs.10000.fa NR_100000 8 6358 90
2 final.contigs.10000.fa NR_100000 4 4290 43
4 final.contigs.10000.fa NR_100000 2 5216 26
8 final.contigs.10000.fa NR_100000 1 7881 18
1 final.contigs.fa NR_100000 8 12751 1105
2 final.contigs.fa NR_100000 4 13847 601
4 final.contigs.fa NR_100000 2 14165 323
8 final.contigs.fa NR_100000 1 *- -

注:

  • “-”表示没有数据,是因为内存超出16GB,任务失败。
  • 切分为1的query指没有进行切分。
  • Mem,运行内存最大值,单位MB。
  • *time: 切分数大于1,总任务时间按照实际运行时间最长的任务算。单位:分钟。

为了确保总任务的总CPU数一致,设置 切分数*threads=8。随着切分的增加,该机器实际消耗的内存也是显著增加。随着切分的增加,该机器运行的总时间也在显著减少。由于BLASTX的任务时间运行,运行内存与query数据和database数据有关。对于长序列比较多的query序列,可能IO或者CPU使用率比较低,多增加并行线程对计算效率有很高的提升空间。对于短序列比较多的query序列,IO和CPU使用率高,多切分可能不会有提升计算效率。在进行安排切分策略前需要对query和database做个简单模拟,看看不切分的任务是否已计算已经饱和。

4 总结

本文只是根据使用BLASTX的经验进行了若干测试,由于测试任务比较少,以上结论可能并不准确,只是向生信人员及其他对BLASTX感兴趣的人员抛出一个优化方案。如需要深入优化BLASTX等软件,需要从blast的源代码入手,结合对CPU运算研究、计算过程中的IO调研并结合任务的预期成本,选择合适的优化策略。

5 致谢

感谢张求学和武雅蓉对本测试提供了优化思路,感谢聚道其他同学帮忙审阅了博文并提供了修改意见。