CN111916153A - 一种并行多重序列比对方法 - Google Patents

Abstract

本发明公开了一种并行多重序列比对方法，该方法包括获取生物学序列数据，K‑mer特征提取，相似度计算，并行聚类分组，构建后缀树，双序列比对，融合比对结果。本发明通过K‑mer特征提取，采用并行聚类算法分组后构建后缀树，并采用Smith‑Waterman算法进行双序列比对，最后融合对比结果，能够较容易的实现多重序列并行化处理，从而快速完成多重序列比对，在序列数量多时相较其它比对方法效率突出，能够适应不同的序列特征，并节省用户大量的时间和精力。

Description

一种并行多重序列比对方法

技术领域

本发明属于计算机技术领域，具体涉及一种并行多重序列比对方法。

背景技术

多重序列比对是对三个以上的生物学序列(如蛋白质序列、DNA序列或RNA序列)所作的序列比对。一般来说，是输入一组假定拥有演化关系的序列。从多重序列比对的结果可推导出序列的同源性，而种系发生关系也可引导出这些序列共同的演化始祖。多重序列比对常用来研究序列的保守性，或是蛋白质结构域的三级结构与二级结构，甚至是个别的氨基酸或核苷酸。传统的多重序列比对方法通常过分依赖于动态规划而导致时间复杂度过高，用户在进行大量同源序列的比对时往往要花费数天甚至数月的时间。

发明内容

针对现有技术中的上述不足，本发明提供了一种并行多重序列比对方法。

为了达到上述发明目的，本发明采用的技术方案为：

一种并行多重序列比对方法，包括以下步骤：

S1、获取生物学序列数据；

S2、利用每个线程对各个序列进行K-mer特征提取；

S3、利用每个线程计算任意两个序列之间的相似度，得到相似度矩阵；

S4、采用并行聚类算法根据相似度矩阵对序列进行分组；

S5、在每个分组中选择与其它序列平均相似度最高的序列作为中心序列，对每个分组的中心序列分配一个线程构建后缀树；

S6、对每个分组或分组内的的一个非中心序列分配一个线程根据后缀树提取中心序列与其它每个序列的同源区段；

S7、对每个分组或分组内的的一个非中心序列分配一个线程采用Smith-Waterman算法对非同源子序列对进行比对，将比对结果和同源区段拼接得到双序列比对结果；

S8、对每个分组或分组内的的一个非中心序列分配一个线程将双序列比对结果融合；

S9、对每个分组分配一个线程根据融合后的比对结果构建每个分组的代表序列，得到多重序列比对结果。

进一步地，所述步骤S2具体包括：

提取每个序列中设定长度的全部连续子序列作为特征，并统计特征频次作为对应序列的特征值。

进一步地，所述步骤S3中计算任意两个序列之间的相似度具体包括：

选取两个序列对应特征中较小的特征值，计算选取的所有特征值的和，得到两个序列之间的相似度。

进一步地，所述步骤S4具体包括：

初始随机生成的k个均值点；

将每个观测分配到聚类中，使得组内平方和达到最小，表示为

其中，

表示第t趟迭代中第i个聚类，

分别表示第t趟迭代中第i和j个聚类对应的聚类中心，x_p表示第p个数据点；

根据得到的每一个聚类，以聚类中观测值的图心，作为新的均值点，表示为

其中，x_j表示第j个数据点。

进一步地，所述步骤S5具体包括：

在每个分组中选择一条序列与其它序列的相似度求和并计算平均值，得到该序列与其它序列的平均相似度；

选择该分组中平均相似度最高的序列作为中心序列，对每个分组的中心序列构建包含设定字符串所有长度后缀作为关键字、其位置作为值的压缩型的树。

进一步地，所述步骤S6具体包括：

从索引0开始，依次查找后缀树中是否存在后缀的前x个元素与当前序列[0:x]相同，若存在，则构成一个同源区段；否则，继续查找。

进一步地，所述步骤S7中采用Smith-Waterman算法对非同源子序列对进行比对具体包括：

设定非同源子序列对为A＝a₁ a₂ ... a_n和B＝b₁ b₂ ... b_n，其中n和m分别为序列A和B的长度，a和b分别表示序列A和B中的元素；

