一种用一维细胞神经网络检测DNA序列相似度的方法

摘要

本发明公开了一种用一维细胞神经网络检测DNA序列相似度的方法，首先设计出一维细胞神经网络基本模型，然后利用这种模型建构一个一维的对偶细胞神经网络；再用两个待检测的DNA序列信息对该网络进行初始化，网络运行过程中，记录各时刻网络中的细胞状态和输出，据此形成最优输出矩阵；再对最优输出矩阵中的元素进行遍历，从而确定最佳的对齐路径；最后根据对齐路径对两个序列进行空格插入操作以便将两个序列全局对齐；序列对齐后，再根据对齐的碱基数量和总的碱基数量来计算其全局相似度。经过测试对比表明，本发明在保证检测准确的基础上，对于长度较长的DNA序列，所需的计算时间比现有方法明显有较大幅度地减少。

说明书

技术领域

本发明属于生物信息学中的DNA序列相似度检测技术领域，更为具体地讲，涉及一种用一维细胞神经网络检测DNA序列相似度的方法，用于对DNA双序列全局相似度的检测。

背景技术

20世纪70年代，DNA测序方法的出现产生出许多生物分子序列数据，这些数据正以几何速度迅速增长，它已成为人类实践产生数据量最大的领域。在人类基因组序列图绘制成功后，人们又相继启动了各种动植物的基因组测序计划。但是，数据并不等于知识和信息，研究和处理这些数据的任务越来越重，我们必须寻找高效地方法来解决这类问题。

DNA一般是通过碱基配对相连接以双链形式存在，而碱基的配对存在特异性，总是一条链上的碱基G与另一条链上的碱基C连接，一条链上的碱基T与另一条链上的碱基A连接。DNA核酸序列就是由这4个基本元素组成的字符序列。因此，DNA序列匹配实际上就是匹配两个由ACGT这4个字符中任意一个字符组成的序列之间的相似度。序列比对就是一个通过某种特定的算法寻找两条或多条序列之间最大匹配。匹配碱基数的过程通过序列比对的方法来挖掘序列之间在结构或功能上的相似性，这对于生物数据库的搜索算法，蛋白质或DNA的结构预测、进化分析和功能分析具有非常重要的实践意义。

根据进行比对的生物序列的个数的不同，序列比对方法可以分为双序列比对方法和多序列比对方法。双序列比对方法又可以分为三种，分别是点阵法、动态规划算法和启发式算法（BLAST算法、FASTA算法等）。多序列比对是一个NP完全问题，是一个尚未解决的难题，它可以分为以下几种：精确比对算法、迭代比对算法、渐进比对算法、启发式算法和基于图论的比对算法等。

双序列比对方法中，点阵法是1970年McIntyre和Gibbs首先提出来的，是最基本的一种可视化的双序列比对方法点阵法的优点是可以直接的发现两个序列间所有可能的匹配，但是它得到的比对结果不够精确，而且只适用于较短的两个序列，面对如今数据量庞大的生物序列数据明显存在着缺陷。动态规划算法的基本思想就是将待求解的问题分解成若干个子问题，先分别把子问题的解求解出来，然后存储子问题的解而避免重复计算，最后通过将子问题的解合并起来就得到了原问题的解。采用动态规划算法求解生物序列比对问题可以在给定的得分系统下得到最优的比对结果，但是如果问题量特别大，那么它的计算速度会非常慢，而且这种方法对参数的选择很敏感，参数的微小改动也会使比对的结果有着较大的变化。求解生物序列比对问题的动态规划算法主要有1970年由Needleman和Wunsch提出的一种全局序列比对算法-Needleman-Wunsch算法（简称NW算法），Smith和Waterman于1981年提出的一种用来解决寻找具有局部相似性区域的Smith-Waterman算法（简称为SW算法），1985年由Pearsom和Lipman首先提出并在1988年进行了改进的一种FASTA算法双序列比对的启发式算法，1990年由Altschul等人提出的一种BLAST算法双序列比对的启发式算法。

而传统的比对算法在解决数据量较大的双序列比对问题时，所需要的时间和存储空间随着序列条数和序列长度的增长呈指数级增长，因此，我们需要研究更好更新的方法来提高算法的搜索速度，减少计算时间。

发明内容

