CN112182497B - 一种基于生物序列的负序列模式的相似性分析方法、实现系统及介质 - Google Patents

Abstract

本发明涉及一种基于生物序列的负序列模式的相似性分析方法、实现系统及介质，包括：(1)数据预处理：将DNA序列中的字母用数字来表示；并分割成若干个块，得到的若干个块作为频繁模式挖掘的数据集；(2)频繁模式挖掘：使用f‑NSP算法来挖掘数据集；(3)对最大频繁正、负序列模式进行图形表示；把最大频繁正、负序列模式转化为数字序列；(4)DNA序列的相似性分析：求取不同DNA序列的相似度，选取相似度最小的对应的DNA序列为待研究的DNA序列。本发明可以有效地对负序列进行有效的表达和分析，并且通过选取不同的最大频繁模式组合，能够得到不同的分析结果，大大节省了计算机的内存和时间的消耗。

Description

一种基于生物序列的负序列模式的相似性分析方法、实现系 统及介质

技术领域

本发明涉及一种基于生物序列的负序列模式的相似性分析方法、实现系统及介质，属于可决策的高效用负序列规则的应用技术领域。

背景技术

近年来，我们获得了海量的生物序列数据，随着DNA及蛋白质测序技术的进步，对解读生物序列数据中所含的各种信息，尤其是DNA序列中的遗传及调控信息、蛋白质序列结构与功能的关系的数据分析工具的需求增加，序列相似性分析得到广泛的应用。每当我们获得一个新的DNA序列的时候，希望通过相似性分析来证明它与某些已知的序列相似，如果和已知的序列具有同源性的话，会大大节省重新测定新序列的功能的时间和精力，而生物序列庞大，这就显得尤为重要了。在生物序列分析中，序列模式挖掘算法有助于识别同时发生的生物序列和发现DNA或蛋白质序列中的关系，因此研究缺失的碱基对序列比单一的去挖掘频繁的序列模式具有更高的意义。在生物信息学研究中，生物序列的相似性分析绝非简单机械的比较，而必然是多种多样的，同时还需要运用许多数学和统计学方法进行辅助分析与评判。序列相似性分析中，比对是最常用和最经典的研究手段。从生物序列的层次分析序列的相似性，推测其结构功能及进化上的联系是基因识别、分子进化、生命起源研究的基础，然而，在进行序列比对时有两方面的问题直接影响相似性分值：取代矩阵和空值罚分，粗糙的比对方法仅适用相同或不同来描述两个碱基的关系。生物序列的相似性分析用于提取储存在蛋白质序列中的信息，为此提出了许多数学方案。生物序列的图形表示可以识别任何序列的信息内容，以帮助生物学家选择另一种复杂的理论或实验方法。图形表示不仅提供了基因数据的可视化定性检查，而且还通过矩阵等对象提供了数学描述。大部分的数学方案是基于2-D和3-D表示的。

关于序列模式挖掘，对于正序列模式(Positive Sequential Pattern,PSP)挖掘仅仅考虑了已经发生了的事件(行为)，不同于传统序列模式挖掘的思路，负序列模式挖掘(Negative Sequential Pattern, NSP)还考虑了未发生的事件(行为)，也就是不存在于序列中的项，这样可以为人类提供更加全面地决策信息，比如，校园中存在的各种现状对学生的学习和生活产生的不同程度的影响；涉嫌医疗欺诈行为的参保人消除不良购药记录；缺失的基因片段可能诱发潜在的疾病等，但是，它们往往容易被人类所忽视，因此，越来越受到从事数据挖掘工作人员的关注。尤其是在生物序列分析中，序列模式挖掘算法有助于识别同时发生的生物序列和发现DNA或蛋白质序列中的关系，因此，研究缺失的碱基对序列比单一的去挖掘频繁的序列模式具有更高的意义。生物数据分析或生物数据挖掘存在一些重要的问题，如寻找共现的生物序列，对生物序列进行有效分类、对生物序列进行聚类分析等。而序列模式挖掘算法有助于识别同时发生的生物序列和发现DNA或蛋白质序列中的关系。生物序列数据往往包含着大量有价值的生物信息，例如，生物序列中频繁出现的基因和蛋白质片段往往含有许多未知的信息，挖掘这些信息具有重要的意义；某些细菌对人体的攻击受其基因中某些片段的影响；一些数目可变的串联重复序列的极度扩张可能会导致相关神经系统方面疾病。此外， DNA序列中的频繁模式的发现将是解释生物遗传特性的一种有效方法，这些频繁模式往往作为生物序列隐含数据的可能趋势和某些事件的相关标记。所以，在蛋白质或DNA等生物序列中频繁模式的挖掘具有重要价值。