创建大小为n+1行m+1列的得分矩阵H，并初始化其首行和首列，表示为

H_k0＝H_0l＝0,(0≤k≤n,0≤l≤m)

其中，H_k0表示矩阵H第k行第0列的值，H_0l表示矩阵H第0行第l列的值，k和l分别表示矩阵H的行号和列号；

按照设定顺序进行打分，填充得分矩阵H剩余部分，表示为

其中，H_uv表示矩阵H第u行第v列的值，H_i-1,j-1+s(a_i,b_i)表示将a_i和b_j比对的相似性得分，H_i-k,j–W_k表示a_i位于一段长度为k的删除的末端的得分，H_i,j-l–W_l表示b_j位于一段长度为l的删除的末端的得分，0表示a_i和b_j到此为止无相似性；

从得分矩阵H中得分最高的元素开始根据得分的来源回溯至上一位置，如此反复直至遇到得分为0的元素。

进一步地，所述步骤S8具体包括：

对于每个非中心序列生成在中心序列中插入空格的数量数组和在当前非中心序列中插入空格的数量数组；

统计所有中心序列对应的空格数组，每个空格位置的数量取它们中的最大值，计算融合后中心序列的空格；

对于每一个非中心序列，计算中心序列对应空格数量数组的变化量，并计算需要在当前非中心序列插入空格数量数组中多加入空格的数量和位置，得到当前非中心序列最终插入空格的结果。

进一步地，所述步骤S9中每个分组的代表序列具体包括：

从融合后的比对结果中选择各个位上统计频数最多的字符构成对应分组的代表序列。

本发明具有以下有益效果：

本发明通过K-mer特征提取，采用并行聚类算法分组后构建后缀树，并采用Smith-Waterman算法进行双序列比对，最后融合对比结果，能够较容易的实现多重序列并行化处理，从而快速完成多重序列比对，在序列数量多时相较其它比对方法效率突出，能够适应不同的序列特征，并节省用户大量的时间和精力。

附图说明

图1为本发明的并行多重序列比对方法流程图；

图2为本发明实施例中后缀树构建过程示意图；

图3为本发明实施例中Smith-Waterman算法示意图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，本发明实施例提供了一种并行多重序列比对方法，包括以下步骤S1至S9：

S1、获取生物学序列数据；

在本实施例中，本发明通过读取生物学序列数据文件，获取待处理的生物学序列数据和对应的序列名；本发明读取的生物学序列数据文件为用于表示核酸序列或多肽序列的格式fasta，其中核酸或氨基酸均以单个字母来表示，且可以在每个序列前设置有序列名及注释。

S2、利用每个线程对各个序列进行K-mer特征提取；

在本实施例中，本发明提取每个序列中设定长度的全部连续子序列作为特征，并统计特征频次作为对应序列的特征值。具体而言，本发明提取每个序列中长度为2的子序列。长度为2的子序列在DNA的情况下一共有16种可能，这16种排列均可以作为序列的特征，其统计频数作为对应序列的特征值。此外，本发明还可以取长度为3或4的子序列以提升相似度的可靠性。

对于步骤S2，本发明为每个序列分配一个线程，并行完成计算。

S3、利用每个线程计算任意两个序列之间的相似度，得到相似度矩阵；

在本实施例中，本发明选取任意的两个序列对应特征中较小的特征值，计算选取的所有特征值的和，得到两个序列之间的相似度。

对于步骤S3，要计算n条序列两两之间的相似度，须进行n*(n–1)次计算，因此本发明为每次计算分配一个线程，并行完成计算。

S4、采用并行聚类算法根据相似度矩阵对序列进行分组；

在本实施例中，本发明可以采用k均值聚类算法利用相似度矩阵进行序列的分组，具体包括：

初始随机生成的k个均值点m⁽¹⁾ ₁,...,m⁽¹⁾ _n；

分配：将每个观测分配到聚类中，使得组内平方和达到最小，表示为

其中，

表示第t趟迭代中第i个聚类，

分别表示第t趟迭代中第i和j个聚类对应的聚类中心，x_p表示第p个数据点；并且每个x^p都只被分配到一个确定的聚类S^t中；

更新：根据得到的每一个聚类，以聚类中观测值的图心，作为新的均值点，表示为

