开源软件的性能
生物信息学中的大数据处理

引言

生物信息学和大数据

生物信息学领域旨在提供工具和分析，以促进对地球生命分子机制的理解，主要通过分析和关联基因组和蛋白质组信息。随着越来越多的基因组信息（包括基因组序列和表达基因序列）变得可用，更高效、更灵敏和更具体的分析变得至关重要。

在DNA测序中，化学和机械过程基本上将DNA和RNA中存在的信息“数字化”。这些序列使用每个核苷酸一个字母的字母表记录。对这些序列数据进行各种分析，以确定其如何构建成更大的构建块以及它如何与其他序列数据相关联。这为生物进化和发育、遗传学以及日益增多的医学研究奠定了基础。

核苷酸链数据来自测序过程，以称为reads的字母串的形式出现。（在生物信息学意义上使用read一词与在计算机科学和软件工程意义上使用该词发生了不幸的冲突。尤其是在读取reads的性能可以调整的情况下，这一点尤其如此，正如我们即将讨论的那样。为了消除这种不幸的冲突，我们指的是来自基因组的序列为基因组reads。）为了分析更大规模的结构和过程，必须将多个基因组reads组合在一起。这种拟合不同于拼图，因为图片通常事先未知，并且碎片可能（并且经常）重叠。另一个复杂之处在于，并非所有基因组reads都具有完美的保真度，并且可能包含各种错误，例如插入或删除字母或用错误的字母替换核苷酸。虽然拥有冗余reads可以帮助组装或拟合拼图碎片，但由于所有现有测序技术中都存在这种不完美的保真度，它也是一个障碍。错误的基因组reads的出现随着数据量的增加而增加，这使得数据的组装变得复杂。

随着测序技术的改进，产生的序列数据量已开始超过采用传统方法分析此类数据的计算机硬件的能力。（测序技术中的许多最先进技术会产生大量的基因组reads，通常为数千万到数十亿个，每个reads的序列为50到100个核苷酸。）预计这种趋势将持续下去，并且是高性能计算(HPC)、分析和信息科学界所知的大数据[1]问题的一部分。随着硬件成为限制因素，人们越来越关注通过软件解决方案来缓解该问题的方法。在本章中，我们介绍了一种这样的软件解决方案，以及我们如何对其进行调整和扩展以处理TB级数据。

我们的研究重点是高效的预处理，其中各种过滤器和分箱方法修剪、丢弃和分箱基因组reads，以改进下游分析。这种方法的好处是限制了需要对下游分析进行的更改，下游分析通常直接使用基因组reads。

在本章中，我们将介绍我们的软件解决方案，并描述我们如何对其进行调整和扩展以有效处理越来越大的数据量。

什么是khmer软件？

Khmer是我们的一套软件工具，用于在使用传统生物信息学工具[2]进行分析之前预处理大量基因组序列数据——与东南亚土著民族没有任何关系。这个名字是通过与k-mer一词的自由联想得来的：作为预处理的一部分，遗传序列被分解成给定长度k的重叠子串。由于许多分子的链通常称为聚合物，因此特定数量分子的链称为k-mer，每个子串代表一条这样的链。请注意，对于每个基因组reads，k-mer的数量将是序列中核苷酸的数量减去k加1。因此，几乎每个基因组reads都将分解成许多重叠的k-mer。

图12.1 - 将基因组序列分解成4-mer。在khmer中，每个k-mer的前向序列和反向互补序列都被散列到相同的值，以识别DNA是双链的。参见未来方向。

由于我们想告诉您我们如何测量和调整这部分开源软件，因此我们将跳过其背后的许多理论。需要说明的是，k-mer计数是其大部分操作的核心。为了紧凑地计算大量k-mer，使用了称为布隆过滤器[3]的数据结构（）。有了k-mer计数，我们就可以将高度冗余的数据排除在进一步处理之外，此过程称为“数字归一化”。我们还可以将低丰度序列数据视为可能的错误，并将其排除在进一步处理之外，这是一种错误修剪方法。这些归一化和修剪过程大大减少了下游分析所需的原始序列数据量，同时在很大程度上保留了感兴趣的信息。

图12.2 - khmer软件的分层视图

