基于生成树的大规模生物数据聚类方法及系统

1.本发明属于数据处理技术领域，尤其涉及基于生成树的大规模生物数据聚类方法及系统。


背景技术：

2.本部分的陈述仅仅是提供了与本发明相关的背景技术信息，不必然构成在先技术。
3.随着基因测序技术的增长和测序成本的降低，生物基因组数据的规模越来越大，并且整体规模呈指数型增长。著名的基因组数据库ncbi的refseq数据库的参考基因组数据规模已经达到了tb级别。在不久的将来，可能会达到pb甚至更高的级别。针对大规模数据集，需要有相对应的算法解决大规模数据问题。
4.聚类算法在生物大数据处理中是一种常用的算法，它的主要途径是将相似度较高的基因组数据聚集到一起，从而达到去重或者揭示物种边界的目的。当前也已经存在很多像cd-hit，uclust等比较常用的生物序列聚类算法工具，虽然很多已有的聚类工具是基于启发式的思想，能够减少一部分计算从而提高工具整体的性能，但是在面对大规模的超长的基因组序列的情况时，这些聚类算法仍然不能在有限的时间和内存占用中完成聚类，这迫切需要新的算法来实现针对大规模基因组数据的聚类。
5.在聚类过程中，基因组序列之间的相似度计算尤为重要。较为传统的生物基因组序列相似度计算算法是局部精确比对算法史密斯-沃特曼算法（smith-waterman algorithm）和全局精确比对算法尼德曼-翁施算法（needleman-wunsch algorithm）。这种经典的基于精确比对的算法的通过计算得分矩阵获取序列间的相似度，计算复杂度非常高，在基因组序列长度为n的情况下，复杂度为o(n2)。在针对较长的基因组序列的情况下，很难在有效的时间内计算出相似度。
6.因此，需要解决的技术问题主要为：生物基因组序列的长度越来越长，传统的序列比对算法不能够在有效的时间内完成对序列相似度的计算。
7.生物序列的数目越来越多，针对全部序列两两之间相似度度量的距离矩阵不能够在有限的内存空间内存储。
8.在当前多核平台处理器上，很多应用不能够充分利用和发挥多核平台并行计算的优势。


技术实现要素：

