CN103150491A - 基于核苷酸位差的频谱3-周期性信噪比计算方法 - Google Patents

Abstract

本发明公开一种基于核苷酸位差的频谱3-周期性信噪比计算方法，属于生物信息学中的基因识别领域。所述信噪比计算方法首先要将DNA序列映射为4个指示序列，然后分别统计并记录4个指示序列里任意两个“1”出现位置之差，并把这些不同的距离值所出现的次数累积存储到一维数组中，再利用余弦函数的周期性，根据诱导公式，最终计算出该DNA序列的信噪比。本发明解决了对于给定的DNA序列，计算其信噪比效率不高的问题。

Description

基于核苷酸位差的频谱3-周期性信噪比计算方法

技术领域

本发明涉及一种基于核苷酸位差的频谱3-周期性信噪比计算方法，包含DNA序列频谱计算、信噪比计算，属于生物信息学中的基因识别技术领域。

背景技术

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

对给定的DNA序列，怎么去识别出其中的编码序列，也称为基因预测，是一个尚未完全解决的问题，也是当前生物信息学的一个最基础、最首要的问题。

基因预测问题的一类方法是基于统计学的。很多国际生物数据网站上也有“基因识别”的算法。比如知名的数据网站http://genes.mit.edu/GENSCAN.html提供的基因识别软件GENSCAN（由斯坦福大学研究人员研发的、可免费使用的基因预测软件），主要就是基于隐马尔科夫链（HMM）方法。但是，它预测人的基因组中有45000个基因，相当于现在普遍认可数目的两倍。另外，统计预测方法通常需要将编码序列信息已知的DNA序列作为训练数据集来确定模型中的参数，从而提高模型的预测水平。但在对基因信息了解不多的情况下，基因识别的准确率会明显下降。

因此在目前基因预测研究中，采用信号处理与分析方法来发现基因编码序列也受到广泛重视。图2是基于序列频谱3-周期性的基因识别方法流程图：在DNA序列研究中，首先需要把A、T、G、C四种核苷酸的符号序列，根据一定的规则映射成相应的数值序列，以便于对其作数字处理。

令I＝{A,T,G,C}，长度（即核苷酸符号个数）为N的任意DNA序列，可表达为

S＝{S[n]|S[n]∈I,n＝0,1,2,…N-1}

即A、T、G、C的符号序列S：S[0],S[1],…,S[N-1]。现对于任意确定的b∈I，令

$u_b\left[n\right]=\left\{\begin{array}{cc}1,& S\left[n\right]=b\\ 0,& S\left[n\right]\≠b\end{array},n=\mathrm{0,1,2},...N-1\right.$

称之为Voss映射，于是生成相应的0-1序列（即二进制序列）{u_b[n]}：

u_b[0],u_b[1],…,，

u_b[N-1](b∈I)，称为指示序列。对指示序列分别做离散Fourier变换(DFT)

$U_b\left[k\right]=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{u}_b\left[n\right]e^{-j\frac{2\mathrm{\πnk}}{N}},k=\mathrm{0,1},...,N-1$   公式1

以此可得到四个长度均为N的复数序列{U_b[k]}，b∈I。计算每个复序列{U_b[k]}的平方功率谱，并相加则得到整个DNA序列S的功率谱序列{P[k]}：

$p\left[k\right]={\left|U_A\right[k\left]\right|}^2+{\left|U_T\right[k\left]\right|}^2+{\left|U_G\right[k\left]\right|}^2+{\left|U_C\right[k\left]\right|}^2,k=\mathrm{0,1},...N-1$   公式2

记DNA序列S的总功率谱的平均值为

$\stackrel{\‾}{E}=\frac{\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}P\left[k\right]}{N}$   公式3

而将DNA序列在特定位置，即

处的功率谱值，与整个序列S的总功率谱的平均值的比率称为DNA序列的“信噪比”（Signal Noise Ratio，SNR），即

$R=\frac{P\left[\frac{N}{3}\right]}{\stackrel{\‾}{E}}$ 公式4

信噪比R大于某个适当选定的阈值R₀（比如R₀＝2），是DNA序列上编码序列片段（外显子）通常满足的特性，而内含子则一般不具有该性质。由此作为判断外显子和内含子的重要区别。