软件的核心是用C++编写的。核心包括数据泵（将数据从在线存储移动到物理RAM的组件）、几种常见格式的基因组reads解析器和几个k-mer计数器。一个应用程序编程接口(API)构建在核心之上。这个API当然可以从C++程序中使用，就像我们对一些测试驱动程序所做的那样，但也作为Python包装器的基础。Python包构建在Python包装器之上。许多Python脚本与包一起分发。因此，khmer软件作为一个整体，是核心组件（用C++编写以提高速度）、高级接口（通过Python公开以方便操作）以及各种工具脚本（提供执行各种生物信息学任务的便捷方法）的组合。

khmer软件支持多个阶段的批处理操作，每个阶段都有单独的数据输入和输出。例如，它可以获取一组基因组reads，计算其中的k-mer，然后可以选择保存布隆过滤器哈希表以供以后使用。稍后，它可以使用保存的哈希表对新的一组基因组reads执行k-mer丰度过滤，并保存过滤后的数据。这种重用早期输出和决定保留什么的灵活性允许用户根据自己的需求和存储约束定制特定过程。

图12.3 - 数据流经khmer软件

大量数据（可能是TB级）必须由软件从磁盘移动到内存。拥有高效的数据泵至关重要，因为从存储到CPU的输入吞吐量可能比从物理RAM到CPU的数据传输吞吐量低三个甚至四个数量级。对于某些类型的数据文件，必须使用解压缩器。在任何一种情况下，解析器都必须有效地处理结果数据。解析任务围绕可变长度行展开，但还必须考虑无效的基因组reads并保留某些生物信息，这些信息可以在后续组装期间利用，例如序列片段末端的配对。每个基因组reads都被分解成一组重叠的k-mer，并且每个k-mer都与布隆过滤器注册或与之比较。如果正在更新或用于比较先前存储的布隆过滤器，则必须将其从存储中加载。如果正在为以后使用或更新创建布隆过滤器，则必须将其保存到存储中。

数据泵始终对文件执行顺序访问，并且可能会被要求一次读取大量数据。考虑到这一点，以下是一些想到的问题

• 我们是否充分利用了数据按顺序访问的事实？

• 是否预取了足够多的数据页到内存中以最大程度地减少访问延迟？

• 是否可以使用异步输入而不是同步输入？

• 我们能否有效地绕过系统缓存以减少内存中的缓冲区到缓冲区的复制？

• 数据泵是否以不会产生任何不必要的访问器或决策逻辑开销的方式向解析器公开数据？

解析器的效率至关重要，因为数据以相当宽松的字符串格式输入，并且必须在进行任何进一步处理之前转换为内部表示形式。由于每个单独的数据记录都相对较小（100-200字节），但记录数却有数百万到数十亿个，因此我们投入了相当多的精力来优化记录解析器。解析器在其核心是一个循环，它将数据流分解成基因组reads并将它们存储在记录中，在此过程中执行一些初始验证。

关于解析器效率的一些考虑因素是

• 我们是否最大程度地减少了解析器访问内存中数据的次数？

• 在从数据流中解析基因组reads时，我们是否最大程度地减少了缓冲区到缓冲区的复制次数？

• 我们是否最大程度地减少了解析循环内的函数调用开销？

• 解析器必须处理杂乱的数据，包括模糊的碱基、过短的基因组reads和字符大小写。这种DNA序列验证是否以尽可能高效的方式完成？

为了迭代基因组reads中的k-mer并对其进行散列，我们可以询问

• k-mer迭代机制是否可以针对内存和速度进行优化？

• 布隆过滤器散列函数是否可以以任何方式进行优化？

• 我们是否最大程度地减少了散列器访问内存中数据的次数？

• 我们能否批量递增散列计数以利用缓存？

分析和测量

仅仅阅读源代码并关注性能就揭示了许多改进的领域。但是，我们希望系统地量化在代码各个部分花费的时间。为此，我们使用了几个分析器：GNU分析器（gprof）和调整和分析实用程序(TAU)。我们还在源代码本身中创建了仪器，允许对关键性能指标进行细粒度的查看。

代码审查

盲目地应用工具来测量系统（软件或其他）很少是一个好主意。相反，在测量系统之前，通常最好先了解一下系统。为此，我们首先通过肉眼审查了代码。