目前存在的相似性分析方法，主要是针对PSP，对于前面我们挖掘出的NSP，尚缺少统一的相似性度量方法。而序列比对有一些缺点，促使人们试图寻找其他方法来比较DNA序列相似性。我们知道NSP在生物数据中的存在是不可避免的，甚至对一些致病基因至关重要。这就迫使我们找到一种方法，对缺失碱基序列的DNA进行相似性分析。

发明内容

针对现有技术的不足，本发明提出了一种基于生物序列的负序列模式的相似性分析方法；

本发明还提出了上述相似性分析方法的实现系统。

为了有效的分析DNA序列的相似性，应考虑解决以下关键问题:(1)如何用数字序列有效地表示 DNA主序列。(2)如何获得并选择合适的可视为DNA序列特征的描述符，根据数字序列对其进行表征。(3)如何有效处理不同长度的DNA序列，并保持其一致性。(4)如何对负序列进行有效的相似性分析。

术语解释：

1、DNA序列，又称基因序列，是使用一串字母表示的真实的或者假设的携带基因信息的DNA 分子的一级结构。

2、f-NSP算法，f-NSP使用位图来存储PSP数据，并通过位操作计算NSC支持度。它为size大于1的PSP创建位图，如果一条正序列被第i条数据序列包含，我们就将这条正序列的位图第i个位置置为1，反之则置为0。每一个位图的长度与数据序列所含有序列数相等。我们采用了新的位图存储结构，就可以使用位或(OR)操作来替代原有的并集操作。每个位图的长度等于数据库中序列的数量。假设s是一条正序列，它的位图使用B(s)来表示，获得位图中“1”的个数使用N(B(s))表示。则给出一个m-size并且n-neg-size的负序列ns，它的支持度是：

如果ns只包含一个负元素，那么序列ns的支持度是：

sup(ns)＝sup(MPS(ns))-sup(p(ns)) (2)

特别地，对于单元素负序列

f-NSP算法包括以下步骤。1.基于GSP算法从序列数据库中查找所有的PSP算法。所有的PSPs 和他们的位图被保存在一个散列表PSPHash中；2.使用NSC(负候选序列)生成方法为每个PSP 生成NSCs；3.使用公式(2)和(3)计算1-neg-size的nsc的支持度。通过公式(1)可以很容易地计算出其他nsc的支持度。具体地说，我们首先在1-negMSS_nsc中得到每个1-neg-MS’的位图。其次，使用OR操作获取位图的并集。然后，根据公式(1)计算nsc的支持度。最后，一个nsc是否是一个NSP是通过比较它与min_sup的支持度来确定的；4.返回结果并结束整个算法。

3、GSP算法，GSP算法是一种基于宽度优先搜索策略的挖掘算法，该算法通过一遍扫描数据库获得该数据库中包含的频繁项集，之后通过相应的连接和剪枝方法生成长度不断递增的候选序列，并基于重复扫描数据库的模式获得候选序列的支持度以判定正序列模式。GSP算法是典型的类 Apriori算法。在Apriori算法的基础上，GSP算法加入了分类层次、时间约束、滑动时间窗技术，使算法整体得到了优化。同时，GSP还对数据集的扫描条件进行了限定，降低了需要扫描的候选序列的数量，减少了无用模式的产生。

4、复平面，又称复数平面，即是z＝a+bi,它对应的坐标为(a,b)，其中，a表示的是复平面内的横坐标，b表示的是复平面内的纵坐标，表示实数a的点都在x轴上，所以x轴又称为"实轴"；表示纯虚数b的点都在y轴上，所以y轴又称为"虚轴"；y轴上有且仅有一个实点即为原点"0"。

5、嘌呤嘧啶图，简单来说就是在平面上绘制向量，将DNA序列中不同的碱基对准确的表示出来。这里我们是在复平面构造一个嘌呤嘧啶图，第一、二象限是嘌呤(A、

G和

)，第四象限是嘧啶(T、

C和

)。表示四个核苷酸A、G、C和及其对应的负序列的单位向量如下。这样不同的碱基对就能别唯一的表示，并且碱基对之间满足共轭关系。这种嘌呤嘧啶图符合DNA序列与它的时间序列是一一对应的特性。

6、DTW，Dynamic Time Warping，它出现的目的也比较单纯，首先在语音识别领域得到了广泛的应用，是一项把时间规划和距离测度结合起来的非线性规划技术，同时用于计算两个时间序列之间的最大相似性，即最小距离。

7、Apriori性质，任一频繁项集的所有非空子集也必须是频繁的。

本发明的技术方案为：