由于这种方法基于DNA序列的功率谱P[k]和信噪比R，因此，研究如何提高计算功率谱P[k]和信噪比R的计算速度，进而提升此类方法的整体效率，将会是一个高度重视的问题。

发明内容

发明目的：本发明针对现有DNA序列方法中计算功率谱和信噪比效率低下的问题，提出了一种基于核苷酸位差的频谱3-周期性信噪比计算方法，通过先将DNA序列进行Voss映射得到4个指示序列，然后分别统计并记录4个指示序列里任意两个“1”出现位置之差，并把这些不同的距离值所出现的次数累积存储到一维数组中，再利用余弦函数的周期性，根据诱导公式，最终计算出该DNA序列的信噪比，同时提高了计算效率。

技术方案：一种基于核苷酸位差的频谱3-周期性信噪比计算方法，主要包括记录各个核苷酸出现位置之差（以下简称位差）并累积统计到数组d[i]中和计算余弦函数周期

并通过取余操作把数组d[i]进行重新归并两部分，其中：所述记录位差包括对DNA序列进行Voss映射、对指示序列中任何两个“1”之间的位差的记录、以及位差值的累积次数的统计，所述计算余弦函数周期t包括通过取余操作使数组d[i]所有下标都落入[1,t]范围内，并把原有的值累积加到取余结果后所得的下标所对应的值当中；

包括如下步骤：

步骤1，记录DNA序列长度N；

步骤2，进行Voss映射，得到四个序列{u_A[n]}、{u_G[n]}、{u_C[n]}、{u_C[n]}，分别记录每两个“1”在该序列中的位置之间的距离，并把所有距离值出现的累积次数保存在数组d[i]中，如d[2]=8表示进行Voss映射之后，得到的四个序列中数字“1”之间距离为2的总次数为8次；

步骤3，计算

当k＝0时直接转到第4步即可），对于数组d[i]的下标只要不落入[1,t]范围内，都要进行取余操作，使数组d[i]所有下标都落入[1,t]范围内，并把原有的值累积加到取余结果后所得的下标所对应的值（称为修正值）当中，例如，经计算t=4，取余前d[2]=8，d[6]=2，进行取余归并操作之后，则变成d[2]=10；

步骤4，k≠0时，构造行向量和列向量

$Q={[\cos \left(\frac{2\mathrm{\π}*k}{N}\right),\cos \left(\frac{2\mathrm{\π}*2k}{N}\right),...\cos \left(\frac{2\mathrm{\π}*\mathrm{tk}}{N}\right)]}^T,$ k＝0时，构造行向量

$L=\left[d\right[1],d[2],...d[N\left]\right]$ 和列向量 $Q={[\mathrm{1,1},...1]}^T;$

步骤5，利用公式P[k]＝N+2LQ计算整个DNA序列S的功率谱序列P[k]（k＝1,…,N-1）；

步骤6，计算总功率谱的平均值为

步骤7，计算信噪比

所述对DNA序列进行Voss映射定义如下：

在DNA序列研究中，首先需要把A、T、G、C四种核苷酸的符号序列，根据一定的规则映射成相应的数值序列，以便于对其作数字处理。

令I＝{A,T,G,C}，长度（即核苷酸符号个数）为N的任意DNA序列，可表达为

S={S[n]|S[n]∈I,n=0,1,2,...,N-1}

即A、T、G、C的符号序列S：S[0],S[1],…,S[N-1]。现对于任意确定的b∈I，令