手动追踪不熟悉代码的执行路径是一个好主意。（其中一位作者 Eric McDonald 在加入项目时不熟悉 khmer 软件，并且他就是这么做的。）虽然确实有分析器（和其他工具）可以生成调用图，但这些图只是抽象的总结。实际遍历代码路径并查看函数调用是一种更具沉浸感和启发性的体验。调试器可用于此类遍历，但并不容易用于探索鲜有问津的代码路径。此外，一步一步地遍历执行路径可能会非常乏味。断点可用于测试代码中的某些点在正常执行期间是否会被命中，但设置断点需要对代码有一些先验知识。作为替代方案，使用具有多个窗格的编辑器效果很好。四个显示窗格通常可以同时捕获一个人在任何给定时刻需要知道的信息——以及能够处理的信息。

代码审查显示了一些问题，其中一些（但不是全部）后来得到了分析工具的证实。我们注意到的一些问题包括：

• 我们预计 k-mer 计数逻辑的流量最高。

• 代码流量最高区域存在对 toupper 函数的冗余调用。

• 基因组读取的输入是逐行进行的，并且按需进行，没有任何预读调整。

• 对每个解析的有效基因组读取结构体进行了按值复制。

尽管上述内容看起来像是相当强烈的自我批评，但我们想强调的是，到目前为止，khmer 更侧重于实用性和正确性。我们的目标是优化现有的、大部分正确的软件，而不是从头重新开发它。

工具

分析工具主要关注在代码的任何特定部分花费的时间量。为了测量这个量，它们在编译时将检测代码注入到代码中。这种检测确实会改变函数的大小，这可能会影响优化期间的内联。检测还会直接给总执行时间带来一些开销；特别是，对代码高流量区域的分析可能会导致相当大的开销。因此，如果您还测量代码的总执行时间，则需要注意分析本身如何影响它。为了衡量这一点，可以使用一个简单的外部数据收集机制，例如 /usr/bin/time，来比较相同的一组优化标志和操作参数下的非分析和分析执行时间。

我们通过测量一系列 k 大小下分析和未分析代码之间的差异来衡量分析的影响——较小的 k 值会导致每个基因组读取更多的 k-mer，从而增加分析器特定的影响。对于 k = 20，我们发现未分析代码的运行速度比分析代码快约 19%，对于 k = 30，未分析代码的运行速度比分析代码快约 14%。

在任何性能调整之前，我们的分析数据显示，正如我们通过观察预测的那样，k-mer 计数逻辑是代码中流量最高的区域。稍微令人惊讶的是，它所占的比例有多大（大约占总时间的 83%），与针对存储的 I/O 操作形成对比（对于一个特定的介质和低带宽争用，大约占总时间的 5%）。鉴于我们的试验数据集约为 500 MB 和 5 GB，我们没有预料到会看到很多缓存效应。¹ 事实上，当我们控制缓存效应时，我们发现它们最多不超过几秒钟，因此与总执行时间的误差范围相比，并不大。这让我们意识到，在代码优化过程的这个阶段，I/O 并不是我们的主要瓶颈。

一旦我们开始对 khmer 软件进行并行化，我们就编写了一些驱动程序，这些程序使用 OpenMP [4] 来测试我们对各种组件的并行化。虽然 gprof 擅长分析单线程执行，但它缺乏在使用多个线程时跟踪每个线程执行的能力，并且它不理解并行化机制，例如 OpenMP。对于 C/C++ 代码，OpenMP 并行化由编译器编译指示确定。GNU C/C++ 编译器在 4.x 系列版本中，如果提供 -fopenmp 开关，则会遵守这些编译指示。当遵守 OpenMP 编译指示时，编译器会在编译指示的位置以及与它们关联的基本块或其他分组周围注入线程处理检测。

由于 gprof 无法轻松地提供我们所需的每个线程报告和 OpenMP 支持，因此我们转向了另一个工具。这是俄勒冈大学领导的合作项目开发的 Tuning and Analysis Utilities (TAU) [5]。市面上有很多并行分析工具——其中许多专注于使用 MPI（消息传递接口）库的程序，这些库在某些类型的科学计算任务中很流行。TAU 也支持 MPI 分析，但由于 MPI 对于 khmer 软件在当前形式中并不是一个真正的选择，因此我们忽略了 TAU 的这方面。同样，TAU 也不是唯一可用于每个线程分析的工具。每个线程分析和能够与 OpenMP 紧密集成相结合是我们选择它的原因之一。TAU 也是完全开源的，并且不依赖于任何一家供应商。