9.为克服上述现有技术的不足，本发明提供了基于生成树的大规模生物数据聚类方法，采用sketch算法估计序列之间的相似度，利用最小生成树从而避免了整个距离矩阵在内存中的存储。
10.为实现上述目的，本发明的一个或多个实施例提供了如下技术方案：
第一方面，公开了基于生成树的大规模生物数据聚类方法，包括：估计基因组序列之间的相似度；基于估计基因组序列之间的相似度流式计算基因组序列之间的距离矩阵和构建最小生成树，通过将距离矩阵划分成子图并且构建子最小生成树来合并生成最终的最小生成树；剪断最小生成树中超过给定阈值长度的边来生成聚类结果。
11.进一步的技术方案，采用sketch算法估计基因组序列之间的相似度。
12.进一步的技术方案，估计基因组序列之间的相似度的具体步骤为：将原始的基因组序列通过滑动窗口的方式分解成长度为k的子序列，称为k-mer，然后保留部分k-mer；以部分k-mer的集合来代表整条序列，从而将直接计算原始基因组序列的相似度转化为用sketch中保存的k-mer的集合之间的相似度计算来估计原始序列的相似度。
13.进一步的技术方案，将原始的生物基因组序列通过滑动窗口的方式产生k-mer，得到k-mer的集合，然后将集合中的k-mer通过哈希函数映射成相应的哈希值，通过minhash方法选取其中最小的固定个数哈希值保存作为sketch，其中固定个数就是sketch的维度。
14.进一步的技术方案，在计算距离矩阵的过程中，将整个距离矩阵所构成的完全图划分为一个个的子图，对每个子图构建子最小生成树，并且在生成距离矩阵的过程中不断更新子最小生成树，最终将子最小生成树进行合并得到最终的最小生成树。
15.进一步的技术方案，所述估计基因组序列之间的相似度之前还包括：采用单个生产者多个消费者的模型进行并行加速读取基因组序列文件并对其进行格式化处理。
16.进一步的技术方案，采用sketch算法估计基因组序列之间的相似度时，在sketch的构建过程中，对于多个基因组sketch的构建是通过多线程并行操作实现。
17.进一步的技术方案，流式计算距离矩阵和构建最小生成树的过程是通过多线程并行实现。
18.第二方面，公开了大规模生物数据聚类系统，包括：相似度估计模块，用于估计基因组序列之间的相似度；最小生成树生成模块，用于基于估计基因组序列之间的相似度流式计算方式来计算基因组序列之间的距离矩阵和构建最小生成树，通过将距离矩阵划分成子图并且构建子最小生成树来生成最小生成树；聚类模块，用于剪断最小生成树中超过给定阈值长度的边来生成聚类结果。
19.进一步的技术方案，所述相似度估计模块中估计基因组序列之间的相似度的具体步骤为：将原始的基因组序列通过滑动窗口的方式分解成长度为k的子序列，称为k-mer，然后保留部分k-mer；以部分k-mer的集合来代表整条序列，从而将直接计算原始基因组序列的相似度转化为用sketch中保存的k-mer的集合之间的相似度计算来估计原始序列的相似度。
20.以上一个或多个技术方案存在以下有益效果：本发明采用sketch算法估计序列之间的相似度，由于sketch中k-mer集合的维度
远小于原始序列的维度，所以，采用sketch算法进行序列相似度分析的计算时间和空间占用远远小于直接进行原始数据精确比对。
21.本发明距离矩阵的计算和子最小生成树的生成是通过流式计算生成的，从而避免了整个距离矩阵在内存中的存储。
22.本发明附加方面的优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
23.构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
24.图1为本发明实施例大规模生物基因组高效聚类工具的整体流程图。
具体实施方式
25.应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
26.需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。
27.在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
28.实施例一本实施例公开了大规模生物数据聚类方法，包括：步骤一：sketch的构建。这一步主要是将原始的生物基因组序列通过滑动窗口的方式产生k-mer得到k-mer的集合，然后将集合中的k-mer通过哈希函数映射成相应的哈希值，通过minhash方法选取其中最小的固定个数（1000）哈希值保存作为sketch，其中这个固定个数就是sketch的维度。由于将k-mer映射成哈希值的哈希函数满足均一性，即k-mer所映射成的哈希值在相应地哈希值空间中是均等分布的，所以，选取的最小的固定个数的哈希值等价于在所有的k-mer中随机抽取的固定个数的k-mer。并且所有的基因组序列采用同样的哈希函数，这样保证同样的k-mer输入有相同的哈希值输出，从而可以通过sketch之间的相似度计算来估计原始序列之间的相似度。序列之间的距离与相似度呈负相关，计算所得的两两序列之间的相似度用于生成包含所有序列之间距离的距离矩阵。其中，两条序列的距离d与他们之间的相似度j之间的关系最简单的可以理解为：d = 1-j。
29.值得注意的是，在sketch构建的过程中，多个基因组序列sketch的构建过程是通过多线程并行实现的。
30.步骤二：流式计算序列之间的距离矩阵和构建最小生成树。距离矩阵的构建和最小生成树的生成过程通过将距离矩阵划分成子图并且构建子最小生成树来进行的。所有的序列需要计算两两之间的相似度，通过对于sketch中选取的哈希值进行排序，然后对于sketch中的哈希值集合求交集和并集，通过交集和并集的来计算两个集合之间的jaccard相似度，两个集合之间的jaccard 相似度定义为：
其中a和b为两个sketch中哈希值所构成集合。另外由于序列变异位点的分布满足泊松分布，所以序列之间的变异距离定义为：其中j代表jaccard相似度，k代表k-mer的大小。所有的序列两两之间的变异距离构成最终的距离矩阵。多序列之间的距离矩阵的计算和子最小生成树的构建是通过多线程并行实现的。
31.本发明中，哈希函数为向量化的murmurhash3，计算距离的过程是集合求交集和并集，从而来度量jaccard相似度，进而得到变异距离（mutation distance）。
32.步骤三：最终聚类结果的生成。通过剪断最小生成树中超过给定阈值长度的边，最小生成树被分割形成由多个子树组成的森林，其中每个联通的子树为一个类，从而生成最终聚类结果，见附图1第三部分。
33.第一部分：输入的数据为标准的fasta格式的生物基因组序列数据文件，通过高效的格式化工具rabbitfx或者kseq将基因组数据格式化。然后将基因组序列碱基内容通过滑动窗口的方式产生k-mer的相关集合，将所有的k-mer通过simd指令集向量化的murmurhash3哈希函数将所有的k-mer转化成相应的哈希值集合，然后从中选取最小的固定数量（1000）的哈希值集合得到sketch，每个基因组得到一个sketch，并对其中保存的哈希值进行排序以便于之后jaccard相似度的计算。
34.第二部分：将所得到的sketch两两配对，对其中的保存的哈希值集合进行集合求交集和并集，得到他们的jaccard相似度和变异距离，所有sketch之间变异距离构成最终的距离矩阵。需要注意的是，sketch之间哈希值集合的变异距离的计算和最小生成树的构建是通过多线程流式进行的，每个线程在完成距离矩阵的t行距离的计算之后，就会生成或者更新子最小生成树，而非完成整个距离矩阵的计算之后才去构建最小生成树，从而避免了保存整体的距离矩阵并且便于多线程并行实现，最终多个线程将各自的子最小生成树合并得到最终的最小生成树。
35.子最小生成树合并最小生成树，利用克鲁斯卡尔算法生成最终最小生成树。
36.第三部分：将得到的最小生成树中超过一定阈值的边剪断，从而最小生成树变成包含多棵子树的森林，每棵子树为一个类，从而获得最终的聚类结果。
37.本发明的聚类方法可以通过一台拥有48核的处理器，2tb固态硬盘的服务器上在24分钟内完成对整个美国国家生物信息中心数据库（ncbi）中的refseq基因组数据集（1.2tb）进行聚类聚类分析，速度较快，并且能够得到同样精确的聚类结果。而对于其他已有的cd-hit和uclust等传统聚类工具，是不能够在有限时间内完成对如此大规模数据的聚类工作的。
38.在处理较长的基因组数据的情况下，基于非精确比对的算法得到发展。其中，sketch算法属于非精确比对算法的一种特殊形式。对于较长的基因组原始序列，sketch算法构建出简明扼要的序列信息保存为sketch，通过对sketch之间的相似度计算来估计原始
序列之间的相似度。由于sketch的大小相对于原始序列来讲非常小，所以，sketch算法可以在有效误差范围内和有限内存空间内高效的完成大规模基因组序列相似度的计算。
39.针对生物数据的聚类算法一般基于所有序列两两之间相似度所构成的距离矩阵来实现的。距离矩阵的维度和基因组序列的数目呈平方的关系，当基因组序列的数目越来越多，达到百万甚至更高级别的情况下，存储整个距离矩阵的内存开销就会变得非常大。最小生成树生成算法是图算法的一种，它的内存占用与输入元素的个数呈线性相关，即与基因组序列的数目呈线性关系，这将大大降低了内存占用。同时，最小生成树可以划分为不同的子树并且进行动态更新与合并，从而可以高效并行，这也将会大大提高程序的执行效率。
40.本发明整个聚类方法的设计与实现都是面向多核平台处理器，在主要的计算部分进行了多线程并行操作。在第一步的读取基因组序列文件并对其进行格式化的过程中，采用单个生产者多个消费者的模型进行并行加速功能，并且在sketch的构建过程中，对于多个基因组sketch的构建是通过多线程并行操作实现；流式计算距离矩阵和构建最小生成树的过程也是通过多线程并行实现。针对sketch中k-mer集合的交集计算，采用单指令多数据流（simd）指令集实现，大大提高了集合求交集的性能。
41.实施例二本实施例的目的是提供了大规模生物数据聚类系统，包括：相似度估计模块，用于估计基因组序列之间的相似度；最小生成树生成模块，用于基于估计基因组序列之间的相似度流式计算基因组序列之间的距离矩阵和构建最小生成树，通过将距离矩阵划分成子图并且构建子最小生成树来生成最小生成树；聚类模块，用于剪断最小生成树中超过给定阈值长度的边来生成聚类结果。
42.本领域技术人员应该明白，上述本发明的各模块或各步骤可以用通用的计算机装置来实现，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。本发明不限制于任何特定的硬件和软件的结合。
43.上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。