本发明的目的在于克服现有技术的不足，提供一种用一维细胞神经网络检测DNA序列相似度的方法，以减少计算时间。

为实现上述发明目的，本发明用一维细胞神经网络检测DNA序列相似度的方法，其特征在于，包括以下步骤：

（1）、设计一维细胞神经网络基本模型

将单细胞进行链状排列，各细胞序号依次用“…、i-1、i、i+1、…”来表示，其中的字母i表示细胞的排列序号；

该基本模型中细胞状态用微分方程组来表示：

$>\left(\begin{array}{c}C\frac{{\partial x}_i\left(t\right)}{\partial t}=-\frac{x_i\left(t\right)}{R_x}+A\otimes Y_i\left(t\right)+B\otimes U_i\left(t\right)+I_i\\ y_i\left(t\right)=f\left(x_i\left(t\right)\right)\end{array}\right)---\left(1\right)$>

其中，方程组（1）中，t表示时间，x_i表示细胞i的状态，A是反馈模板，B是控制模板，I_i、R_x和C分别是三个常量，f(x_i(t))是细胞状态的输出调制函数；Y_i(t)表示细胞i包括自己的邻域输出矩阵，U_i(t)表示细胞i包括自己的邻域输入，分别表示为：

$>\left(\begin{array}{cc}{Y}_i\left(t\right)=\left(\begin{array}{c}{y}_{i-1}\left(t\right)\\ y_i\left(t\right)\\ y_i(t+1)\end{array}\right)& \end{array},U_i\left(t\right)=\left(\begin{array}{c}{u}_{i-1}\left(t\right)\\ u_i\left(t\right)\\ u_{i+1}\left(t\right)\end{array}\right)\right);$>

y_i-1(t)、y_i(t)和y_i+1(t)分别表示细胞i-1、i和i+1的输出，u_i-1、u_i和u_i+1表示细胞i、i-1和i+1接收的细胞输入，“”表示矩阵的卷积运算；

细胞输出调制函数f(x_i(t))的具体形态为：

$>y_i\left(t\right)=f\left(x_i\left(t\right)\right)=\frac{1}{2}\left(\right|x_i\left(t\right)+1|-|x_i\left(t\right)-1\left|\right)---\left(2\right)$>

（2）、构建一维对称细胞神经网络

用步骤（1）设计的一维细胞神经网络模型，先分别生成主子网CNN1和从子网CNN2，再由二者构建一个一维对偶细胞神经网络：

在一维对偶细胞神经网络中，主子网CNN1是固定不动的，而从子网CNN2则是可以沿主子网CNN1平行移动，时间t每增加1，从子网CNN2移动一步，且从子网CNN2每次移动的距离等于主子网CNN1中两个相连细胞之间的距离；主子网CNN1由细胞0、1、2、…、m-1组成，从子网CNN2由细胞0、1、2、…、n组成；

在一维对偶细胞神经网络中，令C=1、R_x=1，则用以表示细胞状态的微分方程简化为：

$>x_i(t+1)=\underset{l\in L\left(i\right)}{\Sigma }A\otimes Y_l\left(t\right)+\underset{l\in L\left(i\right)}{\Sigma }B\otimes U_l\left(t\right)+I_i---\left(3\right)$>

公式(3)中，L(i)表示细胞i在主子网CNN1、从子网CNN2中的细胞邻域即包含细胞i自己、主子网CNN1的前一个细胞i-1、后一个细胞i+1以及细胞i在从子网CNN2中对应的细胞，l则代表细胞i的细胞邻域中的第l个细胞，即l∈L(i)；

时间T=t+1时，细胞i的输出y_i(t+1)相应被重新定义为：

$>\left(\begin{array}{c}{y}_i(t+1)=f\left(x_i(t+1)\right)\\ =\frac{1}{2}\left(\right|x_i(t+1)+1|-|x_i(t+1)-1\left|\right)\end{array}\right)---\left(4\right)$>

（3）、利用步骤（2）构建的一维对偶细胞神经网络，对两个待检测相似度的DNA序列进行全局的碱基对齐；

3.1）、对偶细胞网络的初始化

待匹配的两个DNA碱基序列S₁和S₂的碱基数量分别为K₁和K₂，碱基序列的碱基代码分别表示为S₁(k₁)和S₂(k₂)，且0≤k₁≤K₁-1和0≤k₂≤K₂-1，则主子网CNN1和从子网CNN2的细胞数量分别被初始化为K₁+1和K₂+1，即细胞数量m=K₁+1和n=K₂+1；

