星码：基于全对搜索的序列聚类

2.2 poucet搜索算法

搜索策略可以进一步改进。如果两个连续查询共享一个长度为的前缀k个，计算的连续性k个对于这两个查询，编辑矩阵的第行将完全相同。因此，计算中间体可以存储在trie的节点中，以便下一个trie搜索可以从深度开始k个然而，将编辑矩阵的行存储在节点中会遇到一些困难。的确，在k个第行，对角线右侧的术语取决于两个查询之间不共享的字符。这个问题可以通过在每个节点中存储行和列术语的组合来解决，这些术语形成一个角度形状，看起来像水平翻转的L(图3). 使用此结构，计算中间层存储在深度节点中k个只取决于第一个k个查询的字符。

在单独的窗口中打开

图3。

Poucet搜索算法。该算法遵循与上所示相同的原理图2不同的是，编辑矩阵不是沿行更新，而是沿水平翻转的L更新。随着深度第一次搜索的继续，这些值存储在trie的节点中。由于翻转的L的垂直部分中的值对于节点的每个子节点都是相同的，因此只计算一次（箭头）。搜索路径访问节点时，将计算灰色单元格中的值。将中间产物存储在节点中可以让下一个查询在深度上重新启动k个如果它共享一个通用的长度前缀k个使用当前查询

为了充分利用此属性，输入序列按字母顺序排序，从而最大化连续查询之间的前缀共享。在童话故事《小普塞》中，主人公撒下白色鹅卵石，让哥哥们回家，这让人联想到一个较小的查询（按字母顺序）为下一个查询铺平了道路。因此，我们将此搜索算法称为“poucet”。

2.3无损过滤

当查询没有匹配项时，最好省略trie搜索。为此，starcode使用与下面描述的类似的分区方法Wu和Manber（1992）。查询最初以τ为单位进行分区 + 1段。假设所有片段的长度至少为τ，那么数据库中存在的每个τ-匹配将至少包含一个查询片段的逐字副本。事实上，在查询和要以τ分布的匹配之间最多有τ版本 + 1个区域，因此至少有一个段未修改。由于前面片段中可能存在插入和删除，共享片段可能会从其在查询中的原始位置向上移动到左侧（所有插入）或右侧（所有删除）的τ核苷酸。

这些观察结果是100%灵敏度过滤方法的基础。更准确地说，片段的定义如下：序列的第一个τ核苷酸被删除，序列的其余部分被划分为τ + 1段尺寸相差不超过1（为了一致性，较长的段总是在3'处）。每次向trie中添加序列时，都会对其进行分区，并将其分段添加到τ中 + 1个不同的索引。第一个片段被添加到第一个索引，第二个片段被添加到第二个索引，等等。在搜索之前，以相同的方式对查询进行分区，并在索引中查找其片段。如果没有找到匹配，则此查询在当前数据库中没有τ-匹配，因此省略trie搜索。相反，如果至少找到一个段，则必须执行trie搜索。

如上所述，可以发现查询和τ-匹配之间共享的片段向上移动为τ核苷酸。因此，根据图4，确保不会错过任何匹配：在τ中查找最右边的段 + 在τ-th指数中查找第一个指数、第二个最右边的段和移动了一个核苷酸的相邻段，依此类推，直到在第一个指数中查找移动了最多τ个核苷酸的第一个段及其相邻段。

在单独的窗口中打开

图4。

长度为20且τ为 = 3.删除查询的最后τ核苷酸，其余的被分为四个连续段。每个系列都根据编号为I–IV的不同索引进行查询。例如，针对索引I查询的唯一段是GTTG，而针对索引II查询的段是GCAA、CAAT和AATA。如果在适当的索引中找到任何一个段，则执行trie搜索，否则将忽略它，因为可能没有τ-匹配。无论结果如何，标记为I–IV的段都会添加到相应的相应索引中（即每个索引中只添加一个段）

2.4寻找和构建

为了减少搜索空间的大小，starcode使用了一种动态的“查找和构造”方法，在处理查询的同时构建trie。换句话说，每个序列在插入之前都会与trie进行匹配。如果A和B是相互τ-匹配的，则当B在trie中时将查询A，反之亦然。无论哪种方式，都会发现匹配的A-B。这样可以确保发现每个τ-匹配，同时尽可能保持trie的“瘦”，从而减少搜索时间。算法1和2所示的伪代码总结了整个匹配过程。

算法2 Poucet搜索算法-

1:程序 水龙头(查询,节点,种子,击打,鹅卵石):