虽然 gprof 仅依赖于在编译时注入到源代码中的检测（以及链接的一些额外部分），但 TAU 也提供了此检测和其他检测选项。这些选项是库插桩（/主要用于 MPI 分析）和二进制文件的动态插桩。为了支持这些其他选项，TAU 提供了一个执行包装器，称为 tau_exec。源代码的编译时检测通过一个包装脚本（称为 tau_cxx.sh）来支持。

TAU 需要额外的配置才能支持某些分析活动。例如，要获得紧密的 OpenMP 集成，TAU 需要配置并构建为支持 OPARI。类似地，要使用较新 Linux 内核公开的性能计数器，它需要配置并构建为支持 PAPI。此外，一旦构建了 TAU，您可能希望将其集成到您的构建系统中以方便使用。例如，我们设置了构建系统，以便在需要 TAU 分析时，可以使用 tau_cxx.sh 包装脚本作为 C++ 编译器。如果您尝试构建和使用 TAU，您肯定需要阅读文档。虽然比 gprof 功能强大得多，但它远没有那么容易上手或直观。

手动检测

使用独立的外部分析器检查软件的性能是快速便捷地了解软件各个部分执行时间的一种方法。但是，分析器通常不擅长报告代码在特定函数内的特定自旋锁中花费了多少时间，或者您的数据的输入速率是多少。为了增强或补充外部分析功能，可能需要手动检测。此外，手动检测可能比自动检测侵入性更小，因为您可以直接控制观察的内容。为此，我们创建了一个可扩展的框架，用于在软件本身内部测量吞吐量、迭代次数以及原子或细粒度操作周围的时间等内容。为了确保我们自己诚实，我们在内部收集了一些可以与外部分析器测量的结果进行比较的数字。

对于代码的不同部分，我们需要有不同的指标集。但是，所有不同的指标集都有一些共同点。其中一点是它们大多是时间数据，并且您通常希望在执行期间累积时间数据。另一点是，需要一致的报告机制。考虑到这些因素，我们为所有不同的指标集提供了一个抽象基类 IPerformanceMetrics。IPerformanceMetrics 类提供了一些便捷方法：start_timers、stop_timers 和 timespec_diff_in_nsecs。启动和停止计时器的方法分别测量经过的实际时间和经过的每个线程 CPU 时间。第三种方法计算两个标准 C 库 timespec 对象之间的纳秒级差异，这对于我们的目的来说已经足够精确了。

为了确保手动插入的内部检测的开销在生产代码中不存在，我们小心地将其包装在条件编译指令中，以便构建可以指定将其排除在外。

调整

使软件更有效地工作是一件非常令人欣慰的事情，尤其是在面对数万亿字节数据通过它时。我们的叙述现在将转向我们为提高效率而采取的各种措施。我们将这些措施分为两部分：优化输入数据的读取和解析以及优化布隆过滤器内容的操作。

通用调整

在深入探讨我们如何调整 khmer 软件的一些细节之前，我们想简要介绍一些通用性能调整选项。生产代码通常使用一组安全且简单的优化来构建；这些优化通常可以证明不会改变代码的语义（即引入错误），并且只需要进行一次编译。但是，编译器确实提供了其他优化选项。这些其他选项可以广泛地归类为激进优化（这是编译器文献中一个相当标准的术语）和配置文件引导优化(PGO) [6]。（严格来说，这两个类别并不相互排斥，但通常涉及不同的方法。）

激进优化在某些情况下可能不安全（即引入错误），或者在其他情况下实际上会降低性能。激进优化可能由于各种原因而不安全，包括浮点精度方面的粗心或关于不同操作数与不同内存地址相关联的假设。激进优化也可能特定于特定的 CPU 系列。配置文件引导优化依赖于分析信息，以便在编译和链接期间对如何优化程序做出更明智的猜测。一个经常看到的配置文件引导优化是位置优化——尝试将高度相关的函数作为可执行映像文本段中的邻居放置在一起，以便它们在运行时一起加载到相同的内存页面中。

