CN103218544B - 基于序列相似性与频谱3-周期性的基因识别方法 - Google Patents

Abstract

本发明公开一种基于序列相似性与频谱3-周期性的基因识别方法，属于生物信息学中的基因识别领域。所述基因识别方法利用核苷酸频率特征矩阵来描述DNA序列；对不同序列通过计算核苷酸频率特征矩阵相似度来确定DNA序列的相似度，并根据相似度选定信噪比阈值；通过计算序列核苷酸频谱3-周期性强度以及序列中每个核苷酸位置处的斜率来区分外显子与内含子，并通过修正处理最终确定外显子段。本发明解决了对于给定DNA序列，根据频谱3-周期性基因识别算法中信噪比阈值过于单一，阈值确定过程缺乏考虑基因类别的问题。

Description

基于序列相似性与频谱3-周期性的基因识别方法

技术领域

本发明涉及DNA序列相似性计算、频谱3-周期性强度计算以及根据这两种性质所实现的基因识别方法，属于生物信息学中的基因识别技术领域。

背景技术

DNA是生物遗传信息的载体，是一种长链聚合物，由腺嘌呤（Adenine,A），鸟嘌呤（Guanine,G），胞嘧啶（Cytosine,C），胸腺嘧啶（Thymine,T）这四种核苷酸（nucleotide）符号按一定的顺序连接而成。其中带有遗传讯息的DNA片段称为基因（Gene）（见图1第一行）。在真核生物的DNA序列中，基因通常被划分为许多间隔的片段（见图1第二行），其中编码蛋白质的部分，即编码序列（CodingSequence）片段，称为外显子（Exon），不编码的部分称为内含子（Intron）。外显子在DNA序列剪接（Splicing）后仍然会被保存下来，并可在蛋白质合成过程中被转录（transcription）、复制（replication）而合成为蛋白质。DNA序列通过遗传编码来储存信息，指导蛋白质的合成，把遗传信息准确无误地传递到蛋白质（protein）上去并实现各种生命功能。

对大量、复杂的基因序列的分析，传统生物学解决问题的方式是基于分子实验的方法，其代价高昂。诺贝尔奖获得者W.吉尔伯特（WalterGilbert，1932—，第一个制备出混合脱氧核糖核酸的科学家）1991年曾经指出：“现在，基于全部基因序列都将知晓，并以电子可操作的方式驻留在数据库中，新的生物学研究模式的出发点应是理论的。一个科学家将从理论推测出发，然后再回到实验中去，追踪或验证这些理论假设。”随着世界人类基因组工程计划的顺利完成，通过物理或数学的方法从大量的DNA序列中获取丰富的生物信息，对生物学、医学、药学等诸多方面都具有重要的理论意义和实际价值，也是目前生物信息学领域的一个研究热点。

基因预测，是一个尚未完全解决的问题，也是当前生物信息学的一个最基础、最首要的问题。对于预测，一类方法是基于统计学的，这类方法在对基因信息了解不多的情况下，基因识别的准确率会明显下降。因此在目前基因预测研究中，采用信号处理与分析方法来发现基因编码序列也受到广泛重视。对于基因序列首先要进行数字映射，将其转换为指示序列。然后，对指示序列做离散傅里叶变换，可以计算平方功率谱，求信噪比等。

科学家发现，基因序列存在频谱3-周期特性，这个特性是区分外显子和内含子的重要特征。同时，也有人认为频谱3-周期特性的存在是由于“密码子”使用的偏向性导致的。目前，有一些通过频谱3-周期特性进行基因识别的算法。

但是，在目前的频谱3-周期特性进行基因识别的方法中，所选取的特定基因类型的DNA序列，将其信噪比R的判别阈值取为R₀=2，带有一定的主观性、经验性。对不同的基因类型，所选取的判别阈值应该是不同的。否则，根据某个固定阈值去判断，有可能将本来是外显子的片段误认为是内含子，造成基因识别的准确度下降。因此，我们希望对待测序列先做相似性分析，以与之最相近的基因信噪比阈值作为该序列的判别阈值，然后计算每一位核苷酸的频谱3-周期性强度及斜率，根据所选判别阈值识别出外显子，从而最大程度保证基因识别的精度。