一种基于生物序列的负序列模式的相似性分析方法，包括步骤如下：

(1)数据预处理

对于每个要处理的序列或基因组，在将其进行频繁模式挖掘之前，都要进行预处理。将DNA 序列中的字母用数字来表示；由于DNA序列长度非常长，并将数字表示后的DNA序列分割成若干个块，每个块碱基数量相同，得到的若干个块作为频繁模式挖掘的数据集；

(2)频繁模式挖掘

使用f-NSP算法来挖掘数据集，得到最大频繁正、负序列模式；

(3)对最大频繁正、负序列模式进行图形表示

(4)DNA序列的相似性分析

求取不同DNA序列的相似度，相似度越小，DNA序列越相似。

相似性矩阵可以用来评价DNA相似性分析算法的有效性。它可以从侧面揭示不同物种之间的进化或遗传关系。DNA序列间的距离的计算是DNA相似性分析的基础，欧几里得距离和相关角是最常用的距离计算方法。并且规定序列之间的欧氏距离越小，DNA序列越相似。两个载体之间的相关角度越小，DNA序列越相似。

根据本发明优选的，步骤(2)中，使用f-NSP算法来挖掘数据集，数据集为D，包括步骤如下：

A、使用GSP算法得到所有的正频繁序列，并将每条正频繁序列对应的位图存入哈希表中；包括：

a、扫描数据集得到所有长度为1的序列模式放入原始种子集P₁中；

b、从原始种子集P₁中获取长度为1的序列模式，并将它们通过连接操作生成长度为2的候选序列集C₂；使用Apriori性质对候选序列集C₂进行剪枝，再通过扫描候选序列集C₂确定其中剩余序列的支持度，将支持度高于最小支持度的序列模式保存下来，输出为长度为2的序列模式L₂，并作为长度为2的种子集；用来生成长度递增的候选序列。按照该方法一直输出长度为3的序列模式 L₃、长度为4的序列模式L₄……长度为n+1的序列模式L_n+1，直到挖掘不出新的序列模式，得到序列模式即所有的正频繁序列，最小支持度是人为设置的支持度阈值min_sup；描述为：

L₁→C₂→L₂→C₃→L₃→C₄→L₄……若不能生成L_n+1停止。

B、基于所有正频繁序列生成相应的NSC；

NSC是指负候选序列，正频繁序列统称为正序列，为了从正序列中生成所有非冗余的NSC，生成NSC的关键过程是将具有正模式的非连续元素转换为它们的负伙伴，对于一个k-size的PSP， NSCs是通过将任意m个不相邻的元素改变为它的负数来生成的，用

来表示，

是不小于k/2的最小整数；k-size是指序列的大小为k；比如序列S＝{A T T C C}，其大小为 5-size。NSCs：指的是所有的负候选序列。

例如，<A T C C>的NSC包括：(1)m＝1时，

(2)m＝2时，

这里规定不允许出现两个连续的负项。

C、利用位操作快速计算负侯选序列的支持度；

生成NSCs之后，计算它的支持度，当满足负候选序列的支持度时，得到负频繁序列模式。NSCs的支持度计算如下：给定一个m-size和n-neg-size负序列ns，对于

那么在数据集D中ns的支持度为：

m-size是指序列大小为m；假设ns＝<a₁a₂…a_m>是一个负序列，若ns′仅由ns中所有的正元素组成，则将ns′称作ns的最大正子序列，定义为MPS(ns)；例如，

由这个序列的MPS(ns)和ns中一个负元素a组成的序列称为 1-neg-size最大子序列，定义为1-negMS。例如，

那么它的1-negMS就是

和

通过频繁模式挖掘，获得了12种最大频繁正、负序列模式；

根据本发明优选的，步骤(3)中，对最大频繁正、负序列模式进行图形表示，包括：在复平面构造一个嘌呤嘧啶图，嘌呤嘧啶图中，第一、二象限是嘌呤，包括A、

G和

第三、四象限是嘧啶，包括T、

C和

四个核苷酸A、G、T、C及其对应的负序列的单位向量

如式(Ⅰ)至式(Ⅷ)所示：

(b+di)→A(Ⅰ)

(d+bi)→G(Ⅱ)

(b-di)→T(Ⅲ)

(d-bi)→C(Ⅳ)

式(Ⅰ)至式(Ⅷ)中，b和d是非零的实数，

A和T是共轭的，G和C也是共轭的，即，

A、T、C，G代表现实存在的碱基对，

表示的是DNA序列中本该出现却没有出现的碱基对，又称缺失的碱基对，也叫A、G、T、C及其对应的负序列的单位向量；

通过这种表示方法，将一个DNA序列碱基