其中，x_j表示第j个数据点；并且在对于观测的分配不再变化时收敛。

对于步骤S4，本发明采用并行k均值算法，并行完成计算。

S5、在每个分组中选择与其它序列平均相似度最高的序列作为中心序列，对每个分组的中心序列分配一个线程构建后缀树；

在本实施例中，本发明具体包括：

在每个分组中选择一条序列与其它序列的相似度求和并计算平均值，得到该序列与其它序列的平均相似度；

选择该分组中平均相似度最高的序列作为中心序列，对每个分组的中心序列构建包含设定字符串所有长度后缀作为关键字、其位置作为值的压缩型的树。

下面以"abcabxabcd"为例，对本发明中后缀树的构建方法进行详细说明，如图2所示：

前三个字符从左到右，对逐个字符进行操作。

第1个字符是'a'，创建一条边从根结点到叶结点，树更新如图2(a)。

处理第2个字符'b'，扩展已经存在的边"a"至"ab"；插入一条新边以表示"b"，树更新如图2(b)。

处理第3个字符'c'，重复同样的操作，树更新如图2(c)。

定义活动点为一个三元组，包括(活动结点，活动边，活动长度)；剩余后缀数为一个整数，代表着还需要插入多少个新的后缀。

此时活动点为(root,'c',1)。

处理第4字符'a'，已经存在一条边"abca"的前缀包含了后缀"a"，不再向活动结点插入一条新的边。相反，保留树原来的样子，如图2(d)。

更新活动点为(root,'a',1)，也就是说，活动结点仍为根结点，活动边为'a'，活动长度为1；剩余后缀数自增1得到2。

继续处理下一个字符'b'，保留树原来的样子，如图2(e)。

更新活动点为(root,'a',2)；剩余后缀数自增1得到3。

继续处理下一个字符'x'，此时活动边上接下去一个字符不是'x'，须分裂该边，树更新如图2(f)。

更新活动点为(root,'b',1)；剩余后缀数自减1得到2。

接下来处理后缀"bx"，此时活动边上接下去一个字符不是'x'，须分裂该边，树更新如图2(g)。

更新活动点为(root,'x',0)；剩余后缀数自减1得到1。该新结点不是当前步骤中创建的第一个结点，将先前加入的结点与该新结点通过一个特殊的指针连接，称为后缀连接，如图2(h)。

接下来处理后缀"x"，此时须从活动结点上发出新的分支，内容即为"x"，树更新如图2(i)。

更新活动点为(root,'x',0)；剩余后缀数不变。

继续处理下一个字符'a'，保留树原来的样子。

发现'a'已经存在于活动结点的一个边中，更新活动点为(root,'a',1)；剩余后缀数自增1得到2。

继续处理下一个字符'b'，保留树原来的样子。

更新活动点为(root,'a',2)；剩余后缀数自增1得到3。

继续处理下一个字符'c'，保留树原来的样子。

发现字符'c'存在于分支边上，则更新活动点为(node,'c',1)；剩余后缀数自增1得到4。

继续处理下一个字符'd'，从活动边上分裂出新的分支，内容即为"d"，树更新如图2(j)。

须沿着后缀连接的方向寻找结点，设置该结点为活动结点。(如果不存在，则设置根结点为活动结点。)更新活动点为(node,'c',1)；剩余后缀数自减1得到3。

接下来处理后缀"bcd"，从活动边上分裂出新的分支，内容即为"d"，树更新如图2(k)。

此时从非根结点中新发生一条边，且没有后缀连接，则更新活动点为(root,'c',1)；剩余后缀数自减1得到2。该新结点不是当前步骤中创建的第一个结点，将先前加入的结点与该新结点通过一个后缀连接相连。

接下来处理后缀"cd"，从活动边上分裂出新的分支，内容即为"d"，树更新如图2(l)。

更新活动点为(root,'d',0)；剩余后缀数自减1得到1。新结点不是当前步骤中创建的第一个结点，将先前加入的结点与该新结点通过一个后缀连接相连，如图2(m)。

接下来处理后缀"d"，从活动结点上发出新的分支，内容即为"d"，树更新如图2(n)。

对于步骤S5，本发明为每个组的中心序列分配一个线程，并使用后缀树并行构建算法完成计算。