在项目的这个阶段，我们避免了这两类额外的优化，而是选择了有针对性的算法改进——这些改进可以在许多不同的 CPU 架构上提供好处。此外，从构建系统复杂性的角度来看，激进优化可能会产生可移植性问题，而配置文件引导优化会增加可能失败的活动部件的总数。鉴于我们不为各种架构分发预编译的可执行文件，并且我们的目标受众通常不太了解软件开发或构建系统的复杂性，我们很可能会继续避免这些优化，直到我们认为其好处大于缺点。鉴于这些考虑，我们的主要重点是提高算法的效率，而不是其他类型的调整。

数据泵和解析器操作

我们的测量结果表明，计算 k-mer 所花费的时间超过了从存储中读取输入的时间。鉴于这一有趣的事实，我们似乎应该将所有精力都投入到改进布隆过滤器的性能上。但是，出于几个原因，检查数据泵和解析器还是值得的。一个原因是我们需要修改现有数据泵和解析器的设计，以适应多线程的使用，从而实现可扩展性。另一个原因是我们有兴趣减少内存到内存的复制，这可能会影响布隆过滤器与其与解析器接口的效率。第三个原因是我们希望为积极的预读或数据预取做好准备，以防我们能够提高 k-mer 计数逻辑的效率，从而使输入时间与计数时间具有竞争力。与性能调优无关，在可维护性和可扩展性方面也存在一些问题。

事实证明，上述所有原因都汇聚到了一个新的设计上。我们将在后面更详细地讨论此设计的线程安全方面。现在，我们将重点关注减少内存到内存的复制以及执行相当积极的数据预取的能力。

通常，当程序从块存储设备（例如硬盘）检索数据时，操作系统会缓存一定数量的块，以防再次需要这些块。此缓存活动会带来一些时间开销；此外，无法对要预取到缓存中的数据量进行微调。此外，用户进程无法直接访问缓存，因此必须将其从缓存复制到用户进程的地址空间。这是一个内存到内存的复制。

某些操作系统，例如 Linux，允许对其预读窗口进行一些调整。例如，可以对特定文件描述符调用posix_fadvise(2)和readahead(2)。但是，这些允许的控制相当有限，并且不会绕过缓存。我们有兴趣绕过由操作系统维护的缓存。如果使用O_DIRECT标志打开文件并且文件系统支持该标志，则实际上可以绕过此缓存。使用直接输入并非完全简单，因为来自存储的读取必须是存储介质块大小的倍数，并且必须放置在内存区域中，该区域的基地址是相同块大小的倍数。这需要程序执行文件系统通常会执行的管理工作。我们实现了直接输入，包括必要的管理工作。但是，在某些情况下，直接输入将无法工作或不希望使用。对于这些情况，我们仍然尝试调整预读窗口。我们的存储访问是顺序的，我们可以通过使用posix_fadvise(2)提供提示来告诉操作系统读取比平时更远的内容。

最大程度地减少缓冲区到缓冲区的复制是数据泵和解析器之间共有的挑战。在理想情况下，我们将从存储中读取一次到我们自己的缓冲区，然后为每个基因组读取扫描一次我们的缓冲区，以在缓冲区内用偏移量和长度来界定一个序列。但是，管理缓冲区的逻辑足够复杂，并且解析逻辑（考虑到我们的特定细微差别）也足够复杂，以至于维护一个中间行缓冲区对于程序员理解非常有用。为了减少此中间缓冲区的影响，我们鼓励编译器相当积极地内联代码的这一部分。如果代码此特定区域的性能成为一个足够大的问题，我们最终可能会消除中间缓冲区，但这可能会以牺牲易于理解的软件设计为代价。

布隆过滤器操作

回想一下，我们正在处理由四个字母组成的字母表（A、C、G 和 T）组成的序列，您可能会问这些是大小写字母。由于我们的软件直接操作用户提供的数据，因此我们不能依赖数据始终是大写或小写，因为测序平台和其他软件包可能会更改大小写。虽然对于单个基因组读取来说这很容易解决，但我们需要对数百万或数十亿个读取中的每个碱基重复此操作！