还原为一个数字序列s(n)，如式(Ⅸ)所示：

式(Ⅸ)中，s(0)＝0，其中y(j)满足式(Ⅹ)：

式(Ⅹ)中，j表示序列S中第0,1,2,...,n个位置上的碱基类型，n是被研究的DNA序列的长度；

通过上面的步骤，从“嘌呤嘧啶图”中唯一地获得原始的DNA序列的时间序列；

利用式(Ⅹ)把12种最大频繁正、负序列模式转化为数字序列；比如序列Human1通过式(Ⅸ) -(Ⅹ)得到的复数字序列为s(H1)＝{0.866+0.5i,1.366-0.366i,2.2321+0.134i,3.0981+0.634i,3.5981+1.5i, 4.4641+2i}，模组成的时间序列为S(H1)＝{1.0000,1.4142,2.2361,3.1623,3.8982,4.8916}。通过这样的方法可以得到12种频繁序列模式转化之后的时间序列。

根据本发明优选的，步骤(4)中，求取距离矩阵，距离矩阵用于表示不同DNA序列的相似度。

根据本发明优选的，步骤(4)中，通过DTW算法求取距离矩阵，设转化DNA序列而获得的时间序列为，

其长度分别为m和n；按照它们的时间位置进行排序，构造m×n矩阵A_m×n，矩阵中的每个元素

在矩阵中，把一组相邻的矩阵元素的集合称为弯曲路径，记为W＝w₁,w₂,...,w_k，W的第k个元素w_k＝(a_ij)_k，这条路径满足下列条件：

①max{m,n}≤K≤m+m-1；

②w₁＝a₁₁,w_k＝a_mn；

③对w_k＝a_ij,w_k-1＝a_i'j'必须满足0≤i-i'≤1,0≤j-j'≤1，则

DTW算法为运用动态规划思想寻找一条具有最小弯曲代价的最佳路径，如式(Ⅺ)所示：

其中，i＝2,3,...,m；j＝2,3,...,n。D(m,n)为A_m×n中弯曲路径的最小累加值。

上述相似性分析方法的实现系统，包括依次连接的数据预处理模块、频繁模式挖掘模块、图形表示模块、相似性分析模块；所述数据预处理模块用于执行步骤(1)；所述频繁模式挖掘模块用于执行步骤(2)；所述图形表示模块用于执行步骤(3)；所述相似性分析模块用于执行步骤(4)。

一种计算机可读存储介质，其特征在于，所述计算机可读存储介质中存储有基于生物序列的负序列模式的相似性分析程序，所述基于生物序列的负序列模式的相似性分析程序被处理器执行时，实现任一项所述的基于生物序列的负序列模式的相似性分析方法的步骤。

本发明的有益效果为：

1、本发明可以有效地对负序列进行有效的表达和分析，并且通过选取不同的最大频繁模式组合，能够得到不同的分析结果。

2、本发明选取的是频繁模式来进行的相似性分析，大大节省了计算机的内存和时间的消耗。

附图说明

图1为本发明基于生物序列的负序列模式的相似性分析方法的流程框图；

图2为本发明嘌呤嘧啶图的示意图；

图3为本发明基于生物序列的负序列模式的相似性分析方法的实现系统的结构框图；

图4为实施例中的位或(OR)运算过程示意图；

图5(a)为对最大频繁序列Human1,Opossum2,Rat2 and Chimpanzee2进行相似性分析后绘制的系统发育树示意图；

图5(b)为对最大频繁序列Human2,Opossum1,Rat2,and Chimpanzee1进行相似性分析后绘制的系统发育树示意图；

图6(a)为对最大频繁序列Human2,Opossum2,Rat2 and Chimpanzee1进行相似性分析后绘制的系统发育树示意图；

图6(b)为对最大频繁序列Human3,Opossu3,Rat3 and Chimpanzee3进行相似性分析后绘制的系统发育树示意图；

图7为归一化物种距离示意图。

具体实施方式

下面结合说明书附图和实施例对本发明作进一步限定，但不限于此。

实施例1

一种基于生物序列的负序列模式的相似性分析方法，如图1所示，包括步骤如下：

(1)数据预处理

对于每个要处理的序列或基因组，在将其进行频繁模式挖掘之前，都要进行预处理。将DNA 序列中的字母用数字来表示；由于DNA序列长度非常长，并将数字表示后的DNA序列分割成若干个块，每个块碱基数量相同，得到的若干个块作为频繁模式挖掘的数据集；