S6、对每个分组或分组内的的一个非中心序列分配一个线程根据后缀树提取中心序列与其它每个序列的同源区段；

在本实施例中，本发明采用同源区段提取方法依次处理非中心序列外的其它序列。

同源区段提取方法具体包括：

从索引0开始，依次查找后缀树中是否存在后缀的前x个元素与当前序列[0:x]相同，其中x须大于一定阈值，如序列长度的5％；若存在，则构成一个同源区段；否则，继续查找。

对于步骤S6，本发明为每个组或组内的一个非中心序列分配一个线程，并行计算。

S7、对每个分组或分组内的的一个非中心序列分配一个线程采用Smith-Waterman算法对非同源子序列对进行比对，将比对结果和同源区段拼接得到双序列比对结果；

在本实施例中，本发明根据步骤S6可以得到中心序列和其它序列之一存在x个同源区段，则可将两序列分割为另外的x+1个非同源子序列对，使用Smith-Waterman(史密斯-沃特曼)算法进行比对，将比对结果和同源区段拼接，得到双序列比对结果。

采用Smith-Waterman算法对非同源子序列对进行比对具体包括：

设定非同源子序列对为A＝a₁ a₂ ... a_n和B＝b₁ b₂ ... b_n，其中n和m分别为序列A和B的长度，a和b分别表示序列A和B中的元素；

创建大小为n+1行m+1列的得分矩阵H，并初始化其首行和首列，表示为

H_k0＝H_0l＝0,(0≤k≤n,0≤l≤m)

其中，H_k0表示矩阵H第k行第0列的值，H_0l表示矩阵H第0行第l列的值，k和l分别表示矩阵H的行号和列号；

按照设定顺序进行打分，填充得分矩阵H剩余部分，如图3所示，表示为

其中，H_uv表示矩阵H第u行第v列的值，H_i-1,j-1+s(a_i,b_i)表示将a_i和b_j比对的相似性得分，H_i-k,j–W_k表示a_i位于一段长度为k的删除的末端的得分，H_i,j-l–W_l表示b_j位于一段长度为l的删除的末端的得分，0表示a_i和b_j到此为止无相似性；

从得分矩阵H中得分最高的元素开始根据得分的来源回溯至上一位置，如此反复直至遇到得分为0的元素。

对于步骤S7，本发明为每个组或组内的一个非中心序列分配一个线程，并行计算。

S8、对每个分组或分组内的的一个非中心序列分配一个线程将双序列比对结果融合；

在本实施例中，本发明具体包括：

首先对于每个非中心序列生成两个数组，包括在中心序列中插入空格的数量数组长度为中心序列长度加1)和在当前非中心序列中插入空格的数量数组(长度为当前非中心序列长度加1)；

然后计算融合后中心序列的空格，即统计所有中心序列对应的空格数组，每个空格位置的数量取它们中的最大值；

最后对于每一个非中心序列，计算中心序列对应空格数量数组的变化量，即在每个位置多加入多少空格，并计算需要在当前非中心序列插入空格数量数组中多加入空格的数量和位置，得到当前非中心序列最终插入空格的结果。

对于步骤S8，本发明为每个组或组内的一个非中心序列分配一个线程，并行计算。

S9、对每个分组分配一个线程根据融合后的比对结果构建每个分组的代表序列，得到多重序列比对结果。

在本实施例中，本发明从融合后的比对结果中选择各个位上统计频数最多的字符构成对应分组的代表序列，并将所有分组的代表序列构成一个新的分组，重复步骤S5至S8，得到最终的多重序列比对结果。

对于步骤S9，本发明为每个组分配一个线程，并行计算。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (9)

1.一种并行多重序列比对方法，其特征在于，包括以下步骤：

S1、获取生物学序列数据；

S2、利用每个线程对各个序列进行K-mer特征提取；

S3、利用每个线程计算任意两个序列之间的相似度，得到相似度矩阵；

S4、采用并行聚类算法根据相似度矩阵对序列进行分组；

S5、在每个分组中选择与其它序列平均相似度最高的序列作为中心序列，对每个分组的中心序列分配一个线程构建后缀树；