在对代码进行性能调优之前，它对大小写不敏感，直到它验证 DNA 字符串和生成哈希码的点。在这些点上，它会冗余地调用 C 库的toupper函数以将序列规范化为大写，使用以下宏

#define is_valid_dna(ch) \
    ((toupper(ch)) == 'A' || (toupper(ch)) == 'C' || \
     (toupper(ch)) == 'G' || (toupper(ch)) == 'T')

和

#define twobit_repr(ch) \
    ((toupper(ch)) == 'A' ? 0LL : \
     (toupper(ch)) == 'T' ? 1LL : \
     (toupper(ch)) == 'C' ? 2LL : 3LL)

如果您阅读toupper函数的手册页或检查 GNU C 库的头文件，您可能会发现它实际上是一个区域感知函数，而不仅仅是一个宏。因此，这意味着涉及调用一个可能不简单的函数的开销——至少在使用 GNU C 库时是这样。但是，我们正在使用四个 ASCII 字符的字母表。区域感知函数对于我们的目的来说过于复杂。因此，我们不仅希望消除冗余，而且希望使用更有效的方法。

我们决定在验证序列之前将其规范化为大写字母。（当然，验证发生在尝试将其转换为哈希码之前。）虽然在解析器中执行此规范化可能是理想的，但事实证明，可以通过其他途径将序列引入布隆过滤器。因此，目前，我们选择在验证序列之前立即对其进行规范化。这使我们能够删除序列验证器和哈希编码器中对toupper的所有调用。

考虑到可能会有数 TB 的基因组数据通过序列规范化器，因此我们有兴趣尽可能地对其进行优化。一种方法是

#define quick_toupper( c ) (0x60 < (c) ? (c) - 0x20 : (c))

对于每个字节，上述内容应执行一次比较、一次分支和可能一次加法。我们能做得更好吗？事实证明，可以。请注意，每个小写字母都有一个 ASCII 代码，该代码比其大写对应字母大 32（十六进制 20），并且 32 是 2 的幂。这意味着 ASCII 大写和小写字符仅在一个位上有所不同。这个观察结果表明“位掩码！”

c &= 0xdf; // quicker toupper

上述内容有一个按位运算，没有比较，也没有分支。大写字母原样通过；小写字母变成大写。完美，正如我们所期望的。由于我们的努力，我们在整个过程的运行时获得了大约 13% 的加速（！）

我们的布隆过滤器的哈希表……“庞大”。要增加特定 k-mer 的哈希码的计数，意味着要访问接近 N 个不同的内存页面，其中 N 是分配给过滤器的哈希表数量。在许多情况下，需要为下一个 k-mer 更新的内存页面与当前 k-mer 的内存页面完全不同。这可能导致从主内存中循环访问大量内存页面，而无法利用缓存的优势。如果我们有一个具有 79 个字符长序列的基因组读取，并且正在扫描长度为 20 的 k-mer，并且我们有 4 个哈希表，那么最多可能触及 236（59 * 4）个不同的内存页面。如果我们正在处理 5000 万个读取，那么很容易看出这有多么昂贵。该怎么办呢？

一个解决方案是批量更新哈希表。通过累积大量不同 k-mer 的哈希码，然后定期使用它们以逐表的方式递增计数，我们可以大大提高缓存利用率。这方面初步的工作看起来很有希望，并且希望在您阅读本文时，我们将已将此修改完全集成到我们的代码中。尽管我们在前面讨论测量和分析时没有提到它，但cachegrind（Valgrind [7] 发行版中的一部分开源程序）是一个非常有用的工具，可用于衡量此类工作的有效性。

并行化

随着当今世界多核架构的激增，尝试利用它们非常诱人。但是，与许多其他问题领域（如计算流体动力学或分子动力学）不同，我们的“大数据”问题依赖于数据的高吞吐量处理——它必须在并行化的某个点之后基本上成为 I/O 绑定。超过此点，向其投入更多线程无济于事，因为对存储介质的带宽已饱和，线程最终只会增加阻塞或 I/O 等待时间。也就是说，利用一些线程可能很有用，特别是如果要处理的数据存储在物理 RAM 中，因为 RAM 通常比在线存储具有更高的带宽。如前所述，我们已将预取缓冲区与直接输入结合使用。多个线程可以使用此缓冲区；下面将详细介绍。I/O 带宽不是多个线程必须共享的唯一有限资源。用于 k-mer 计数的哈希表是另一个。下面也将讨论对这些哈希表的共享访问。