2:计算 节点-NW后面的特定列◃图1

​三：3:为所有人 小孩中的节点节点 做

4:计算 小孩-NW后面的特定行◃图1

​5:5:计算使用行和列的中心值◃图1

在单独的窗口中打开

图5。

可扩展性。(一)运行时间的对数与待聚类序列数的对数。(b条)运行时间是聚类距离的函数。(c（c）)运行时间与输入序列的长度。(d日)不同数量的并行线程的相对性能提高

​6:6:如果中心值>τ 然后◃ 超过不匹配。

在单独的窗口中打开

图6。

已知集群结构的人工数据集的基准测试结果（见正文）。(一)通过已识别簇的数量来测量精确度。星码识别正确的簇数，而种子和cd-hit识别每个真阳性约40个假阳性。绘制了第一个平分线，并指示了完美的结果。(b条)通过识别对的数量测量精度。Slidesort识别的对比星码少5-10%。水平线表示比率为1。(c（c）)不同工具的运行时间。随着集群的大小，星码的运行时间增加，但仍具有竞争力。(d日)不同工具的内存使用情况。星码的内存使用量随着集群大小的增加而减少。

7:持续与下一个小孩

8:结束 如果

9:如果 节点深度 = 高度 然后◃ 找到命中。

10:节约 节点中的顺序击打

11:持续与下一个小孩

12:结束 如果

13:如果 节点深度≤秒e（电子）e（电子）d日 然后

14:节约 节点在里面鹅卵石在当前深度

15:结束 如果

16:呼叫水龙头(查询,小孩,种子,击打,鹅卵石)

17:结束 对于

18:结束 程序

2.5并联