S6、对每个分组或分组内的的一个非中心序列分配一个线程根据后缀树提取中心序列与其它每个序列的同源区段；

S7、对每个分组或分组内的的一个非中心序列分配一个线程采用Smith-Waterman算法对非同源子序列对进行比对，将比对结果和同源区段拼接得到双序列比对结果；

S8、对每个分组或分组内的的一个非中心序列分配一个线程将双序列比对结果融合；

S9、对每个分组分配一个线程根据融合后的比对结果构建每个分组的代表序列，得到多重序列比对结果。

2.根据权利要求1所述的并行多重序列比对方法，其特征在于，所述步骤S2具体包括：

提取每个序列中设定长度的全部连续子序列作为特征，并统计特征频次作为对应序列的特征值。

3.根据权利要求1所述的并行多重序列比对方法，其特征在于，所述步骤S3中计算任意两个序列之间的相似度具体包括：

选取两个序列对应特征中较小的特征值，计算选取的所有特征值的和，得到两个序列之间的相似度。

4.根据权利要求1所述的并行多重序列比对方法，其特征在于，所述步骤S4具体包括：

初始随机生成的k个均值点；

将每个观测分配到聚类中，使得组内平方和达到最小，表示为

其中，

表示第t趟迭代中第i个聚类，

分别表示第t趟迭代中第i和j个聚类对应的聚类中心，x_p表示第p个数据点；

根据得到的每一个聚类，以聚类中观测值的图心，作为新的均值点，表示为

其中，x_j表示第j个数据点。

5.根据权利要求1所述的并行多重序列比对方法，其特征在于，所述步骤S5具体包括：

在每个分组中选择一条序列与其它序列的相似度求和并计算平均值，得到该序列与其它序列的平均相似度；

选择该分组中平均相似度最高的序列作为中心序列，对每个分组的中心序列构建包含设定字符串所有长度后缀作为关键字、其位置作为值的压缩型的树。

6.根据权利要求1所述的并行多重序列比对方法，其特征在于，所述步骤S6具体包括：

从索引0开始，依次查找后缀树中是否存在后缀的前x个元素与当前序列[0:x]相同，若存在，则构成一个同源区段；否则，继续查找。

7.根据权利要求1所述的并行多重序列比对方法，其特征在于，所述步骤S7中采用Smith-Waterman算法对非同源子序列对进行比对具体包括：

设定非同源子序列对为A＝a₁ a₂ ... a_n和B＝b₁ b₂ ... b_n，其中n和m分别为序列A和B的长度，a和b分别表示序列A和B中的元素；

创建大小为n+1行m+1列的得分矩阵H，并初始化其首行和首列，表示为

H_k0＝H_0l＝0,(0≤k≤n,0≤l≤m)

其中，H_k0表示矩阵H第k行第0列的值，H_0l表示矩阵H第0行第l列的值，k和l分别表示矩阵H的行号和列号；

按照设定顺序进行打分，填充得分矩阵H剩余部分，表示为

其中，H_uv表示矩阵H第u行第v列的值，H_i-1,j-1+s(a_i,b_i)表示将a_i和b_j比对的相似性得分，H_i-k,j–W_k表示a_i位于一段长度为k的删除的末端的得分，H_i,j-l–W_l表示b_j位于一段长度为l的删除的末端的得分，0表示a_i和b_j到此为止无相似性；

从得分矩阵H中得分最高的元素开始根据得分的来源回溯至上一位置，如此反复直至遇到得分为0的元素。

8.根据权利要求1所述的并行多重序列比对方法，其特征在于，所述步骤S8具体包括：

对于每个非中心序列生成在中心序列中插入空格的数量数组和在当前非中心序列中插入空格的数量数组；

统计所有中心序列对应的空格数组，每个空格位置的数量取它们中的最大值，计算融合后中心序列的空格；

对于每一个非中心序列，计算中心序列对应空格数量数组的变化量，并计算需要在当前非中心序列插入空格数量数组中多加入空格的数量和位置，得到当前非中心序列最终插入空格的结果。

9.根据权利要求1所述的并行多重序列比对方法，其特征在于，所述步骤S9中每个分组的代表序列具体包括：

从融合后的比对结果中选择各个位上统计频数最多的字符构成对应分组的代表序列。

Images