本发明中，每个序列首先被分割成几个块，每个块由相同数量的连续基组成。这些块是相互独立的，并且块的大小可以在实践中改变。注意，如果最后一个块的大小小于指定的块大小，那么这个块将被丢弃。为了更清楚，下面是一个关于分割块的例子。在本例中，有两个序列S₁和S₂。假设块大小为15，这两个序列分别被分成2个和3个块。最后一个大小为3的块被丢弃。其中每分一个块都被用曲线和直线标记出来。这也叫序列阻塞，它是一个重要的步骤，它带来两个主要的优点。首先，可以捕获序列的细粒度信息，包括位置信息和排序信息。其次，即使对于长序列，阻塞也可以减少序列处理的内存和时间消耗。

目前可以用于序列相似性研究的DNA序列很少，并且寻找更合适的DNA序列仍然是一个问题。来自15个物种的红蛋白基因的三个外显子序列是最常用的DNA序列。这三个基因序列包括第一、第二和第三外显子，序列的平均长度分别为92个碱基、222个碱基和114个碱基。其中， 11个不同物种的β基因的第一个外显子是应用最广泛的DNA序列数据。

所选数据集来自四个物种的β-蛋白基因的第一个外显子，如表1所示：

表1

(2)频繁模式挖掘

使用f-NSP算法来挖掘数据集，得到最大频繁正、负序列模式；

(3)对最大频繁正、负序列模式进行图形表示

(4)DNA序列的相似性分析

求取不同DNA序列的相似度，相似度越小，DNA序列越相似。

相似性矩阵可以用来评价DNA相似性分析算法的有效性。它可以从侧面揭示不同物种之间的进化或遗传关系。DNA序列间的距离的计算是DNA相似性分析的基础，欧几里得距离和相关角是最常用的距离计算方法。并且规定序列之间的欧氏距离越小，DNA序列越相似。两个载体之间的相关角度越小，DNA序列越相似。

实施例2

根据实施例1所述的一种基于生物序列的负序列模式的相似性分析方法，其区别在于：

步骤(2)中，使用f-NSP算法来挖掘数据集，数据集为D，包括步骤如下：

A、使用GSP算法得到所有的正频繁序列，并将每条正频繁序列对应的位图存入哈希表中；包括：

a、扫描数据集得到所有长度为1的序列模式放入原始种子集P₁中；

b、从原始种子集P₁中获取长度为1的序列模式，并将它们通过连接操作生成长度为2的候选序列集C₂；使用Apriori性质对候选序列集C₂进行剪枝，再通过扫描候选序列集C₂确定其中剩余序列的支持度，将支持度高于最小支持度的序列模式保存下来，输出为长度为2的序列模式L₂，并作为长度为2的种子集；用来生成长度递增的候选序列。按照该方法一直输出长度为3的序列模式 L₃、长度为4的序列模式L₄……长度为n+1的序列模式L_n+1，直到挖掘不出新的序列模式，得到序列模式即所有的正频繁序列，最小支持度是人为设置的支持度阈值min_sup；描述为：

L₁→C₂→L₂→C₃→L₃→C₄→L₄……若不能生成L_n+1停止。

使用图4来解释说明位或运算(OR)。序列S如果sup(s)≥min_sup，则称为频繁(正)序列模式，而如果sup(s)<min_sup，则称为不频繁序列模式。假设一条正频繁序列为<G C TA>，并且sup(C A)＝5，则根据负候选生成方法，一条负候选序列ns为

则相应的，MPS(ns) ＝<CA>，P(1-negMS₁)＝<GCA>，P(1-negMS₂)＝<C TA>。假设B(<G CA>)＝|1|0|0|1|0|，B(<C TA>)＝ |1|1|0|1|0|。那么B(<GCA>)ORB(<CTA>)的位图如图4所示。因此可以轻易的得到N(unionbitmap)＝4，之后由公式1得到

C、基于所有正频繁序列生成相应的NSC；

NSC是指负候选序列，正频繁序列统称为正序列，为了从正序列中生成所有非冗余的NSC，生成NSC的关键过程是将具有正模式的非连续元素转换为它们的负伙伴，对于一个k-size的PSP， NSCs是通过将任意m个不相邻的元素改变为它的负数来生成的，用

来表示，

是不小于k/2的最小整数；k-size是指序列的大小为k；比如序列S＝{A T T C C}，其大小为 5-size。NSCs：指的是所有的负候选序列。

例如，<A T C C>的NSC包括：(1)m＝1时，

(2)m＝2时，

这里规定不允许出现两个连续的负项。

C、利用位操作快速计算负侯选序列的支持度；

生成NSCs之后，计算它的支持度，当满足负候选序列的支持度时，得到负频繁序列模式。NSCs的支持度计算如下：给定一个m-size和n-neg-size负序列ns，对于