用u¹(i)和u²(j)表示主子网CNN1的第i个细胞输入和从子网CNN2的第j个细胞输入，则满足0≤i≤K₁且0≤j≤K₂，主子网CNN1和从子网CNN2中各细胞的细胞输入分别按公式（5）和公式（6）进行赋值：

其中，符号“*”表示细胞的输入u设置为空值；

主子网CNN1中的另一个常量参数初始化赋值为I_i=2；主子网CNN1中使用到的反馈模板Α和控制模板B分别初始化为下列两个常数矩阵：

A=[0 1 0]和B=[0 1 -1]

此外，还要将主子网CNN1中细胞i，i=0,1,..,K₁，的初始状态即t=0时均分别设置为x_i(0)=0、y_i(0)=0；主子网CNN1的第0个细胞和从子网CNN2第K₂个细胞对齐；

3.2）、迭代地计算主子网CNN1中细胞在各时刻的状态和输出

时间t每增加1，从子网CNN2沿主子网CNN1的排列需要增加方向移动一步；

对主子网CNN1，如果细胞i的正下方的那个细胞j_L存在，则选取其3个邻域细胞即细胞i-1,i以及从子网CNN2中正处于i正下方的那个细胞j_L；在时间t，t=1,2,…,m+n-1时，当时间t和细胞序号i同时满足条件1≤t≤m+n-1和1≤i≤m+1时，分别计算各细胞的最优状态和最优输出而如果细胞i的正下方的那个细胞j_L不存在，则不计算细胞最优状态的和最优输出值；

所述的最优状态和最优输出分别按下列公式计算：

$>\overline{x_i}\left(t\right)=\max \{x_{i-1}(t-2)+{2I}_i,x_{i-1}(t-1)-I_i,x_i(t-1)-I_i\}---\left(7\right)$>

其中，函数max(...)表示求取输入参数中的最大值，x_i-1(t-2)、x_i-1(t-1)和x_i(t-1)均按公式(3)进行计算；

3.3）、形成细胞的最优输出矩阵

根据步骤3.2）计算得到主子网CNN1的所有细胞在各时刻的最优状态和最优输出，然后按照第1列为细胞1从t=1到n时刻最优输出、第2列为细胞2从t=2到1+n时刻最优输出、…、第m列为细胞m从t=m到m+n时刻最优输出得到主网络CNN1的最终的细胞最优输出矩阵S_y；

3.4）、对两个DNA序列的碱基进行全局对齐

根据步骤3.3）得到的最优输出矩阵S_y,从矩阵左上角的元素开始，从左至右、从上至下遍历矩阵，确定最优输出矩阵S_y中值为1的矩阵元素位置，并将确定好的各元素按顺序连接以形成碱基的对齐路径P；

根据已确定好的碱基对齐路径P，分三种情况对DNA碱基序列S₁和S₂进行操作：从第一个元素1开始，如果下一个1位于其下，则在序列S₁的当前位置插入符号“*”；如果下一个元素1位于其右侧，则在序列S₂的当前位置插入符号“*”；如果下一个1刚好位于其右下位置，则不对S₂和S₂的当前位置做任何操作。

处理完成S_y的第一个元素后，按前述的三种情况继续处理第二个元素，直到输出矩阵S_y的全部值为1的元素全部都已经被处理完毕，此时S₂和S₂已按碱基的排列顺序完成了全局对齐；

（4）、计算两个DNA碱基序列序列S₂和S₂的全局相似度

定义序列S₂和S₂的全局相似度为SC(S₁,S₂)，则这两个DNA碱基序列的全局相似度按如下公式进行计算：

$>\mathrm{SC}(S_1,S_2)=\frac{2\times N_{\mathrm{match}}}{\mathrm{Len}\left(S_1\right)+\mathrm{Len}\left(S_2\right)}\times 100\%---\left(9\right)$>

其中，符号N_match表示两个DNA碱基序列S₂和S₂经全局序列对齐以后，匹配成功的碱基对数量，Len(S₁)和Len(S₂)分别表示序列S₁和S₂的实际长度。

本发明的发明目的是这样实现的：