线程安全和线程

在继续详细介绍之前，澄清一些关于术语的要点可能会有所帮助。人们经常将某物是线程安全的概念与某物是多线程的概念混淆。如果某物是线程安全的，那么多个线程可以同时访问它，而无需担心获取或存储损坏。如果某物是多线程的，那么它将由多个执行线程同时操作。

作为我们并行化工作的一部分，我们重构了 C++ 核心实现的部分内容，使其成为线程安全的，而无需对特定的线程方案或库做出任何假设。因此，Python threading模块可用于使用核心实现周围的 Python 包装器的脚本，或者核心周围的 C++ 驱动程序可以使用更高级别的抽象（如我们前面提到的 OpenMP）或显式地使用 pthreads 实现线程，例如。实现这种线程方案的独立性并保证线程安全，同时又不破坏对 C++ 库的现有接口，这是一个有趣的软件工程挑战。我们通过让作为线程安全公开的 API 部分维护自己的每个线程状态对象来解决此问题。这些状态对象在一个 C++ 标准模板库 (STL) map中查找，其中线程标识号是键。特定线程的标识号是通过让该线程本身通过系统调用查询操作系统内核来找到的。此解决方案确实会引入少量开销，因为线程在每次进入通过 API 公开的函数时都会查询其标识号，但它巧妙地避免了破坏现有接口的问题，这些接口是在考虑单个执行线程的情况下编写的。

数据泵和解析器操作

在高性能计算（HPC）领域遇到的多核机器可能有多个内存控制器，其中一个控制器比另一个CPU更靠近（就信号传输距离而言）一个CPU。这些是非一致内存访问（NUMA）架构。使用这种架构的机器带来的一个影响是，内存获取时间可能因物理地址而异。由于生物信息学软件通常需要大量的内存才能运行，因此通常会在这些机器上运行。因此，如果使用多个线程（可能固定到各种NUMA节点），则必须考虑物理RAM的局部性。为此，我们将预取缓冲区划分为多个段，段数等于执行线程数。每个执行线程负责分配其特定的缓冲区段。缓冲区段通过一个状态对象进行管理，该状态对象在每个线程的基础上维护。

布隆过滤器操作

Bloom过滤器哈希表消耗了大部分主内存（参见图12.1），因此无法在线程之间有效地拆分为单独的副本。相反，所有线程必须共享一组表。这意味着线程之间会争夺这些资源。需要内存屏障[8]或某种形式的锁定来防止两个或多个线程同时尝试访问同一内存位置。我们使用原子加法运算来递增哈希表中的计数器。这些原子操作[9]在许多平台上都受到多个编译器套件的支持，其中包括GNU编译器，并且不依赖于任何特定的线程方案或库。它们在要更新的操作数周围建立内存屏障，从而为特定操作添加线程安全性。

一个我们没有解决的性能瓶颈是，在k-mer计数完成后将哈希表写入存储的时间。我们不认为这是一个优先级很高的任务，因为对于给定的Bloom过滤器大小，写入时间是恒定的，并且不依赖于输入数据量。对于我们用于基准测试的特定5 GB数据集，我们发现k-mer计数花费的时间是哈希表写入时间的六倍多。对于更大的数据集，这种比率会变得更加明显。也就是说，我们最终也希望在这里提高性能。一种可能性是在操作的k-mer计数阶段将写入成本分摊。URL缩短网站bit.ly有一个名为dablooms[10]的计数Bloom过滤器实现，它通过将其输出文件内存映射到哈希表内存来实现这一点。采用他们的想法，结合哈希表的批量更新，将有效地为我们提供在进程生命周期内以突发形式进行的异步输出，并从执行结束时切掉整个写入时间。但是，我们的输出不仅仅是计数表，还包括包含一些元数据的标题；考虑到这一点来实现内存映射是一项需要仔细思考和谨慎处理的工作。

扩展性

使khmer软件具有可扩展性是否值得我们付出努力？是的。当然，我们没有实现完美的线性加速。但是，对于每个核心数量翻倍，我们目前大约可以获得1.9倍的加速。