那么在数据集D中ns的支持度为：

m-size是指序列大小为m；假设ns＝<a₁a₂…a_m>是一个负序列，若ns′仅由ns中所有的正元素组成，则将ns′称作ns的最大正子序列，定义为MPS(ns)；例如，

由这个序列的MPS(ns)和ns中一个负元素a组成的序列称为 1-neg-size最大子序列，定义为1-negMS。例如，

那么它的1-negMS就是

和

通过频繁模式挖掘，获得了12种最大频繁正、负序列模式；

最大频繁序列模式。给定一个DNA序列S，该序列为碱基序列，S＝<s₁ s₂...s_n>，其中s_i(1≤i≤n) 是一个字符的字符集Ω＝{A、T、C、G}。如果一个模式<s_k s_k+1...s_m>(1≤k≤m≤n)的支持度不小于最小支持，这个序列就是频繁序列。最大频繁模式是指它的超序列都不频繁的模式。设 min_sup＝0.3，获得多种最大频繁序列模式。选取其中12种频繁序列模式来作为序列模式分析的数据集。这12种频繁序列模式如下表2：

表2.

实施例3

根据实施例1所述的一种基于生物序列的负序列模式的相似性分析方法，其区别在于：步骤(3) 中，对最大频繁正、负序列模式进行图形表示，包括：在复平面构造一个嘌呤嘧啶图，嘌呤嘧啶图中，第一、二象限是嘌呤，包括A、

G和

第三、四象限是嘧啶，包括T、

C和

(b+di)→A(Ⅰ)

(d+bi)→G(Ⅱ)

(b-di)→T(Ⅲ)

(d-bi)→C(Ⅳ)

四个核苷酸A、G、T、C及其对应的负序列的单位向量

如式(Ⅰ)至式(Ⅷ) 所示：

式(Ⅰ)至式(Ⅷ)中，b和d是非零的实数，

A和T是共轭的，G和C也是共轭的，即，

A、T、C，G代表现实存在的碱基对，

表示的是DNA序列中本该出现却没有出现的碱基对，又称缺失的碱基对，也叫A、G、T、C及其对应的负序列的单位向量；如图2所示。

通过这种表示方法，将一个DNA序列碱基

还原为一个数字序列s(n)，如式(Ⅸ)所示：

式(Ⅸ)中，s(0)＝0，其中y(j)满足式(Ⅹ)：

式(Ⅹ)中，j表示序列S中第0,1,2,...,n个位置上的碱基类型，n是被研究的DNA序列的长度；

通过上面的步骤，从“嘌呤嘧啶图”中唯一地获得原始的DNA序列的时间序列；

利用式(Ⅹ)把12种最大频繁正、负序列模式转化为数字序列；比如序列Human1通过式(Ⅸ) -(Ⅹ)得到的复数字序列为s(H1)＝{0.866+0.5i,1.366-0.366i,2.2321+0.134i,3.0981+0.634i,3.5981+1.5i, 4.4641+2i}，模组成的时间序列为S(H1)＝{1.0000,1.4142,2.2361,3.1623,3.8982,4.8916}。通过这样的方法可以得到12种频繁序列模式转化之后的时间序列。

实施例4

根据实施例1所述的一种基于生物序列的负序列模式的相似性分析方法，其区别在于：

步骤(4)中，通过DTW算法求取距离矩阵，距离矩阵用于表示不同DNA序列的相似度。

设转化DNA序列而获得的时间序列为，

其长度分别为m和n；按照它们的时间位置进行排序，构造m×n矩阵A_m×n，矩阵中的每个元素

在矩阵中，把一组相邻的矩阵元素的集合称为弯曲路径，记为 W＝w₁,w₂,...,w_k，W的第k个元素w_k＝(a_ij)_k，这条路径满足下列条件：

①max{m,n}≤K≤m+m-1；

②w₁＝a₁₁,w_k＝a_mn；

③对w_k＝a_ij,w_k-1＝a_i'j'必须满足0≤i-i'≤1,0≤j-j'≤1，则

DTW算法为运用动态规划思想寻找一条具有最小弯曲代价的最佳路径，如式(Ⅺ)所示：

其中，i＝2,3,...,m；j＝2,3,...,n。D(m,n)为A_m×n中弯曲路径的最小累加值。

通过对12种频繁序列转化之后的时间序列进行DTW距离测量，分别获得8种PSPs和4种NSP之间的距离矩阵，分别如表3、表4所示：

表3

表4

