引言 生命科学方兴未艾， 从食品工业中的细菌培养鉴定到癌症快速诊断，基于 DNA 分析的应用不断出现，但同时基因分析应用也面临着很大挑战；许多新技术、新方法被应用到基因序列分析应用中，包括 Spark、FPGA 以及 GPU 协处理器加速等，这些技术的应用不仅能够使大部分生命科学领域的应用，包括开源和 ISV 软件，在不需要复杂的 MPI 编程情况下实现并行化处理，同时 Spark 内存内计算技术也能够提高分析效率，加速工作流程， 缩短分析时间，从而有更多新的发现。本文将介绍如何利用 Spark 技术运行常用的基因序列分析应用，包括在 Spark 不同模式下的运行方法， 运行过程以及运行结果分析，并比较在不同运算平台以及不同运行参数情况下的性能和加速比。
1. 基因序列分析工作流 基因序列分析工作流以 GATK 的最佳实践为标准。它以最初的 FASTQ 文件为输入，从 BWA-mem 测序到 GATK 的 HaploTyperCaller，完成对整个样板的测序分析。
 图 1、GATK 最佳实践
在测序工作流的第一阶段，BWA-mem 对输入文件 FASTQ 执行比对，生成序列比对和映射文件 SAM，然后通过 SortSam 生成一个经过排序的 BAM 文件，实际上，BAM 文件是 SAM 文件的二进制形式，此后的处理均基于 BAM 二进制文件。
BAM 文件传送给 Picard 工具 MarkDuplicates, 去除重复的片段，并生成一个合并的、去除重复片段的 BAM 文件。以下的几步，RealignerTargetCreator、IndelRealigner、BaseRecalibrator、PrintReads 和 HaplotypeCaller 都是 GATK 的一部分，是对高吞吐序列数据进行分析的软件包。
 图 2、基因测序工作流分解
序列分析的主要工作是数据前处理，经过处理的数据可以为后续的分析工作所调用。前处理阶段，比对和排序是计算密集且比较耗费时间的过程，尽管通过多处理器或多线程的方式可以提高效率，但是在实际工作中由于计算方法的复杂程度以及需要分析的数据量迅速增加，当前一个分析过程仍然可能花费超过 1 天时间，费用从 200 美元到 600 美元不等。Spark 技术可以将串行的分析并行化，将数据分段优化并进行动态负载均衡以提高效率。
GATK4 就是 Broad 推出的基于 Spark 技术的基因序列分析软件包。GATK4 数据前处理的流程是：
 图 3、基于 Spark 技术数据前处理
- 合并输入文件和参考文件
- 分成数据块
- 数据块的数量取决于集群大小和可用资源
- 工作流分为多级，在数据处理之前只划分一次
- 分级类似于 Mapreduce
- TaskManager 分配任务给 executor
- BlockManager 利用 spark.broadcast.blockSize 设置块的每一个片的大小；数值太大，在广播过程中会减小并发（使运行变慢），但是，数值太小时，BlockManager 的性能会受影响，缺省是 4M
- 剩余内存空间不断减少，当剩余内存太小时，任务会中断并被踢出。Spark 会尝试重启任务，当超过设定的重启次数仍无法成功时，作业就非正常结束。
GATK4 仍处在不断开发、不断完善的过程中，其所提供的工具和工作流也在不断增加，目前最新版本的 GATK4 提供的工作流包括：
- BQSRPipelineSpark 在 Spark 上执行 BQSR 的两个步骤(BaseRecalibrator 和 ApplyBQSR)
- BwaAndMarkDuplicatesPipelineSpark 以名称排序的文件为输入运行 BWA 和 MarkDuplicates.
- ReadsPipelineSpark 以 BWA 的校正片段为输入运行 MarkDuplicates 和 BQSR，其输出用于后续的分析
图 4 是 ReadsPipelineSpark 工作流示意图。
 图 4、 ReadsPipelineSpark 工作流
2. 序列分析优化方法 序列分析中的不同应用对系统资源有同步的需求。从图 5 可以看到，有的应用占用 CPU 资源比较高，如 BWA 不仅占用大量处理器资源，且运行时间长，而有的应用则需要大量内存，处理时间同样比较长，如 HaploTyperCaller。
 图 5、不同应用对系统资源的需求
一般地，有四种方法对基因处理过程进行优化和加速：
- -nt 在引擎 engine 级别进行并行，并行处理基因序列的不同部分
- -nct 在 walker 级别进行并行，加速处理基因序列每个单独区域
- MapReduce 同时生成许多实例，每个实例处理基因序列的不同的（任意的）部分
- 利用科学库优化
在 GATK 工作流中可以通过设置-nt 和-nct 参数，提高作业运行效率。
 GATK4 是 GATK 基于 Spark 开发的版本，它有很多可以在 Spark 环境中运行的工具和工作流，它采用分级的方式运行作业，其工作工程类似于 Mapreduce 。有 3 种运行模式：
- None-spark standalone 模式
- Spark standalone 模式
- Spark cluster 模式
在输入数据和参考文件设置正确的情况下，大部分 GATK4 工具都可以在 Spark 集群模式下成功运行。 有些应用在集群模式下的运行结果可以得到显著的提升，如 CountReadsSpark，有的应用，特别是当工作流则需要更多的系统资源时，在 spark standalone 模式下无法运行，会报告"Not enough space to cache RDD in memory"错误，而在 Spark Cluster 模式下则能顺利运行， 如 CollectInsertSizeMetricsSpark 。
 图 6、CollectInsertSizeMetricsSpark 结果
对序列分析应用加速的另一种方法是提供基于 POWER8 处理器的优化科学库。以 HaplotypeCaller 分析为例，它在分析过程中占用大量的内存，运行时间最长。不同厂家也在开发基于自己软件栈的加速库，如 Intel GKL 基因内核库。
 IBM 提供一个在 POWER8 系统上优化的 PairHMM 算法，它充分发挥了 POWER8 系统上新的软件、硬件特性，目前，该优化科学库可以运行在 POWER8 Ubuntu14 和 RHEL7 操作系统上。
最新版本的科学库利用与 POWER8 上 Java 相同的浮点精度对 HaplotyperCaller 进行加速，同时它充分利用并发多线程 SMT 以及向量指令集，对 HaplotypeCaller 的加速性能超过以前的版本，特别是在单线程模式（即 –nct 选项未指定）。在单线程模式下，利用 PairHMM 加速，HaplotypeCaller 消耗的时间只有一半，加速比达到 1.88 倍。
在 P8 系统上调用 PairHMM：
 加速基因测序应用的其他方法还包括 FPGA、GPGPU 计算以及全硬件加速的 Edico Dragon Solution 等，不在本文档讨论范围。 |