查询在连续的块中进行排序和分区。然后，匹配步骤分两个阶段进行。在构建阶段，根据上述算法，从每个块的序列构建不同的trie。在第二种情况下，针对所有其他尝试查询所有序列块。如果查询在中进行分区N个块，第一阶段包括N个寻找和构建工作，而第二个工作包括N个(N个−1）/2个查询作业。在每个阶段，作业之间没有相互依赖性，因此匹配算法可以有效地并行化N个大于独立线程的数量。

2.6聚类

星码的默认聚类算法是为了纠正排序错误而设计的。此方法使用消息传递(麦凯，2002年)识别和统计“规范”序列（在聚类术语中也称为质心）。默认情况下，每个序列都会将其读取计数转移到最接近的τ-匹配，前提是后者的计数至少多出五倍。如果不满足条件，则不会进行转移。如果序列有多个相等接近的τ-匹配，则计数在它们之间平均分配。该过程以递归方式重复，从读取计数最低的序列开始。在过程结束时具有正读取计数的序列被认为是规范的。簇由所有将其读取计数转移到同一规范序列的序列组成（将其读取数转移到不同规范的序列被丢弃）。请注意，簇的半径可以大于用于匹配的最大距离。

由于没有任何测序技术的错误率高于20%，预计测序错误产生的序列的读取计数将始终是标准序列的五倍或更低。否则，序列很可能是不相关的，或者两者都是从同一个规范序列派生的。可以使用命令行选项修改此行为群集比率允许更灵活或更严格的聚类，例如将唯一的输入序列聚类在一起，群集比率必须设置为1。

对于其他序列聚类问题，starcode实现了一种称为“球体聚类”的多用途算法。在球体聚类中，序列按出现频率排序。从最常见的序列开始，每个序列都成为规范序列并声明其所有τ-匹配，从而形成半径τ的簇（因此而得名）。声明的序列会立即删除，因此它们只能属于一个簇。

3.2基准

我们根据序列聚类算法slidesort对星码进行了基准测试(清水和津田，2011年)，种子(包et（等） 阿尔。, 2011)，彩虹(Chong（冲）et（等） 阿尔。, 2012)和cd-hit(傅et（等） 阿尔。, 2012). 尽管slidesort是一种全对搜索算法，但它包含在基准测试中，因为序列比较是序列聚类问题中计算量最大的步骤。Rainbow只在配对读取上运行，而其他工具在单个读取上运行。因此，所有工具都不能在同一个数据集上运行。

序列聚类算法的性能对数据集中簇的大小很敏感，在许多应用中，簇的大小是未知的。因此，我们在人工数据集上建立了一个基准，以测试工具在已知簇结构上的准确性和伸缩性。我们生成了四个数据集，共有100万50万人，分布在1到1000个集群中。每个星团由相同质心序列的100个重复组成，再加上通过合并三个误差（最多一个indel）从质心衍生的卫星。每个集群的卫星数量从999颗不等 900至900。Rainbow没有在此基准测试中进行测试，因为生成的数据是单次读取。此外，评估slidesort的准确性存在问题，因为每个数据集中3个匹配的数量未知（同一集群中的卫星对可能是3个匹配，也可能不是）。由于这个原因，我们只比较了星码找到的对的数量和slidesort找到的对数量。测试结果总结如下图6.

虽然starcode在所有四个数据集上都实现了完美的聚类，但通过seed和cd-hit实现的聚类是不完整的。这两种工具在所有数据集的每个真实集群上都可以识别大约40个假集群(图6a） 。我们还观察到，slidesort在所有数据集中找到的3个匹配项比starcode少5-10%(图6b） ●●●●。我们对这个结果感到惊讶，因为slidesort被认为是一个精确的算法。然而，当我们在较小的数据集上运行额外的测试时，情况显然不是这样，在这些数据集上，朴素的两两比较是可行的（有关星码存储库的更多信息，请参阅http://github.com/gui11aume/starcode/tree/master/misc).

不同工具的运行时间与集群大小的函数如所示图6c.slidesort和starcode的运行时间分别显示线性和亚线性趋势。种子和cd-hit大约在恒定时间内运行，与集群大小无关。尽管有这样的结果，星码的性能仍然具有竞争力，即使对于100万个序列的集群也是如此。内存使用情况如所示图6c.通过slidesort和cd-hit实现最小内存占用，与其他工具的最大差异为一个数量级。请注意，与slidesort进行比较并不完全公平，因为它没有在内存中保存集群所需的完整图形。对于1000大小的集群，星码的内存使用率最高，但随着集群大小的增加，星码内存使用率会降低，并低于种子的内存使用量。总之，星码是以增加内存占用为代价在这些数据集上实现完美精度的唯一工具。考虑到输出的准确性，星码在运行时间方面保持了具有竞争力的性能。

人工数据的性能并不总是与实验数据集的性能一致。典型的实验带来了额外的困难。例如，集群的大小可能不均匀，读取可能包含近似恒定的区域，这些区域通常会降低基于过滤器的算法的性能。我们对序列聚类算法进行了基准测试，以解决将数千名记者聚集在并行（TRIP）条码中的问题(阿赫塔et（等） 阿尔。, 2013). 简言之，TRIP的原理是用随机条形码标记报告转录本，并测量RNA中条形码的丰度，作为基因表达的代理。由于标签序列未知，因此没有参考来匹配异常条形码。

用于基准测试的数据集的基本属性总结如下表3。已对数据集1（SRR950457）进行预处理，以提取条形码并删除读取的常量部分。文件中仅包含15到17个核苷酸之间的条形码。数据集2（PRJEB7686）由原始Illumina单次读取组成。这些数据集因读取大小、总读取计数和经验簇大小而异。根据starcode的输出，数据集1的最大集群包含大约70个 000个序列，而数据集2包含四个具有超过100万个序列的簇。数据集3（SRR950477）已被包括在内，用于在成对聚类模式下根据彩虹对星码进行基准测试。

starcode、seed、slidesort、rainbow和cd-hit的运行时间总结如下表4。我们调整了第一个数据集的距离阈值，以补偿减少的序列长度。starcode和slidesort都是使用选项“-d日 2'而cd-hit的标识设置为'-c（c） 0.85’. 由于最小序列长度的限制，我们无法在数据集1上运行seed。Starcode在所有数据集上都比其他工具快得多。Seed和cd-hit位居第二，在数据集1和2上的运行时间分别约为35倍和20倍。Rainbow在成对读取的集群工作中几乎慢了一个数量级。我们没有记录过去10天的准确运行时间，因为这比星际码的运行时间高几个数量级。

相同数据集上不同工具的内存占用如所示表5。这些值表示在上述数据集上运行的整个过程中的峰值内存使用量。在短读（数据集1）中，starcode利用trie压缩优于其他工具。在数据集3上，星码的内存使用量明显大于彩虹。Starcode和cd-hit在数据集2上使用了相似的内存量。两者所需的内存是slidesort的两倍，slidesoort的优点是在全对比较期间不存储完整的图形。

我们要感谢玛丽亚·查卓（Maria Chatzou）对初稿的宝贵反馈，并感谢陈恒昌（Heng-Chang Chen）执行了果蝇属TRIP实验。