据了解，Human、Chimpanzee属于灵长类动物，Rat属于啮齿类动物，Opossum属于后兽动物。本发明方法的整体变化与其分类一致，因此，本发明提出的方法是有效可行的。并且提出的方法对于短序列和长序列都是有效的，由于本发明使用的数据是挖掘之后的频繁模式，用于比较的序列的长度普遍减短，且保留了原序列的特性，所以计算非常简单，节省了计算机的内存消耗。通过四个物种之间的相似性比较，可以知道不同的模式组合得到不同的结果，这些结果可能在不同的考虑下有用。

随机选取某几个最大频繁序列，序列的距离矩阵(如表3、表4所示)，表3、4中列出的不同数据组的相似度，如果能够合理进行聚类的话，利用本发明方法构建系统发育树。分子进化遗传学分析MEGA5是一个用户友好的软件，用于建立序列比对和系统发育树。系统发生树是一个树状分支图，总结了各种生物的遗传或进化关系。图5(a)为对最大频繁序列Human1,Opossum2,Rat2 and Chimpanzee2进行相似性分析后绘制的系统发育树示意图；图5(a)为对最大频繁序列Human2, Opossum1,Rat2,and Chimpanzee1进行相似性分析后绘制的系统发育树示意图；图6(a)为对最大频繁序列Human2,Opossum2,Rat2 and Chimpanzee1进行相似性分析后绘制的系统发育树示意图；图6(a)为对最大频繁序列Human3,Opossu3,Rat3 and Chimpanzee3进行相似性分析后绘制的系统发育树示意图；本发明选取四种频繁模式的组合，便得到了四种不同的分类结果，这都符合物种的进化规律。

通过归一处理数据，让本发明的结果和其他方法进行比较。图7为归一化物种距离示意图。其中，纵坐标为归一化距离。图7展示了本方法与两种比较方法的结果及MEGA结果之间的Pearson相关系数。表5详细说明了四种方法与其他物种和人类之间的距离。

表5

表5中，括号中的值是归一化到0到1之间的真实距离。Ref.【1】参见ZhiyiMo,WenZhu,Yi Sun,Qilin Xiang,MingZheng,MinChen,ZejunLi.One novel representationof DNA sequence based on the global and local position information.[J].Scientific reports,2018,8(1).Ref.【2】参见Yu Hong-Jie,Huang De-Shuang.Graphical representation for DNA sequences via joint diagonalizationof matrix pencil.[J]. IEEE Journal of Biomedical&Health Informatics,2013,17(3):503-511.计算了本方法与两种比较方法结果之间的Pearson相关系数。

可以看出，本发明方法与MEGA的相关系数最高，说明本发明方法能够更准确地计算出DNA 序列之间的相似性。此外，从图7可以看出，本发明方法与MEGA计算的曲线更加接近，这再次说明本发明方法与MEGA的相关性最高。

比对可知，通过这种方法，可以有效地对负序列进行有效的表达和分析，并且通过选取不同的最大频繁模式组合，能够得到不同的分析结果。由于所选取的是频繁模式来进行的相似性分析，大大节省了计算机的内存和时间的消耗。此方法也和MEGA具有最高的相关性。

实施例5

根据实施例1-4任一所述的一种基于生物序列的负序列模式的相似性分析方法的实现系统，如图3所示，包括依次连接的数据预处理模块、频繁模式挖掘模块、图形表示模块、相似性分析模块；数据预处理模块用于执行步骤(1)；频繁模式挖掘模块用于执行步骤(2)；图形表示模块用于执行步骤(3)；相似性分析模块用于执行步骤(4)。

实施例6

一种计算机可读存储介质，其特征在于，计算机可读存储介质中存储有基于生物序列的负序列模式的相似性分析程序，基于生物序列的负序列模式的相似性分析程序被处理器执行时，实现实施例1-4任一所述的基于生物序列的负序列模式的相似性分析方法的步骤。

Claims (4)

1.一种基于生物序列的负序列模式的相似性分析方法，其特征在于，包括步骤如下：

(1)数据预处理

将DNA序列中的字母用数字来表示；并将数字表示后的DNA序列分割成若干个块，每个块碱基数量相同，得到的若干个块作为频繁模式挖掘的数据集；

(2)频繁模式挖掘

使用f-NSP算法来挖掘数据集，得到最大频繁正、负序列模式；

(3)对最大频繁正、负序列模式进行图形表示；包括：在复平面构造一个嘌呤嘧啶图，嘌呤嘧啶图中，第一、二象限是嘌呤，包括A、

G和

第三、四象限是嘧啶，包括T、

C和

四个核苷酸A、G、T、C及其对应的负序列的单位向量

如式(Ⅰ)至式(Ⅷ)所示：

(b+di)→A(Ⅰ)

(d+bi)→G(Ⅱ)

(b-di)→T(Ⅲ)

