基因识别问题及其算法实现.doc

基因识别问题及其算法实现摘 要：本文运用信号系统理论、相关性分析、多目标规划、fisher判别等知识研究了基因识别问题，提出了一种基于频谱峰值位置的基因识别方法。针对问题一，总功率谱的平均值通过Parseval定理在时域中计算，功率谱和信噪比转换为矩阵乘法运算，但矩阵乘法运算量也很大。统计1、2、3的倍数位置上各碱基数目然后与相乘可减去计算量。利用NC_012920_1基因数据求功率谱和信噪比，三种方法所需时间分别是3.7s、2.9s、0.36s。利用Z-curve映射和Voss映射的关系式，求得Z-curve映射的频谱是Voss的4倍，信噪比为4/3倍，同时发现Z-curve映射具有更强的生物意义。对实数映射进行相关分析得到外显子和内含子频谱图明显不同，可用来识别两者。针对问题二，利用外显子的频谱3-周期性，以保守的阈值判别基因序列，低于阈值的基因，以很大概率为内含子，直接去除，剩下序列中外显子比例提高，增强3-周期性；该方法采用自适应的阈值，具有更强的适应性并通过实验得到验证。我们选取了四个性能指标：敏感度、特异度、近似相关系数、相关系数，通过加权将多指标转换为单一的综合指标。考虑到评价近似相关系数、相关系数为综合指标，将四者权值分别设为40%、40%、10%和10%。针对问题三，创新性地提出了一种基于频谱峰值位置的基因识别方法，该方法有较高的性能指标，对部分基因的外显子端点能作出较为准确的识别，分析了影响该方法性能的因素，并对未注释的DNA序列进行了编码区域的预测。提出了一种通过添加外显子基因的方式改善某些序列3-周期性质不显著的问题。针对问题四，考虑到基因序列的其他特性，如碱基组成成分，碱基位置相关性，密码子使用偏好性等，利用前面三个特征再结合3-周期性构造统计特征量；使用fisher判别的方法，对原数据进行坐标变换，再借助方差分析的思想构造一个判别函数。对于可能突变了的基因，若有其正常序列，进行序列比对即可；若无则借助频谱和其他特征量进行综合判别。最后，对本文建立的几个模型进行评价，提出了推广的方向。关键词：峰值位置识别 3-周期性 Parseval定理 自适应阈值 fisher判别- 27 -目 录1、问题重述- 1 -1.1、相关材料- 1 -1.2、问题提出- 1 -2、问题分析- 2 -3、模型假设- 3 -4、符号说明- 3 -5、模型的建立与求解- 4 -5.1、问题1的模型建立及求解- 4 -5.1.1、功率谱与信噪比的快速计算方法- 4 -5.1.2、Z-curve映射Voss映射的关系- 8 -5.1.3、实数映射的关系- 11 -5.2、问题2的模型建立及求解- 12 -5.2.1、阈值的确定- 12 -5.2.1、评价标准- 13 -5.3、问题3的模型建立及求解- 15 -5.3.1、基于频谱峰值相对位置的基因识别方法的提出- 15 -5.3.2、峰值位置识别方法的实验结果- 17 -5.3.3、识别性能的影响因素- 18 -5.3.4、确定端点位置- 21 -5.3.5、 对未被注释的DNA序列编码区域的预测- 22 -5.4、问题4的模型建立及求解- 23 -5.4.1、识别基因编码序列的其它特征- 23 -5.4.2、识别基因编码序列的基因突变- 25 -6、模型的评价与改进- 25 -6.1、模型的评价- 25 -6.2、模型的推广- 25 -参考文献- 25 -1、问题重述DNA是生物遗传信息的载体，其由腺嘌呤（Adenine,A），鸟嘌呤（Guanine,G），胞嘧啶（Cytosine,C），胸腺嘧啶（Thymine,T）这四种核苷酸按一定的顺序连接而成。这些长链上不仅包含制造人类全部蛋白质的信息，还有按照特定的时空模式把这些蛋白质装配成生物体的四维调控信息（三维空间和一维时间）。如何找到这些信息的编码方式和调节规律，是了解遗传本质的核心，也是医学研究的重要方向之一。1.1、相关材料其中带有遗传讯息的DNA片段称为基因（见图1第一行）。其他的DNA序列片段，有些直接以自身构造发挥作用，有些则参与调控遗传讯息的表现。在真核生物的DNA序列中，基因通常被划分为许多间隔的片段（见图1第二行），其中编码蛋白质的部分称为外显子，不编码的部分称为内含子。DNA序列外显子(Exon) 内含子(Intron)基因（Gene） 图1 真核生物DNA序列（基因序列）结构示意图对给定的DNA序列，怎么去识别出其中的编码序列（即外显子），也称为基因预测，是一个尚未完全解决的问题，也是当前生物信息学的一个最基础、最首要的问题。基因预测问题的一类方法是基于统计学的。然而统计预测方法通常需要将编码序列信息已知的DNA序列作为训练数据集来确定模型中的参数，从而提高模型的预测水平。但在对基因信息了解不多的情况下，基因识别的准确率会明显下降。因此在目前基因预测研究中，采用信号处理与分析方法来发现基因编码序列也受到广泛重视。1.2、问题提出基因识别研究中，首先需要把A、C 、G、T四种核苷酸的符号序列，根据一定的规则映射成相应的数值序列，以便于对其作数字处理。在使用Voss映射时，发现DNA序列中的外显子有明显的3-周期性，而内含子没有。为了得到很长DNA序列的功率谱或信噪比时，需要进行离散Fourier变换(DFT)，其总体计算量很大，会影响到所设计的基因识别算法的效率，需要提出一种基于此的快速算法。除了Voss映射外，实际上人们还研究过许多不同的数值映射方法。如著名的Z-curve映射，试探讨Z-curve映射的频谱与信噪比和Voss映射下的频谱与信噪比之间的关系。而对于实数映射，如：，是否也有功率谱与信噪比的快速计算公式？对不同的基因类型，所选取的判别阈值也许应该是不同的。通过有代表性的基因序列，探求每类基因的阈值确定方法和阈值结果。分析编码与非编码区间分类方法的有效性和分类错误。目前基因识别方面的多数算法结果还不是很充分。如通过频谱或信噪比特征进行基因识别的算法，由于DNA序列随机噪声的影响等原因，还很难“精确地”确定基因外显子区间的两个端点。对此，设计的基因识别算法并对其准确率做出适当评估，并对附件中给出的6个未被注释的DNA序列的编码区域进行预测。采用频谱或信噪比这样单一的判别特征，将影响、限制基因识别正确率。总结或提出一些识别基因编码序列的其它特征指数，并对此做相关的分析？能否利用频谱或信噪比方法去发现基因编码序列可能存在的突变呢？2、问题分析对于问题一，传统方法中傅里叶变换需同时计算N个点的频谱，如果窗口滑动P次，则要重复计算N个点的频谱P次，计算量非常大。事实上功率谱只是Voss映射下三分之一处的功率谱，而总功率谱的平均值可以通过Parseval定理在时域中计算，那么功率谱和信噪比可以转换为矩阵乘法运算。但矩阵乘法运算的计算量也很大，通过统计在1、2、3的倍数位置上A、G、T、C的数目然后与相乘即可去计算所需的功率谱，进而求信噪比。改进后的快速算法只需进行3次乘法和2次加法即可求得。而求得功率谱只需进行3*4=12次乘法和2*4+3=11次加法。探求Z-curve映射的频谱与信噪比和Voss映射下的频谱与信噪比之间的关系可以从Z-curve映射和Voss映射之间的关系式入手，利用每一位置只能取A、G、T、C中的一个核苷酸性质求得。对于实数映射，实际上属于一维映射，它包含了DNA序列中的所有信息，可对其进行相关分析探求序列内在关系。对于问题二，我们主要从外显子在频谱3-周期性角度考虑，基因的编码采取Voss映射，以极保守的阈值，阈值以下的基因，以很大可能性是内含子，可以去除掉。而在其他的位置上，也会因为噪声或计算不准确等原因，外显子并不能如愿的凸现出来。保守阈值以下的基因我们都认为是内含子，将他们的Voss映射的四个指示函数的相应位置上全部置为1。其他位置的编码保持不变，或加强。避免了对不同的基因，确定不同阈值繁琐的计算。对以上步骤执行若干次迭代，观察训练效果。我们选取了四个性能指标：敏感度、特异度、近似相关系数、相关系数，通过加权将多指标转换为单一的综合指标。对于问题三，基于固定长度滑动窗口上频谱曲线的基因识别方法和基于DNA序列上“移动序列”信噪比曲线的基因识别方法，从问题说明文件可以看出，它们判别外显子的判据十分模糊，对外显子端点的界定十分粗糙，并且有经验型的判别阈值，在识别基因时会遇到诸多困难。基于统计的方法，我们认为应该提出新的区别于上述方法的识别判据，经过一些实验我们发现，基因序列频谱峰值的相对位置是一个性能良好的判据。对于问题四，采用频谱或信噪比单一判别特征识别短的编码序列效果并不理想，考虑到基因序列还有其他特性，如碱基组成成分，碱基位置相关性，密码子使用偏好性等，利用前面三个特征再结合3-周期性构造统计特征量提取序列中的四种特征量；对于一个基因，将特征参数代入求出y值，然后与判别阈值R进行比较，就可判别其为外显子还是内含子。对于可能突变了的基因，若有其正常序列，进行序列比对即可；若有其相似类种基因，利用基因间的相似性，大段基因突变通过频谱或者信噪比检验，一个或几个发生突变，借助其他特征量进行综合判别；若无该基因的任何先验信息，则应从其它方面去考虑识别基因突变。3、模型假设假设一：样本中的100个人和鼠类的，以及200个哺乳动物类的基因序列真实，来源可靠，能够作为检验模型准确性的样本； 假设二：序列中个别异常点不影响序列的外显子和内含子分界点；假设三：预测中外显子预测成内含子和内含子预测成外显子的代价相同；假设四：外显子总是存在偏向性；4、符号说明：核苷酸的Voss映射：固定窗口的长度： A的Voss映射：的傅里叶变换： C的Voss映射：的傅里叶变换： G的Voss映射：的傅里叶变换： T的Voss映射：的傅里叶变换： Z-curve映射之一：的傅里叶变换： Z-curve映射之一：的傅里叶变换： Z-curve映射之一：的傅里叶变换：Voss映射下的功率谱：Z-curve映射下的功率谱：Voss映射下的信噪比： Z-curve映射下的信噪比：Voss映射下总功率谱的平均值： Z-curve映射下总功率谱的平均值：被正确预测外显子碱基数目：被错误预测外显子碱基数目：被正确预测外显子碱基数目：被错误预测外显子碱基数目：预测评价指标敏感度：预测评价指标特异度：预测评价指标近似相关系数：预测评价指标相关系数：基因序列中碱基A的个数：基因序列中碱基A的个数：基因序列中碱基C的个数：基因序列中碱基C的个数：基因序列中碱基G的个数：基因序列中碱基G的个数：基因序列中碱基T的个数：基因序列中碱基T的个数5、模型的建立与求解5.1、问题1的模型建立及求解对于很长的DNA序列，在计算其功率谱或信噪比时，离散Fourier变换(DFT)的总体计算量仍然很大，会影响到所设计的基因识别算法的效率。5.1.1、功率谱与信噪比的快速计算方法用传统方法研究DNA编码序列（外显子）的特性时，对指示序列分别做离散Fourier变换(DFT) (5.1.1)以此可得到四个长度均为N的复数序列，。计算每个复序列的平方功率谱，并相加则得到整个DNA序列S的功率谱序列： (5.1.2)由于外显子序列的功率谱曲线在频率处，具有较大的频谱峰值，而内含子则没有类似的峰值。这种统计现象被称为碱基的3-周期。 记DNA序列的总功率谱的平均值为 (5.1.3)将DNA序列在特定位置，即处的功率谱值，与整个序列S的总功率谱的平均值的比率称为DNA序列的“信噪比”，即 (5.1.4)DNA序列的信噪比值的大小，既表示频谱峰值的相对高度，也反映编码或非编码序列3-周期性的强弱。信噪比R大于某个适当选定的阈值（比如），是DNA序列上编码序列片段（外显子）通常满足的特性，而内含子则一般不具有该性质。用传统预测方法对基因组序列的编码区进行预测和定位时，傅里叶变换是同时计算N个点的频谱，如果窗口滑动P次，则要重复计算N个点的频谱P次，计算量非常大1。即使是采用快速傅里叶变换，完成一个DNA序列编码区的预测也需要相当长的时间。重复计算不需要的N-1个点的频谱是一种资源浪费。在对Voss映射，求总功率谱时有(5.1.5)在信号处理中， Parseval定理2告诉我们，信号在时域的总能量等于其在频域的总能量。即在离散傅里叶变换中有 (5.1.6) 公式(5.1.4)对应的离散傅里叶变换为 (5.1.7)在这里，我们需要对Parseval定理进行修正，修正后的为 (5.1.8)在DNA序列中，每个位置只能是A、C 、G、T中的一个，那么对于Voss映射，中只有一个为1,其余都为0，故DNA序列的任意位置有 (5.1.9)根据总功率谱的计算公式，总功率为(5.1.10)Voss映射下的总功率谱的平均值为 (5.1.11)接下来只需求Voss映射下三分之一处的功率谱即可求得信噪比，由的定义式有从数字(5.1.12)看出，采用该快速算法，不需要重复计算N个点的频谱，只需对点计算频谱就能实现预测。定义则(5.1.13)在计算时，先将矩阵W存储在内存中，需要计算时，利用这个矩阵运算式子即可求得。故每计算求得一个，只需进行四次矩阵乘法和三次加法即可。相对于传统算法计算所有点的离散傅里叶，这里只需进行简单矩阵运算即可。在实际计算中，矩阵乘法运算的计算量也很大，这从原理上说明其不可能对程序运行速度改进很好，很难成为真正的“快速算法”。分析式子(5.1.13)可知，中只含有0,1两种数值，而矩阵w中只有三种取值，更有趣的是是重复规律出现，这样只需统计中在1、2、3的倍数位置上A的数目然后与w中的三种数值相乘之和的模平方即可。定义在1、2、3的倍数位置上A的数目为，则有(5.1.14)从式(5.1.14)得知，改进后的快速算法只需进行3次乘法和2次加法即可求得。而求得功率谱只需进行3*4=12次乘法和2*4+3=11次加法。利用题目给的NC_012920_1(人线粒体全基因组) 基因数据中3000-7000位置处对该方法进行了验证。利用matlab3中tic和toc函数得传统方法程序运行所需时间为“Elapsed time is 3.752879 seconds”，而该快速算法程序运行所需时间为“Elapsed time is 2.992869 seconds”，改近后的快速算法程序运行所需时间为“Elapsed time is 0.363487 seconds”。快速算法需进行四次矩阵乘法和三次加法，显然矩阵乘法制约其优势，相对于传统方法，其速度上并没有太大优势，只是在内存方面有所节约，但这不是我们所关心的。而改进后的算法则只需统计在1、2、3的倍数位置上A、C 、G、T的数目即可，不需要进行傅里叶变换和矩阵乘法运算，相对于传统算法，其性能提高十倍以上。若窗口滑动距离很短，所需进行的傅里叶变换越多，则该算法的快速性更能体现。图1为三种算法得到的频谱图，从图中看出三幅图像的总体趋势相同，只是细节上有所区别，这说明改进算法得到的功率谱图具有实际意义。图a 传统算法得到的功率谱图图c 改进后的快速算法得到的功率谱图图b 快速算法得到的功率谱图 图5-1 三种算法得到的频谱图5.1.2、Z-curve映射Voss映射的关系由附件一中资料可得，Z-curve 映射有表达式4: (5.1.15)傅里叶变换具有线性叠加性5，时域中的式(5.1.15)在频域中有 (5.1.16)根据的定义式子有 (5.1.17)将式 (5.1.16)代入(5.1.17)得(5.1.18)因为在序列的任意位置有 (5.1.19)所以由此得到为序列1在值的傅里叶变换值，而由本题傅里叶变换式(5.1.1)可知序列为常数1的傅里叶变换值为，故在k（）值的傅里叶变换值为0。重写式子(5.1.18)有从式(5.1.21)上面结果看出，Z-curve映射的频谱与Voss映射下的频谱的关系为： (5.1.22)利用题目给的NC_012920_1(人线粒体全基因组) 基因数据求得和的比值如图5-2所示。图5-2 和的比值下面求Z-curve映射的信噪比与Voss映射下的信噪比之间的关系。要求信噪比，需先求得总功率谱的平均值，即 (5.1.23)对于Voss映射，总功率谱有 (5.1.24)在离散傅里叶变换中，有修正后的Parseval定理式(5.1.8)。而在时域方面中只有一个为1,其余都为0，联系是(5.1.9)得总功率谱为(5.1.25)同理，对于Z-curve映射，总功率谱有 (5.1.26)考虑到基因序列A、C、G、T四种核苷酸对应Z-curve映射如表5-1。表5-1 核苷酸对应的Z-curve映射AC GT1-11-111-1-11-1-11不管基因序列取A、C、G、T中的哪一个，都有 (5.1.27)再考虑到离散傅里叶中修正后的Parseval定理式(5.1.8)，得故对于信噪比有 (5.1.29)从上面结果看出，Z-curve映射的信噪比与Voss映射下的信噪比之间的关系为。在计算功率谱和信噪比中，由表5-1发现对于DNA序列第n个位置，若存在碱基A，则；若存在碱基T，则，；若存在碱基C，则；若存在碱基G，则，即的取值只能是-1和1。Z-curve在三个映射序列上的投影有着明确的生物学意义。其中分量表征着嘌呤/嘧啶碱基沿序列的分布；表征着氨基/酮基碱基沿序列的分布；分量表征着强氢键/弱氢键碱基沿序列的分布，这三种独立的分布完整的描述了所研究的DNA序列。可见，Z-curve不仅可以把DNA的字符序列转换成数字序列，而且其三个映射序列都有明确的生物学意义。5.1.3、实数映射的关系在Voss映射中，我们使用的是四维映射，即一个DNA序列映射出四组0-1序列，而在Z-curve映射中则是采用三维映射。对于实属映射，如：，我们可得到其一维映射。一维映射与多维不同，它与核苷酸序列是一一对应关系。同时根据有关研究发现，基因外显子片段上，四种核苷酸符号在序列的三个子序列上分布的“非均衡性”使得基因序列具有偏向性。在这我们可以理解为外显子的核苷酸序列具有某种内在的联系，而内含子则无此特性。外显子的这种内在联系可以利用通过自相关序列进行反映出来，这里定义的自相关为 (5.1.30)将时域中的自相关转换到频域中观察，发现3-周期现象依然明显，因此可以采用自相关法来计算实数序列的自相关，然后再利用快速傅立叶变换来观察其频谱在三分之一处是否有一频谱峰值。利用题目给的NC_012920_1(人线粒体全基因组) 基因数据对该方法进行了验证。基因序列在3300处属于外显子，其相关函数的傅立叶变换如图所示；基因序列在4300处属于内含子，其相关函数的傅立叶变换如图所示。图5-3 外显子的功率谱图5-4 内含子的功率谱从上图看出，外显子在三分之一处有一频谱峰值，而内含子没有，这样就可以通过该性质来区分外显子和内含子。5.2、问题2的模型建立及求解我们主要从外显子在频谱3-周期性角度考虑，基因的编码采取Voss映射，以极其保守的阈值判别基因序列，低于阈值的基因，以很大概率为内含子，直接去除，剩下序列中外显子比例提高，增强了序列的3-周期性。适当提高阈值再进行以上步骤，得到外显子区间；该方法采用自适应的阈值，具有更强的适应性并通过实验证明了该方法的有效性。我们选取了四个性能指标：敏感度、特异度、近似相关系数、相关系数，通过加权将多指标转换为单一的综合指标。5.2.1、阈值的确定在保守阈值Thred以下的基因我们都认为是内含子，将他们的Voss映射的四个指示函数的相应位置上全部置为1。其他位置的编码保持不变，或加强。对以上步骤执行若干次迭代，观察训练效果。我们认为本节所提出的方法主要有以下创新：首先，不改变基因的相对位置，去除内含子可能存在的噪声。如果将这个过程比作一个滤波器，每次迭代都将底部的确定无疑的内含子信号滤掉，同时过滤掉的还有内含子所包裹的噪音信号。其次，不改变外显子上的碱基出现频率的偏向性。这就保证了不改变基因的3-周期特性。为本方法的实现提供了理论上实现的可能。最后，本方法不需要样本数据，易于实现。其结果使得外显子在频谱图上，不断地凸现出来，这对于阈值的选择就有比较大的空间。即使在选择阈值时没有理论上的依据，同样可以凭借经验得出较好的阈值选择。同时，避免了对不同的基因，确定不同阈值繁琐的计算。5.2.1、评价标准为了对不同的基因预测算法进行评价，我们在采取以下评价四个指标：敏感度、特异度、近似相关系数、相关系数敏感度。四个评价指标计算公式如下 (5.2.1) (5.2.2) (5.2.3) (5.2.4) (5.2.4)其中相关系数用来衡量序列全局准确性(即包括编码区,也包括非编码区的准确度)，但是它在一些极端情况下不可计算，所以在这我们采取了一个新的度量标，即近似相关系数。和都是作为核苷酸水平上全局准确性的度量标准，在大部分对算法进行评价的文章中都同时使用这两个度量6。本文求算功率谱采取固定窗口方法，所采取的“保守的阈值”为0.15，在进行了迭代加强之后，判定外显子的阈值可以定在0.2到0.35之间，在这个范围之间的阈值的不同，并不能导致性能结果有很大的不同。下面运用此方法随机在众多的基因库中随机地选择数据进行试验验证。图5-5到图5-8为两种方法得到的频谱图；表5-2为两种方法得到的性能指标。图a 固定阈值的频谱曲线图b 自适应阈值的频谱曲线图5-5 两种方法计算genes200.gene(1,33)的频谱曲线图a 固定阈值的频谱曲线图b 自适应阈值的频谱曲线图5-6 两种方法计算genes200.gene(1,59)的频谱曲线图a 固定阈值的频谱曲线图b 自适应阈值的频谱曲线图5-7 两种方法计算genes200.gene(1,87)的频谱曲线图a 固定阈值的频谱曲线图b 自适应阈值的频谱曲线图5-8 两种方法计算genes200.gene(1,138)的频谱曲线表5-2 两种方法得到的性能指标数据来源计算方法33号基因固定阈值0.731050.338380.0862260.0862160.4450自适应阈值0.785340.346060.188370.186760.490159号基因固定阈值0.470020.0971260.0935880.0780760.2440自适应阈值0.911270.0991910.189170.156110.438787号基因固定阈值0.608660.307930.163830.163220.3993自适应阈值0.308120.392870.392870.383620.6009138号基因固定阈值0.365230.368360.484720.470540.6429自适应阈值0.845470.38060.226070.22600.5356从上面结果看出，采用自适应阈值判别基因序列的方法，图像的总体趋势保持不变，证明在模型上大的假设成立。虽然在个别指标上（如87号基因上的敏感度指标）劣于固定门限，但从综合评价指标来看，都高于固定阈值判别方法。说明改善了方法性能。这是由于低于阈值的基因，以很大概率为内含子，直接被去除，剩下序列中外显子比例提高，增强了序列的3-周期性。运用对序列的处理， 扩大了阈值选择的范围，增强了阈值选择的适应性。5.3、问题3的模型建立及求解5.3.1、基于频谱峰值相对位置的基因识别方法的提出从前面可以看出，把频谱曲线和信噪比曲线作为判据，判断外显子和内含子不是十分明确，这里提出使用频谱最大值位置作为判据。由资料可知，基因外显子序列的功率谱曲线，在频率处具有较大的频谱峰值，即是，一段外显子序列的功率谱曲线的最大值，一般出现在频率处，内含子则没有这种性质；我们把最大值位于整段谱线的相对位置作为区分外显子和内含子的判据。对一个DNA序列和它的指示序列，。取长度M（通常取为3的倍数，例如M=99, 129, 255, 513等）作为固定窗口长度。对任意n（），在以n为起始的长度为M的序列片段上，作四个指示序列的离散Fourier变换(DFT) (5.3.1)相加则得到整个DNA序列的功率谱序列： (5.3.2)这里，求出它在固定时的最大值位置，即解，得到，定义函数，并画出其频谱曲线。下图为基于固定长度滑动窗口上频谱曲线的基因识别方法，基于固定长度滑动窗口上频谱峰值相对位置的基因识别方法，对NC_012920_1(人线粒体全基因组)基因序列的处理结果对比。图5-9 频谱法和峰值位置法对比上方为频谱方法的曲线以及表示外显子位置的线段，下方为峰值位置方法的曲线。从图中可以看出，频谱曲线方法比较难以对外显子进行识别，频谱峰值位置的方法识别较为明显，可以看出值为1/3处的曲线和外显子区间有比较明显的对应关系。这里因为峰值的位置几乎不会等于1/3，所以应设定一个判定的阈值，即。我们得到满足阈值的一个点的位置集合，因为外显子一般是连续的，先排除掉孤立点，定义孤立点满足且（因为外显子的3-周期性质，可能平移1、2个序列不表现性质），剩下一个新的位置集合，认为所有位置为外显子。定义一个与序列对应的函数，时表示点为外显子；时表示点为内含子。令所有满足的，其余位置，即完成识别算法。5.3.2、峰值位置识别方法的实验结果在阈值，窗口（不变情况下不再提及）的情况下，用峰值位置识别方法对NC_012920_1(人线粒体全基因组)的识别结果如下：图5-10 峰值位置法对NC_012920_1的处理结果值为1处为识别的外显子位置，0处为识别的内含子位置，0.5处为正确的外显子位置，以后的图类似，不再赘述。从性能评价指标和图，可以看出对NC_012920_1基因的识别较好。该方法对genes200中的11号基因的处理结果：图5-11 峰值位置法对genes200中的11号基因的处理结果5.3.3、性能影响因素我们认为，对识别性能的影响因素主要有：基因本身的外显子内含子分布，截取窗口长度，设定的阈值。A、基因本身的影响基因本身的影响因素，我们认为主要表现在外显子和内含子的长度上，一般来说，外显子、内含子越长，性能越好；外显子、内含子越短，性能越差。Gene100中的31号基因，外显子长度平均为116，最小为33；图5-12 峰值位置法对genes100中的31号基因的处理结果Gene100中的90号基因，外显子长度平均为648.5，最小为591；图5-13 峰值位置法对genes100中的31号基因的处理结果我们认为，3-周期性质是一种统计规律，一般要在样本较大时才能体现得较为明显，在外显子（或者内含子）较短的情况下，3-周期性质可能不能很好地体现。B、窗口长度的影响窗口长度的影响较为复杂，与外显子的长度有关。对Gene200中的77号基因，其外显子平均长度为155，最短为37。当时，识别如下：图5-14 峰值位置法对genes200中的77号基因的处理结果当时，识别如下：图5-15 峰值位置法对genes200中的77号基因的处理结果可以看出窗口长度增加了，性能降低了。对Gene200中的99号基因，其外显子平均长度为396.5，最短为355。当时，识别如下：图5-16 峰值位置法对genes200中的99号基因的处理结果当时，识别如下：图5-17 峰值位置法对genes200中的100号基因的处理结果可以看出窗口长度增加了，性能稍有提高。结合2个基因的长度，可以得出一个定性的结论：在基因外显子较短时，窗口长度增加，性能降低；在基因外显子较长时，窗口长度增加，性能提高。这里可以通过窗口长度与研究外显子基因平均长度的比值对性能的影响做进一步的试验研究。C、阈值对性能的影响在阈值，窗口下，对genes200中的11号基因的处理结果：图5-18 峰值位置法对genes200中的11号基因的处理结果性能几乎没有影响。在阈值下，对genes200中的11号基因的处理结果：图5-19 峰值位置法对genes200中的11号基因的处理结果可以看到性能有很大降低。当外显子的3-周期性质表现明显时，频谱的峰值位置精确地在1/3处，所以当阈值减小时，性能变化不明显；但是当阈值增大时，有可能导致内含子误判为外显子。5.3.4、确定端点位置对于已经识别完成的外显子标志函数，只需在且出标记为外显子开端，在且出标记为外显子末端即可。对之前识别过的genes200中的11号基因，其外显子端点序列为识别的外显子端点序列为平均偏差为16，平均绝对值偏差为，方差为972。这是比较少见的识别端点较为准确的情况，一般情况下端点的识别偏差非常大。识别genes200中的55号基因，其外显子端点序列为识别的外显子端点序列为图5-20 峰值位置法对genes200中的55号基因的端点识别可以看到，端点比较准确，但是少识别了一段外显子。5.3.5、预测DNA编码区域使用峰值位置方法进行编码区预测，对genes6中的1号基因预测如下： 图5-21 峰值位置法对genes6中的1号基因的预测外显子端点序列为对所有基因的外显子预测在附件中。5.3.6、基因添加方法在外显子较短的情况下，序列很可能表现不出3-周期性，因为这是一种统计规律，所谓的“大数效应”，在样本较少的情况下，统计规律可能并不明显，这是使用统计方法识别基因必然会遇到的问题。我们考虑这样一种方法，在一段残缺的基因序列后认为地添加一段基因序列，以使得它表现出3-周期性。首先实验验证其可行性。我们在gene100文件中收集了338个外显子和内含子序列（见数据文件gene100外显子序列.mat，gene100内含子序列.mat中储存的ex，in数据），并且准备了一段待添加的基因ex(259)；我们对这个基因序列进行了峰值位置判定，再把这些基因加上基因ex(259)后，进行峰值位置判定，结果如下：添加前有106个外显子、2个内含子满足判据（即符合判定为外显子的条件），其差为104；添加后有186个外显子、21个内含子满足判据，其差为165；这里举一个典型的例子：图5-22 ex(7)添加基因前图5-23 ex(7)添加基因后这里的ex(7)外显子基因长度为73，属于较短的外显子基因；可以看出e(7)经过基因添加后，峰值明显地移到了1/3处。经过基因添加后，满足判据的外显子明显增多，而满足判据的内含子变化较小，识别方法的性能得到了一定的提高。由此可以看出，外显子基因添加方法能提高识别方法的性能，可能是一种改善短外显子3-周期性质不显著问题的有效途径；但是在这里，因为时间关系，如何挑选添加基因，添加基因的影响等问题，我们没有进一步地深入研究。5.4、问题4的模型建立及求解5.4.1、识别其它特征人们发现，对某些DNA序列而言，其部分编码序列（外显子），尤其是短的（长度小于100bp）的编码序列，就可能不具有频谱或者信噪比显著性。采用频谱或信噪比这样单一的判别特征，也许是影响、限制基因识别正确率的一个重要原因。考虑到基因序列的多种特性，如碱基组成成分，碱基位置相关性，密码子使用偏好性等7，利用基因序列的这些特性，再结合3-周期性来提高较短基因的识别精度。下面介绍四个统计特征量的构造。密码子使用频率是基因编码区的一个重要统计特征。基因在编码生成蛋白质的过程中，每三个碱基组成一个密码子，转换成个氨基酸。碱基包含A、C、G、T四种，氨基酸包含有20种，即共有4x4x4=64种不同的密码子组合对应20种不同的氨基酸。对于一条基因编码序列，假设长度L为3的整数倍，将每三个相连的碱基作为个密码子，统计出各个密码子的数量，比如密码子AAA的个数记作，则该密码子的使用频率为，定义密码子使用的频率特征向量为，其中D属于密码子集合，。蛋白质编码区和非编码区所包含的碱基成分是不相同的，根据生物学理论可知蛋白质编码区中的C、G含量较高，而非编码区中的A、T含量较高。故碱基组成成分是一个重要的统计特征。四为种碱基的个数分别，则四种碱基对应的含量可表示为，碱基组成成分的统计特征量为。各种碱基的成分作为一种统计特征只反映了碱基的百分含量，而没有体现碱基出现的周期性，即碱基出现位置的相关性。碱基A出现的位置分别记做，取相连的两个位置之间的差，这个差值体现了该碱基出现位置的相关性。位置差值组成的序列为计算序列的二阶中心矩 (5.4.1)其中，由统计学知识，能体现该碱基出现的周期性，波动幅度越小，说明位置的差值变化越小，即碱基出现的周期性越强；反之，则说明该碱基出现的周期性越弱。定义碱基位置的相关性特征向量。考虑到该特征值相对以上两种特征量数值大很多，所以在算法中对每个向量都进行了归化。3-周期性前面第一问已经详细说明，考虑到该功率谱数值较大，为了和其他几个特征量结合使用，在算法中对特征量取自然对数，即。选取基因识别的特征量后，把DNA序列分成两类，即外显子和内含子，分别记为P和N。提取序列中的四种特征量，其中包含64个参数，包含4个参数，包含4个参数，包含4个参数，使用fisher判别的方法，对原数据进行坐标变换，寻求能将总体尽可能分开的方向。再借助方差分析的思想构造一个判别函数 (5.4.2)其中向量C确定的原则是使得两个总体组间的区别最大，同时每个组内部的离差最小。确定了判别函数