发明内容

发明目的针对现有利用频谱3-周期性进行基因识别的方法中阈值选取的不足，本发明提出了一种基于序列相似性与频谱3-周期性的基因识别方法，通过利用核苷酸频率特征矩阵来描述DNA序列，对不同序列通过计算核苷酸频率特征矩阵相似度来确定DNA序列的相似度，并根据相似度选定信噪比阈值；计算序列核苷酸频谱3-周期性强度以及序列中每个核苷酸位置处的斜率来区分外显子与内含子，并修正处理最终确定外显子段，从而提高基因识别的准确度。

技术方案：一种基于序列相似性与频谱3-周期性的基因识别方法，主要包括DNA特征值计算和基因识别算法两部分，其中：所述DNA特征值计算包括构建核苷酸频率特征矩阵、计算DNA序列相似度、频谱3-周期性强度计算及每一位核苷酸的斜率计算，所述基因识别算法包括利用序列相似度确定信噪比阈值、判断每一位核苷酸是否属于外显子以及最后对外显子段进行修正；

步骤1，对待识别DNA序列构建其核苷酸频率特征矩阵，根据DNA序列相似度计算公式计算该DNA序列与已知物种序列的相似度，以相似度最近已知物种的信噪比阈值作为所述DNA序列识别时所用的信噪比阈值R₀；

步骤2，对待测DNA序列的每一位核苷酸计算其频谱3-周期性强度以及该核苷酸位置处的斜率，用步骤1所得的信噪比阈值R₀来区分外显子与内含子；

步骤3，对外显子段进行修正，对于步骤2识别的结果，若外显子长度小于预设值，则应将其归为内含子；若内含子长度小于预设值，则应归为外显子；最终得到基因识别结果。

所述核苷酸频率特征矩阵是三个四阶矩阵，如下所示：