本发明用一维细胞神经网络检测DNA序列相似度的方法，首先设计出一维细胞神经网络基本模型，然后利用这种模型建构一个一维的对偶细胞神经网络；再用两个待检测的DNA序列信息对该网络进行初始化，网络运行过程中，记录各时刻网络中的细胞状态和输出，据此形成最优输出矩阵；再对最优输出矩阵中的元素进行遍历，从而确定最佳的对齐路径；最后根据对齐路径对两个序列进行空格插入操作以便将两个序列全局对齐；序列对齐后，再根据对齐的碱基数量和总的碱基数量来计算其全局相似度。经过测试对比表明，本发明在保证检测准确的基础上，对于长度较长的DNA序列，所需的计算时间比现有方法明显有较大幅度地减少。

附图说明

图1是本发明涉及到的一维细胞神经网络基本模型示意图；

图2是图1所示的一维细胞神经网络基本模型中单个细胞的结构图；

图3是本发明构建的一维对偶细胞耦合神经网络示意图；

图4是本发明中全局的碱基对齐流程图；

图5是本发明所述方案在初始状态(t=0)时两个子网的位置关系图；

图6是本发明所述的细胞状态矩阵和输出矩阵的串接规则图；

图7是三种方法的“碱基总数-计算时间”曲线对比图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

在本实施列中，将单细胞进行链状排列，设计的一维细胞神经网络基本模型如图1所示。各细胞序号依次用“…、i-1、i、i+1、…”来表示，其中的字母i表示细胞的排列序号。

一维细胞神经网络基本模型中单个细胞的结构如图2所示，其中用x_i表示细胞i的状态，y_i、y_i-1和y_i+1分别表示从细胞i从其邻域细胞i、i-1和i+1中接收到的反馈输出，u_i、u_i+1和u_i-1表示细胞i从其邻域细胞i、i-1和i+1接收的细胞输入，A是反馈模板，B是控制模板，I_i、R_x和C分别是三个常量，f(x_i)是细胞状态的输出调制函数。

在本实施例中，基于如图1所示的一维细胞神经网络模型，先分别生成主子网CNN1和从子网CNN2，再由二者构建一个如图3所示的一维对偶细胞神经网络。

如图3所示，一维对偶细胞神经网络中，主子网CNN1是固定不动的，而从子网CNN2则是可以沿CNN1平行移动的。当t=1时，从子网CNN2开始移动，t每增加1，从子网CNN2移动一步，且从子网CNN2每次移动一步的距离等于主子网CNN1中两个相连细胞之间的距离。CNN1由细胞0、1、2、…、m-1组成，用i表示细胞序号，则有0≤i≤m-1。CNN2由细胞0、1、2、…、n组成，j表示细胞序号，则有0≤j≤n-1。图3中字符u和y（图中的下标为细胞序号）仍然表示细胞的输入和输出，且实线箭头表示细胞之间存在连接输入，虚箭头表示细胞之间不存在连接输入。

如图3所示，公式(3)中，L(i)表示细胞i在主子网CNN1、从子网CNN2中的细胞邻域即包含细胞i自己、主子网CNN1的前一个细胞i-1、后一个细胞i+1以及细胞i在从子网CNN2中对应的细胞，在图3中即为加黑的四个细胞。其中，l则代表细胞i的细胞邻域中的第l个细胞，即l∈L(i)；

在本实施例中，利用步骤（2）构建的一维对偶细胞神经网络，对两个待检测相似度的DNA序列进行全局的碱基对齐，相应的流程如图4所示。为一个时刻输出生成最优输出矩阵和全局对齐的过程，与发明内容相同，不再赘述。

在对其过程中，对偶细胞网络的初始化时，主子网CNN1和从子网CNN2的关系如图5所示，主子网CNN1的第0个细胞和从子网CNN2第K₂个细胞对齐。

根据步骤3.2）计算得到主子网CNN1的所有细胞在各时刻的状态和输出，按如图6所示所示的顺序连接，得到主子网CNN1的最终的细胞最优输出矩阵S_y。

细胞最优输出矩阵S_y串接序列是：$>\overline{y_1}\left(1\right),\overline{y_1}\left(2\right),...,\overline{y_m}(m+n-2),$>$>\overline{y_{m-1}}(m+n-2),\overline{y_m}(m+n-1).$>

实例