图12.4 - 从1到8个CPU核心的加速因子

在并行计算中，必须注意Amdahl定律[11]和收益递减规律。在并行计算的上下文中，Amdahl定律的常用公式为$S(N) = \frac{1}{(1 - P) + \frac{P}{N}}$，其中S是在给定N个CPU核心时实现的加速，而P是并行化的代码的比例。对于$\lim_{N\to\infty} S = \frac{1}{(1 - P)}$，这是一个常数。软件使用的存储系统的I/O带宽是有限的且不可扩展的；这导致了非零(1 − P)。此外，并行部分中对共享资源的争用意味着$\frac{P}{N}$实际上是$\frac{P}{N^l}$，其中l < 1，而理想情况是l = 1。因此，即使在有限数量的核心上，收益也会更快地递减。

使用更快的存储系统，例如固态驱动器（SSD）而不是硬盘驱动器（HDD），可以提高I/O带宽（从而减少(1 − P)），但这超出了软件的范围。虽然我们无法对硬件做太多的事情，但我们仍然可以尝试改进l。我们认为我们可以进一步改进围绕共享资源（例如哈希表内存）的访问模式，并且可以更好地简化每个线程状态对象的用法。做这两件事可能会让我们改进l。

结论

khmer软件是一个不断变化的目标。定期向其中添加新功能，我们正在努力将其集成到生物信息学界使用的各种软件堆栈中。像学术界中的许多软件一样，它最初是一个探索性的编程练习，并发展成为研究代码。正确性是并且一直是该项目的主要目标。虽然性能和可扩展性不能真正被视为事后想法，但它们优先于正确性和实用性。也就是说，我们在可扩展性和性能方面做出的努力已经取得了良好的成果，包括单线程执行的加速以及通过使用多个线程显着减少总执行时间的能力。考虑性能和可扩展性问题导致了数据泵和解析器组件的重新设计。展望未来，这些组件应该不仅能够从可扩展性中受益，而且还能够提高可维护性和可扩展性。

未来方向

展望未来，一旦我们解决了基本的性能问题，我们主要感兴趣的是扩展程序员的API，提供经过良好测试的用例和文档，以及提供用于集成到更大管道中的特征明确的组件。更广泛地说，我们希望利用低内存数据结构理论的进步来简化某些用例，并且我们也对研究某些更难处理的数据集挑战的分布式算法感兴趣，这些挑战在我们不久的将来将会出现。

khmer开发面临的一些其他问题包括扩展哈希选项以允许对单链DNA使用不同的哈希函数，以及添加滚动哈希函数以允许k  >  32。

我们期待继续开发此软件，并希望对分子生物学家和生物信息学家面临的大数据问题产生影响。我们希望您喜欢阅读有关科学中使用的一些高性能开源软件。

致谢

感谢Alexis Black-Pyrkosz和Rosangela Canino-Koning的评论和讨论。

参考文献

[1]Various，“大数据”。http://en.wikipedia.org/w/index.php?title=Big_data&oldid=521018481。

[2]C. T. Brown等人，“khmer：基因组数据过滤和分区软件”。http://github.com/ged-lab/khmer。

[3]Various，“布隆过滤器”。http://en.wikipedia.org/w/index.php?title=Bloom_filter&oldid=520253067。

[4]O. members，“OpenMP”。http://openmp.org。

[5]A. D. Malony等人，“TAU：调整和分析实用程序”。http://www.cs.uoregon.edu/Research/tau/home.php。

[6]Various，“概要引导优化”。http://en.wikipedia.org/w/index.php?title=Profile-guided_optimization&oldid=509056192。

[7]J. Seward等人，“Valgrind”。http://valgrind.org/。

[8]Various，“内存屏障”。http://en.wikipedia.org/w/index.php?title=Memory_barrier&oldid=517642176。

[9]Various，“原子操作”。http://en.wikipedia.org/w/index.php?title=Linearizability&oldid=511650567。

[10]b. software developers，“dablooms：一个可扩展的计数布隆过滤器”。http://github.com/bitly/dablooms。

[11]Various，“Amdahl定律”。http://en.wikipedia.org/w/index.php?title=Amdahl%27s_law&oldid=515929929。