(d-bi)→C(IV)

式(Ⅰ)至式(Ⅷ)中，b和d是非零的实数，

A和T是共轭的，G和C也是共轭的，即，

A、T、C，G代表现实存在的碱基对，

表示的是DNA序列中本该出现却没有出现的碱基对，又称缺失的碱基对，也叫A、G、T、C及其对应的负序列的单位向量；

通过这种表示方法，将一个DNA序列碱基

还原为一个数字序列s(t)，如式(Ⅸ)所示：

式(Ⅸ)中，1≤t≤n，s(0)＝0，其中y(j)满足式(X)：

式(X)中，j表示序列S中第0,1,2,...,n个位置上的碱基类型，n是被研究的DNA序列的长度；

利用式(X)把12种最大频繁正、负序列模式转化为数字序列；

(4)DNA序列的相似性分析

求取不同DNA序列的相似度，相似度越小，DNA序列越相似；

求取距离矩阵，距离矩阵用于表示不同DNA序列的相似度；

通过DTW算法求取距离矩阵，设转化DNA序列而获得的时间序列为，

其长度分别为m和n；按照它们的时间位置进行排序，构造m×n矩阵A_m×n，矩阵中的每个元素

在矩阵中，把一组相邻的矩阵元素的集合称为弯曲路径，记为W＝w₁,w₂,...,w_K，W的第k个元素w_k＝(a_ij)_k，这条路径满足下列条件：

①max{m,n}≤K≤m+m-1；

②w₁＝a₁₁,w_K＝a_mn；

③对w_k＝(a_ij)_k,w_k-1＝(a_i'j')_k-1必须满足0≤i-i'≤1,0≤j-j'≤1，则

DTW算法为运用动态规划思想寻找一条具有最小弯曲代价的最佳路径，如式(Ⅺ)所示：

式(Ⅺ)中，i＝2,3,...,m；j＝2,3,...,n，D(m,n)为A_m×n中弯曲路径的最小累加值。

2.根据权利要求1所述的一种基于生物序列的负序列模式的相似性分析方法，其特征在于，步骤(2)中，使用f-NSP算法来挖掘数据集，数据集为D，包括步骤如下：

A、使用GSP算法得到所有的正频繁序列，并将每条正频繁序列对应的位图存入哈希表中；包括：

a、扫描数据集得到所有长度为1的序列模式放入原始种子集P₁中；

b、从原始种子集P₁中获取长度为1的序列模式，并将它们通过连接操作生成长度为2的候选序列集C₂；使用Apriori性质对候选序列集C₂进行剪枝，再通过扫描候选序列集C₂确定其中剩余序列的支持度，将支持度高于最小支持度的序列模式保存下来，输出为长度为2的序列模式L₂，并作为长度为2的种子集；按照该方法一直输出长度为3的序列模式L₃、长度为4的序列模式L₄……长度为n+1的序列模式L_n+1，直到挖掘不出新的序列模式，得到序列模式即所有的正频繁序列，最小支持度是人为设置的支持度阈值min_sup；

B、基于所有正频繁序列生成相应的NSC；

NSC是指负候选序列，正频繁序列统称为正序列，对于一个k-size的PSP，NSCs是通过将任意m个不相邻的元素改变为它的负数来生成的，用

来表示，

是不小于k/2的最小整数；k-size是指序列的大小为k；NSCs是指所有的负候选序列；

C、利用位操作快速计算负侯选序列的支持度；

NSCs的支持度计算如下：给定一个m-size和n-neg-size负序列ns，对于

那么在数据集D中ns的支持度为：

m-size是指序列大小为m；假设ns＝<a₁a₂…a_m>是一个负序列，若ns′仅由ns中所有的正元素组成，则将ns′称作ns的最大正子序列，定义为MPS(ns)；由这个序列的MPS(ns)和ns中一个负元素a组成的序列称为1-neg-size最大子序列，定义为1-negMS；

通过频繁模式挖掘，获得了12种最大频繁正、负序列模式。

3.权利要求1或2所述一种基于生物序列的负序列模式的相似性分析方法的实现系统，其特征在于，包括依次连接的数据预处理模块、频繁模式挖掘模块、图形表示模块、相似性分析模块；所述数据预处理模块用于执行步骤(1)；所述频繁模式挖掘模块用于执行步骤(2)；所述图形表示模块用于执行步骤(3)；所述相似性分析模块用于执行步骤(4)。

4.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质中存储有基于生物序列的负序列模式的相似性分析程序，所述基于生物序列的负序列模式的相似性分析程序被处理器执行时，实现权利要求1或2所述基于生物序列的负序列模式的相似性分析方法的步骤。

Images