$M_{\mathrm{sx}}=\left(s_{\mathrm{xij}}\right),s_{\mathrm{xij}}=\left\{\begin{array}{cc}1,& f_{\mathrm{xi}}>f_{\mathrm{xj}}\\ 0,& f_{\mathrm{xi}}=f_{\mathrm{xj}}\\ -1& f_{\mathrm{xi}}<f_{\mathrm{xj}}\end{array}\right.i=\mathrm{1,2,3,4},j=\mathrm{1,2,3,4}$ 公式1.1

$M_{\mathrm{sy}}=\left(s_{\mathrm{yij}}\right),s_{\mathrm{yij}}=\left\{\begin{array}{cc}1,& f_{\mathrm{yi}}>f_{\mathrm{yj}}\\ 0,& f_{\mathrm{yi}}=f_{\mathrm{yj}}\\ -1& f_{\mathrm{yi}}<f_{\mathrm{yj}}\end{array}\right.i=\mathrm{1,2,3,4},j=\mathrm{1,2,3,4}$ 公式1.2

$M_{\mathrm{sz}}=\left(s_{\mathrm{zij}}\right),s_{\mathrm{zij}}=\left\{\begin{array}{cc}1,& f_{\mathrm{zi}}>f_{\mathrm{zj}}\\ 0,& f_{\mathrm{zi}}=f_{\mathrm{zj}}\\ -1& f_{\mathrm{zi}}<f_{\mathrm{zj}}\end{array}\right.i=\mathrm{1,2,3,4},j=\mathrm{1,2,3,4}$ 公式1.3

其中，设DNA序列S的长度为N，核苷酸编号从0到N-1，公式中M_sx、M_sy、M_sz分别表示DNA序列S的在x、y、z位置上出现核苷酸频率的特征值，f_xi表示第i种核苷酸出现在0，3，6，…，N-3位置处的频率，f_xj表示第j种核苷酸出现在0，3，6，…，N-3位置处的频率；f_yi表示第i种核苷酸出现在1，2，7，…，N-2位置处的频率，f_yj表示第j种核苷酸出现在1，2，7，…，N-2位置处的频率；f_zi表示第i种核苷酸出现在2，5，8，…，N-1位置处的频率，f_zj表示第j种核苷酸出现在2，5，8，…，N-1位置处的频率。这里，核苷酸的顺序取A，T，G，C，即第一种核苷酸为A，第二种为T，以此类推；

所述DNA序列相似度的计算定义为计算序列所对应的核苷酸频率特征矩阵的相似度。对于两个DNA序列S₁与S₂，分别对其求得核苷酸频率特征矩阵M_sx ¹、M_sy ¹、M_sz ¹、M_sx ²、M_sy ²、M_sz ²。对两个基因序列在x，y，z位置上的相似性d(S_x ¹,S_x ²)，d(S_y ¹,S_y ²)，d(S_z ¹,S_z ²)定义如下：

$\left\{\begin{array}{c}d({s_x}^1,{s_x}^2)=\frac{1}{2*4(4-1)}{\left|\right|{M_{\mathrm{sx}}}^1-{M_{\mathrm{sx}}}^2\left|\right|}_F=\frac{1}{24}\sqrt{\underset{i,j=1}{\overset{4}{\mathrm{\&Sigma;}}}{|{s_{\mathrm{xij}}}^1-{s_{\mathrm{xij}}}^2|}^2}\\ d({s_y}^1,{s_y}^2)=\frac{1}{2*4(4-1)}{\left|\right|{M_{\mathrm{sy}}}^1-{M_{\mathrm{sy}}}^2\left|\right|}_F=\frac{1}{24}\sqrt{\underset{i,j=1}{\overset{4}{\mathrm{\&Sigma;}}}{|{s_{\mathrm{yij}}}^1-{s_{\mathrm{yij}}}^2|}^2}\\ d({s_z}^1,{s_z}^2)=\frac{1}{2*4(4-1)}{\left|\right|{M_{\mathrm{sz}}}^1-{M_{\mathrm{sz}}}^2\left|\right|}_F=\frac{1}{24}\sqrt{\underset{i,j=1}{\overset{4}{\mathrm{\&Sigma;}}}{|{s_{\mathrm{zij}}}^1-{s_{\mathrm{zij}}}^2|}^2}\end{array}\right.$ 公式2

两个DNA序列的相似性d(S¹,S²)定义为：

$d(s^1,s^2)=\sqrt{d{({s_x}^1,{s_x}^2)}^2+d{({s_y}^1,{s_y}^2)}^2+d{({s_z}^1,{s_z}^2)}^2}$ 公式3

d(s¹,s²)越小，就认为两个序列S₁、S₂越相似，否则越不相似；

所述频谱3-周期性强度计算需要计算DNA序列在位置处的总功率谱值，然后计算DNA序列背景噪声的三周期信号比。在长度为N的DNA序列S=S₀S₁S₂…S_n-1中，将核苷酸符号b∈I=﹛A,T,G,C﹜出现在该序列的0,3,6,...N－3与1,4,7,…N－2以及2,5,8,…N－1等位置上的频数分别记为x_b，y_b和z_b，则得到序列在处的总功率谱值如下：

$P\left[\frac{N}{3}\right]=\underset{b\&Element;I}{\mathrm{\&Sigma;}}(x_b^2+y_b^2+z_b^2-x_b{y}_b-x_b{z}_b-y_b{z}_b),I=\{A,T,G,C\}$ 公式4

DNA序列背景噪声的三周期信号比定义为：

$S\left[N\right]=\frac{P\left[\frac{N}{3}\right]}{N}$ 公式5

S[N]就作为每个核苷酸的频谱3-周期性强度；

所述每一位核苷酸的斜率计算如下：

$\mathrm{slope}\left(i\right)=\frac{s\left[i\right]-s[i-50]}{50},i=\mathrm{50,51}...N-1$ 公式6

其中i表示DNA序列第i个核苷酸，由于大部分外显子和内含子的长度大于50，故间隔50位计算。

所述外显子段修正为：如果我们判断出一段区域是外显子，且其长度小于50，则应对其修正，将其归为内含子；同样的，如果我们判断出一段区域是内含子，且其长度小于50，则应对其修正，将其归为外显子。

本发明采用上述技术方案，具有以下有益效果：本发明解决了对于给定DNA序列，根据频谱3-周期性基因识别算法中信噪比阈值过于单一，阈值确定过程缺乏考虑基因类别的问题，可以提高基因识别的准确度。

附图说明

图1为现有技术中真核生物DNA序列（基因序列）结构示意图；

图2为本发明实施例的方法流程图；

图3为本发明实施例的中外显子识别的流程图；

图4为本发明实施例的外显子内含子修正示意图。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图2所示，基于序列相似性与频谱3-周期性的基因识别方法分为三个步骤：

步骤1，对待识别DNA序列构建其核苷酸频率特征矩阵，根据相似度计算公式计算该DNA序列与已知物种序列的相似度，以相似度最近物种的信噪比阈值作为该序列识别时所用的信噪比阈值R₀。

假设给定一段DNA序列片段S=AGCGTACCGCTA，则各核苷酸出现频率为：

$\begin{array}{cccc}{f}_{x1}=\frac{1}{4}& f_{x2}=0& f_{x3}=\frac{1}{4}& f_{x4}=\frac{1}{2}\\ f_{y1}=0& f_{y2}=\frac{1}{2}& f_{y3}=\frac{1}{4}& f_{y4}=\frac{1}{4}\\ f_{z1}=\frac{1}{2}& f_{z2}=0& f_{z3}=\frac{1}{4}& f_{z4}=\frac{1}{4}\end{array}$

由公式1.1、1.2、1.3计算出M_sx、M_sy、M_sz三个矩阵为：

$M_{\mathrm{sx}}=\left[\begin{array}{cccc}0& 1& 0& -1\\ -1& 0& -1& -1\\ 0& 1& 0& -1\\ 1& 1& 1& 0\end{array}\right]$ $M_{\mathrm{sy}}=\left[\begin{array}{cccc}0& -1& -1& -1\\ 1& 0& 1& 1\\ 1& -1& 0& 0\\ 1& -1& 0& 0\end{array}\right]$ $M_{\mathrm{sz}}=\left[\begin{array}{cccc}0& 1& 1& 1\\ -1& 0& -1& -1\\ -1& 1& 0& 0\\ -1& 1& 0& 0\end{array}\right]$

对于M_sx这个矩阵，若某一行-1很多，则这一行所对应的核苷酸在0，3，6，…，N-3位置处出现的频率较低，反之，则较高；对于M_sy这个矩阵，若某一行-1很多，则这一行所对应的核苷酸在1，4，7，…，N-2位置处出现的频率较低，反之，则较高；对于M_sz这个矩阵，若某一行-1很多，则这一行所对应的核苷酸在2，5，8，…，N-1位置处出现的频率较低，反之，则较高。

我们取10条已知基因序列与未知序列做相似性计算，得到结果如表1：

表1序列相似比较及阈值选取

通过对比，我们发现第9号（No.9）已知基因与待测基因相似度最高，故选取其信噪比阈值作为基因识别中的阈值R₀=1.3。

步骤2，对待测DNA序列的每一位计算其频谱3-周期性强度以及该核苷酸位置处的斜率，用步骤1所得的信噪比阈值R₀来区分外显子与内含子，其算法流程如图3所示，具体步骤如下；

步骤2-1，初始化位置游标i=1；

步骤2-2，对DNA序列第i个位置处的核苷酸计算其在位置处的总功率谱值以及背景噪声的三周期信号比S[i]（此时N=i）；

步骤2-3，计算第i个位置处的核苷酸的斜率slope[i]；

步骤2-4，若S[i]≥R₀，且slope[i]＞0成立，则该位置处的核苷酸属于外显子；否则属于内含子；

步骤2-5，令位置游标i=i+1；

若i未超过DNA序列长度，则转向步骤2-2；否则结束步骤2。

对于序列前100位，我们所得背景噪声的三周期信号比、斜率、外显子识别情况如表2：

表2外显子识别情况

步骤3，对步骤2识别的结果进行处理，若外显子长度小于50，则应将其归为内含子；若内含子长度小于50，则应归为外显子。最终得到基因识别结果。

如图4所示，第一条DNA识别结果是经过步骤2得到的识别结果序列。序列中存在部分外显子和内含子段长度小于50，需要进行修正。第一段长度小于50的序列是第145到178位，这一段原来为内含子，现在将其改为外显子，得到序列如图2中的第二条DNA识别结果；其中的第263到281位长度小于50，这一段原来为外显子，现在改为内含子，得到序列如图2中的第三条DNA识别结果；其中的第435到459位长度小于50，这一段原来为外显子，现在改为内含子，得到序列如图2中的第四条DNA识别结果。此时，已经没有外显子或内含子段长度小于50，修正结束。

Claims (4)

1.一种基于序列相似性与频谱3-周期性的基因识别方法，其特征在于，包括DNA特征值的求解和基因识别两部分，其中：所述DNA特征值的求解包括构建核苷酸频率特征矩阵、计算DNA序列相似度、频谱3-周期性强度计算及每一位核苷酸的斜率计算；所述基因识别包括利用序列相似度确定信噪比阈值、判断每一位核苷酸是否属于外显子以及对外显子段进行修正；

步骤1，对待识别DNA序列构建其核苷酸频率特征矩阵，根据DNA序列相似度计算公式计算该DNA序列与已知物种序列的相似度，以相似度最近已知物种的信噪比阈值作为所述DNA序列识别时所用的信噪比阈值R0；

步骤2，对待测DNA序列的每一位核苷酸计算其频谱3-周期性强度以及该核苷酸位置处的斜率，用步骤1所得的信噪比阈值R₀来区分外显子与内含子；

步骤3，对外显子段进行修正，对于步骤2识别的结果，若外显子长度小于预设值，则应将其归为内含子；若内含子长度小于预设值，则应归为外显子；最终得到基因识别结果；

所述核苷酸频率特征矩阵是三个四阶矩阵，如下所示：

M_sx＝(s_xij), $s_{xij}=\left\{\begin{array}{ccc}\begin{array}{cc}1,& f_{xi}>f_{xj}\\ 0,& f_{xi}=f_{xj}\\ -1& f_{xi}<f_{xj}\end{array}& i=1,2,3,4& j=1,2,3,4\end{array}\right.$ 公式1.1

M_sy＝(s_yij), $s_{yij}=\left\{\begin{array}{ccc}\begin{array}{cc}1,& f_{yi}>f_{yj}\\ 0,& f_{yi}=f_{yj}\\ -1& f_{yi}<f_{yj}\end{array}& i=1,2,3,4& j=1,2,3,4\end{array}\right.$ 公式1.2

M_sz＝(s_zij), $s_{zij}=\left\{\begin{array}{ccc}\begin{array}{cc}1,& f_{zi}>f_{zj}\\ 0,& f_{zi}=f_{zj}\\ -1& f_{zi}<f_{zj}\end{array}& i=1,2,3,4& j=1,2,3,4\end{array}\right.$ 公式1.3

其中，设DNA序列S的长度为N，核苷酸编号从0到N-1，公式中M_sx、M_sy、M_sz分别表示DNA序列S的在x、y、z位置上出现核苷酸频率的特征值，f_xi表示第i种核苷酸出现在0，3，6，…，N-3位置处的频率，f_xj表示第j种核苷酸出现在0，3，6，…，N-3位置处的频率；f_yi表示第i种核苷酸出现在1，4，7，…，N-2位置处的频率，f_yj表示第j种核苷酸出现在1，4，7，…，N-2位置处的频率；f_zi表示第i种核苷酸出现在2，5，8，…，N-1位置处的频率，f_zj表示第j种核苷酸出现在2，5，8，…，N-1位置处的频率；这里，核苷酸的顺序取A，T，G，C，即第一种核苷酸为A，第二种为T，以此类推。

2.如权利要求1所述的基于序列相似性与频谱3-周期性的基因识别方法，其特征在于，所述步骤2的具体实施步骤如下：

步骤2-1，初始化核苷酸位置游标i＝1；

步骤2-2，对DNA序列第i个位置处的核苷酸计算其在位置处的总功率谱值以及背景噪声的三周期信号比S[i]，此时N＝i；

步骤2-3，计算第i个位置处的核苷酸的斜率slope[i]；

步骤2-4，若S[i]≥R₀，且slope[i]＞0成立，则该位置处的核苷酸属于外显子；否则属于内含子；

步骤2-5，令位置游标i＝i+1；

若i未超过DNA序列长度，则转向步骤2-2；否则结束步骤2。

3.如权利要求2所述的基于序列相似性与频谱3-周期性的基因识别方法，其特征在于，所述DNA序列相似度的计算定义为计算序列所对应的核苷酸频率特征矩阵的相似度；对于两个DNA序列S₁与S₂，分别对其求得核苷酸频率特征矩阵M_sx ¹、M_sy ¹、M_sz ¹、M_sx ²、M_sy ²、M_sz ²；对两个基因序列在x，y，z位置上的相似性d(S_x ¹,S_x ²)，d(S_y ¹,S_y ²)，d(S_z ¹,S_z ²)定义如下：

$\left\{\begin{array}{c}d\left({s_x}^1,{s_x}^2\right)=\frac{1}{2*4\left(4-1\right)}\left|\right|{M_{sx}}^1-{M_{sx}}^2|{|}_F=\frac{1}{24}\sqrt{\underset{i,j=1}{\overset{4}{\&Sigma;}}|{s_{xij}}^1-{s_{xij}}^2{|}^2}\\ d\left({s_y}^1,{s_y}^2\right)=\frac{1}{2*4\left(4-1\right)}\left|\right|{M_{sy}}^1-{M_{sy}}^2|{|}_F=\frac{1}{24}\sqrt{\underset{i,j=1}{\overset{4}{\&Sigma;}}|{s_{yij}}^1-{s_{yij}}^2{|}^2}\\ d\left({s_z}^1,{s_z}^2\right)=\frac{1}{2*4\left(4-1\right)}\left|\right|{M_{sz}}^1-{M_{sz}}^2|{|}_F=\frac{1}{24}\sqrt{\underset{i,j=1}{\overset{4}{\&Sigma;}}|{s_{zij}}^1-{s_{zij}}^2{|}^2}\end{array}\right.$ 公式2

两个DNA序列的相似性d(S¹,S²)定义为：

$d(s^1,s^2)=\sqrt{d{({s_x}^1,{s_x}^2)}^2+d{({s_y}^1,{s_y}^2)}^2+d{({s_z}^1,{s_z}^2)}^2}$ 公式3

d(s¹,s²)越小，就认为两个序列S₁、S₂越相似，否则越不相似。

4.如权利要求3所述的基于序列相似性与频谱3-周期性的基因识别方法，其特征在于，所述频谱3-周期性强度计算需要计算DNA序列在位置处的总功率谱值，然后计算DNA序列背景噪声的三周期信号比；在长度为N的DNA序列S＝S₀S₁S₂…S_n-1中，将核苷酸符号b∈I＝﹛A,T,G,C﹜出现在该序列的0,3,6,...N－3与1,4,7,…N－2以及2,5,8,…N－1等位置上的频数分别记为x_b，y_b和z_b，则得到序列在处的总功率谱值如下：

$P\&lsqb;\frac{N}{3}\&rsqb;=\underset{b\&Element;I}{\&Sigma;}(x_b^2+y_b^2+z_b^2-x_b{y}_b-x_b{z}_b-y_b{z}_b),I=\{A,T,G,C\}$ 公式4

DNA序列背景噪声的三周期信号比定义为：

$S\&lsqb;N\&rsqb;=\frac{P\&lsqb;\frac{N}{3}\&rsqb;}{N}$ 公式5

S[N]就作为每个核苷酸的频谱3-周期性强度；

所述每一位核苷酸的斜率计算如下：

$slope\left(i\right)=\frac{s\&lsqb;i\&rsqb;-s\&lsqb;i-50\&rsqb;}{50},i=50,51...N-1$ 公式6

其中i表示DNA序列第i个核苷酸；

所述外显子段修正为：如果判断出一段区域是外显子，且其长度小于50，则应对其修正，将其归为内含子；同样的，如果判断出一段区域是内含子，且其长度小于50，则应对其修正，将其归为外显子。