下面结合两个分别取NCBI数据库中真实存在的DNA碱基序列来对本发明的实施过程作进一步的具体描述。

选取的DNA碱基序列识别号分别是S62051和NM_008134。为了展示方便，只选取其中的两个序列片段来进行实施过程说明，这两个序列片段的详细情况如表1所示：

DNA碱基序列碱基数序列碱基片段代码S₁8AAGCTCTGS₂6CAGCAT

表1

如表1所示的两个DNA序列片段S₁和S₂，它们包含的碱基数量分别为K₁=8和K₂=6，分别按步骤3下属的4个子步骤对S₁和S₂进行对齐，分别如下：

①按步骤（3）的步骤3.1），对设计的一维对偶细胞神经网络进行初始化：

m=8+1=9，n=6+1=7，CNN1中C=1,Rx=1,I_i=2。

u¹(0)=*，u¹(1)=S1(0)=A、u¹(2)=S1(1)=A，…,u¹(8)=S1(7)=G；

u²(0)=*，u²(1)=S2(0)=C、u²(2)=S2(1)=A,…,u²(6)=S2(5)=T；

②按步骤（3）的步骤3.2），对偶细胞神经网络迭代运行，根据公式(7)和(8)分别计算当时间t=1，2，3，…，14时，主子网CNN1中各细胞的最优状态和最优输出；

③按步骤（3）的步骤3.3），将各时刻的细胞最优输出串接而成最优输出矩阵S_y，得到的S_y如表2所示。

100000000100000000011000000000100000000010000000010000000001111

表2

④按步骤（3）的步骤3.4），对DNA碱基序列S₁和S₂的碱基进行全局对齐,全局对齐以后的两个碱基序列如表3所示。

*AAGC*TCTGCA*GCAT***

表3

⑤按步骤（4）中的公式(9),计算序列S₁和S₂的全局相似度。根据序列的全局对齐结果(表3所示)，可得到N_match=4，Len(S₁)=8，且Len(S₁)=8，则这两个序列的全局相似度SC(S₁,S₂)=(2×4)÷(8+6)=57.14%

在本实施例中，本发明的方法，还在NCBI数据库中的大量真实DNA序列上进行了实施验证，并分别与现有技术的MILP方法和SPA方法进行了对比。实施过程中对比的主要指标是各方案的计算时间和方案所获得的序列全局相似度，详细的对比情况如表4和表5所示。

表4

表5

其中，表4是本发明与现有技术方案的序列相似度对比表，表5是本发明与现有技术方案的计算时间对比（单位：毫秒）

如表4和表5所示，从左到右各列对应的两个DNA碱基序列长度之和是逐渐增加的，两个表分别展示了来自于NCBI数据库中的6组DNA碱基序列用不同方法计算得到的相似度和所需的计算时间。从表4中显示的相似度数据库可以看出对于部分短序列对（如S62051：长度226与NM_008134：长度625），本发明中介绍的计算方法得到相似度与其它两种方法一致，但随着序列总和长度的增加（如NG_009301:长度42028与NM_000405：长度3690），本发明中的方法获得的相似度比其它两种方法略高。其主要原因就是本发明的方法在对齐DNA序列时，能够比其它两种方法获得更多的对齐碱基对数目。从表5中显示的计算时间数据也可以看出，对于短序列对（如S62051：长度226与NM_008134：长度625），三种方法的计算时间没有显著差别，本发明的方法仅比其它两种方法少用10毫秒左右。但随着序列总和长度的增加（如NG_009301:长度42028与NM_000405：长度3690），本发明所需的计算时间约是SPA方法的33%，也仅是MILP方法所需时间的45%左右。

为了展示计算时间与序列长度和之间的变化规律，分别绘制了这三种方法各自的“碱基总数-计算时间”曲线。图7中展示的“碱基总数-计算时间”对比曲线表明，当两个序列的长度总和较小时（如小于5000），三条曲线基本重合，意味着此时三种方法所需的计算时间并无显著差异。当两个DNA序列的长度总和继续增加时，SPA和MILP方法对应的时间曲线会陡峭爬升，与此同时本发明的时间曲线也会爬升，但其爬升的陡峭程度比另两个曲线平缓很多。由此可见，当两个序列的总长度和较大时，本发明所需的计算时间比SPA和MILP方法明显有较大幅度地减少。

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