$u_b\left[n\right]=\left\{\begin{array}{cc}1,& S\left[n\right]=b\\ 0,& S\left[n\right]\≠b\end{array},n=\mathrm{0,1,2},...N-1\right.$

称之为Voss映射，于是生成相应的0-1序列（即二进制序列）{u_b[n]}：

u_b[0],u_b[1],…,，u_b[N-1](b∈I)，也称为指示序列。

所述对指示序列中任何两个“1”之间的位差进行记录为分别在4个序列{u_A[n]}、{u_G[n]}、{u_C[n]}、{u_C[n]}中进行；

所述位差值的累积次数的统计定义为把4个指示序列里所有距离值出现的累积次数保存在数组d[i]中；

所述通过取余操作把使数组d[i]所有下标都落入[1,t]范围内，并把原有的值累积加到取余结果后所得的下标所对应的值当中定义为：对于数组d[i]，

，只要，则令i₂＝i₁modt，d[i₂]＝d[i₂]+d[i₁]，例如，经计算t=4，取余前d[2]=8，d[6]=2，进行取余归并操作之后，则变成d[2]=8+2=10。

本发明采用上述技术方案，具有以下有益效果：对于给定的DNA序列，采用核苷酸位置差，避免了DFT运算中的复数运算；利用余弦函数的周期性，大大减少了计算量，减少了计算该DNA序列的功率谱和信噪比所需时间，从而提高了计算效率。

附图说明

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

图2为现有技术中基于序列频谱3-周期性的基因识别方法

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

图4为基于核苷酸位差的频谱3-周期性信噪比算法与信噪比FFT算法时间性能对比图（其中，*表示FFT算法，o表示核苷酸位差的频谱3-周期性信噪比算法）。

具体实施方式

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

如图2所示，基于核苷酸位差的频谱3-周期性信噪比计算方法分为七个步骤：

步骤1，记录DNA序列长度N；

步骤2，进行Voss映射，得到四个序列{u_A[n]}、{u_G[n]}、{u_C[n]}、{u_C[n]}，分别记录每两个“1”在该序列中的位置之间的距离，并把所有距离值出现的累积次数保存在数组d[i]中，如d[2]=8表示进行Voss映射之后，得到的四个序列中数字“1”之间距离为2的总次数为8次；

步骤3，计算当k＝0时直接转到第4步即可），对于数组d[i]的下标只要不落入[1,t]范围内，都要进行取余操作，使数组d[i]所有下标都落入[1,t]范围内，并把原有的值累积加到取余结果后所得的下标所对应的值（称为修正值）当中，例如，经计算t=4，取余前d[2]=8，d[6]=2，进行取余归并操作之后，则变成d[2]=10；

步骤4，k≠0时，构造行向量

和列向量 $Q={[\cos \left(\frac{2\mathrm{\π}*k}{N}\right),\cos \left(\frac{2\mathrm{\π}*2k}{N}\right),...\cos \left(\frac{2\mathrm{\π}*\mathrm{tk}}{N}\right)]}^T,$ k＝0时，构造行向量 $L=\left[d\right[1],d[2],...d[N\left]\right]$ 和列向量 $Q={[\mathrm{1,1},...1]}^T;$

步骤5，利用公式P[k]＝N+2LQ计算整个DNA序列S的功率谱序列P[k]（k＝1,…,N-1）；

步骤6，计算总功率谱的平均值为

步骤7，计算信噪比

下面举例（为方便叙述，截取一段人类线粒体基因NC_012920_1.fasta上的长度为200的DNA序列进行）说明以上7个步骤：

步骤1，记录N＝200；

步骤2，对DNA序列进行Voss映射，得到4个指示序列{u_A[n]}、{u_G[n]}、{u_C[n]}、{u_C[n]}，分别在每个序列里求任意两个“1”的位置差，并把所有距离值出现的次数累积存储到d[i]中，如下表所示

表14个指示序列中所有距离值累积次数

步骤3，k＝1时，t＝N/k＝200，对于d[i]，

，i∈[1,200]，所以无需进行取余操作。而当k＝2时，有t＝N/k＝100，对于d[i]，

，

，

则进行令j＝imod100d[j]＝d[j]+d[i]，其结果如下表所示

表2k＝2时4个指示序列中所有距离值累积次数修正值

而当k＝3,4…N-1时，情况与上类似；

步骤4，k＝1时，构造行向量L＝[47,55,50,51…,0,0,0]和列向量

$Q={[\cos \left(\frac{2\mathrm{\π}*1}{200}\right),\cos \left(\frac{2\mathrm{\π}*2*1}{N}\right),...\cos \left(\frac{2\mathrm{\π}*200*1}{N}\right)]}^T={[\cos \left(\frac{\mathrm{\π}}{100}\right),\cos \left(\frac{\mathrm{\π}}{50}\right),...\cos \left(2\mathrm{\π}\right)]}^T,$

k＝2时，构造行向量L＝[74,78,70,67…,30,23,31]和列向量

$Q={[\cos \left(\frac{2\mathrm{\π}*2}{200}\right),\cos \left(\frac{2\mathrm{\π}*2*2}{N}\right),...\cos \left(\frac{2\mathrm{\π}*100*2}{N}\right)]}^T={[\cos \left(\frac{\mathrm{\π}}{50}\right),\cos \left(\frac{\mathrm{\π}}{25}\right),...\cos \left(2\mathrm{\π}\right)]}^T,$

而当k＝3,4…N-1时，情况与上类似；

步骤5，k＝1时，P[1]＝N+2LQ＝10226，k＝2时，P[1]＝N+2LQ＝164.0716，

对于k＝3,4…N-1时，情况类似；

步骤6，计算总功率谱的平均值为

步骤7，计算信噪比

如图4所示，为使用本计算方法与使用信噪比FFT算法运行效率对比图（其中，*表示FFT算法，o表示本算法），实验结果表明，本方法比传统的FFT算法所需时间较少，而且DNA长度越大，本方法在时间效率上的优势越明显。

Claims (3)

1.一种基于核苷酸位差的频谱3-周期性信噪比计算方法，其特征在于：主要包括记录各个核苷酸出现位置之差并累积统计到数组d[i]中和计算余弦函数周期

并通过取余操作把数组d[i]进行重新归并两部分，其中：所述记录位差包括对DNA序列进行Voss映射、对指示序列中任何两个“1”之间的位差的记录、以及位差值的累积次数的统计，所述计算余弦函数周期t包括通过取余操作把使数组d[i]所有下标都落入[1,t]范围内，并把原有的值累积加到取余结果后所得的下标所对应的值当中；

具体包括如下步骤：

步骤1，记录DNA序列长度N；

步骤2，对DNA序列进行Voss映射，得到四个指示序列{u_A[n]}、{u_G[n]}、{u_C[n]}、{u_C[n]}，分别记录每两个“1”在所述序列中的位置之间的距离，并把所有距离值出现的累积次数保存在数组d[i]中；

步骤3，计算

当k＝0时直接转到第4步即可），对于数组d[i]的下标只要不落入[1,t]范围内，都要进行取余操作，使数组d[i]所有下标都落入[1,t]范围内，并把原有的值累积加到取余结果后所得的下标所对应的值当中；

步骤4，k≠0时，构造行向量L＝[d[1],d[2],…d[t]]和列向量 $Q={[\cos \left(\frac{2\mathrm{\π}*k}{N}\right),\cos \left(\frac{2\mathrm{\π}*2k}{N}\right),...\cos \left(\frac{2\mathrm{\π}*\mathrm{tk}}{N}\right)]}^T,$ k＝0时，构造行向量 $L=\left[d\right[1],d[2],...d[N\left]\right]$ 和列向量 $Q={[\mathrm{1,1},...1]}^T;$

步骤5，利用公式P[k]＝N+2LQ计算整个DNA序列S的功率谱序列P[k]（k＝1,…,N-1）；

步骤6，计算总功率谱的平均值为

步骤7，计算信噪比

2.如权利要求1所述的基于核苷酸位差的频谱3-周期性信噪比计算方法，其特征在于：所述对DNA序列进行Voss映射定义如下：

在DNA序列研究中，首先需要把A、T、G、C四种核苷酸的符号序列，映射成相应的数值序列；

令I＝{A,T,G,C}，长度为N的任意DNA序列，可表达为

S＝{S[n]|S[n]∈I,n＝0,1,2,…N-1}

即A、T、G、C的符号序列S：S[0],S[1],…,S[N-1]；现对于任意确定的b∈I，令

$u_b\left[n\right]=\left\{\begin{array}{cc}1,& S\left[n\right]=b\\ 0,& S\left[n\right]\≠b\end{array},n=\mathrm{0,1,2},...N-1\right.$

称之为Voss映射，于是生成相应的0-1序列{u_b[n]}：u_b[0],u_b[1],…,，u_b[N-1](b∈I)，也称为指示序列；

所述对指示序列中任何两个“1”之间的位差进行记录为分别在4个序列{u_A[n]}、{u_G[n]}、{u_C[n]}、{u_C[n]}中进行；

所述位差值的累积次数的统计定义为把4个指示序列里所有距离值出现的累积次数保存在数组d[i]中。

3.如权利要求2所述的基于核苷酸位差的频谱3-周期性信噪比计算方法，其特征在于：所述通过取余操作使数组d[i]所有下标都落入[1,t]范围内，并把原有的值累积加到取余结果后所得的下标所对应的值当中定义为：对于数组d[i]，只要则令i₂＝i₁modt，d[i₂]＝d[i₂]+d[i₁]。

Images