基于Bootstrap估计的基因识别方法

【摘要】随着1980年第一个完整的生物基因序列的发表，生命科学的研究进入了新纪元。人类基因组计划的启动、大肠杆菌基因组的测序、人类基因组序列图的绘制成功以及人类基因组测序工作的完成等重大科研成果的相继出现，正引领人们逐步了解生物基础物质中承载和包含的信息。同时基因研究的成果也正应用到疾病攻克、刑事侦查、食品安全和防止物种入侵等多个社会生活领域。在海量生物基因数据中，寻找和发现蕴含的深层信息，已成为当前研究的重点领域。将DNA序列信息作为分析的源头，找到编码区域，对于帮助我们认识遗传、代谢、发育、分化和进化等重要规律有重大意义。本文将数学模型应用于基因识别及算法实现问题中，运用计算机编程手段，对基因识别相关问题进行分析。

【关键词】Bootstrap 基因识别 阈值

一、采用信噪比进行DNA识别的阈值确定方法

（一）信噪比的计算

信噪比的计算由公式

构成，计算过程由MATLAB软件实现。

（二）阈值确定方法——基于Bootstrap抽样法

DNA序列的编码区识别研究中，作为一个广泛使用的特征指数，信噪比的阈值的选取一直是众多学者所分析和阐述的重点，不同物种类型基因、不同序列长度的编码区间，阈值的大小往往不同。本文采用统计学中非参数思想方法，Bootstrap有放回重复抽样方法，对所给的100个人和鼠类的以及200个哺乳动物类的基因序列进行阈值的分析和确定，得到比较合理科学的结论。

在基因外显子序列的识别过程中也同样伴随着先验信息缺乏和样本量不足的问题，信噪比的形式和计算方法也是影响阈值和最终结果的重要因素，多重因素相互影响，导致不同物种基因识别时阈值的确定成为一大难题。本文将Bootstrap方法引入到信噪比阈值的确定过程中，借助有放回重复抽样的非参数统计思想来提高阈值确定的精准性。

本文主要以100个人和鼠类的基因序列信息为例，附带分析200个哺乳动物类基因序列，采用bootstrap方法为该类生物推定最佳阈值，MATLAB软件和R软件结合使用，确定一个较为准确的阈值。

具体算法步骤设计和流程图如下：

1.选取基因序列。将100个人和鼠类的外显子序列作为基础数据，选取全部438个外显子序列；

2.计算原始样本。对于每一段DNA的外显子序列，运用MATLAB软件FFT函数计算各序列信噪比，即每个外显子序列判别时最大阈值，则长度为438的所有外显子序列的信噪比组成原始样本，且该样本含有待估计的最佳阈值的分布特性。

3.抽取Bootstrap样本。基于原始样本，设置Bootstrap抽样次数，利用R软件进行统计模拟，随机有放回的进行抽样，获得1000个Bootstrap样本（一般认为200次已经足够进行参数的分布、数字特征的近似估计，1000次以上认为估计精度等十分理想）。

4.Bootstrap抽样结果。从原始样本中做438次有放回的随机抽样，获得一个无序采样集，称其为一次Bootstrap抽样，计算的均值，则得到一个阈值

的估计参数，进行1000次抽样，得到估计参数的Bootstrap

抽样结果。

5.阈值确定。为待估参数计算的置信

区间（本文取），计算此置信区间的均值，

为选定的生物获得一个预测外显子编码区的最佳阈值。

为了验证基于Bootstrap的阈值确定方法，本文通过上述步骤计算了100个人和鼠类的基因序列，共438个外显子基因片段的信噪比，信噪比的计算过程通过MATLAB软件实现，Bootstrap重复抽样过程通过统计软件R实现，可得到每个片段信噪比及每个基因的平均信噪比图像。

选取a=0.1的显著性水平，进行基于Bootstrap的1000次重复抽样运算，得到90%的信噪比置信区间，并得到Bootstrap抽样结果的直方图，取区间的中点作为最优信噪比阈值见图1。

从表1和表2结果可以看出，经过1000次重复抽样，得到的信噪比阈值接近正态分布，具有很好的数学性质和统计意义，最终选定100个人和鼠类基因最优信噪比阈值为1.33，200个哺乳类动物基因最优信噪比阈值为1.55。

二、结论

目前的生物信息学所面临的三个基本挑战有：寻找基因、定位基因的外显子编码区以及预测识别基因。其中，基因预测作为生物信息学领域中的一个重要研究方向，近些年受到多数学者的青睐，取得了很多的研究成果。

对于识别不同物种类型基因编码区的阈值选取问题，本文引入Bootstrap有放回重复抽样方法，着重分析了100个人和鼠类的基因序列的识别阈值，在a=0.1的显著性水平下，取90%的信噪比置信区间中点作为最优信噪比阈值，选定100个人和鼠类基因最优信噪比阈值为1.33，200个哺乳类动物基因最优信噪比阈值为1.55。

对于频谱和信噪比区分编码区和非编码区的有效性问题，通过分析发现，其有效性是与被识别基因片段本身的长度相关的，越长的序列片段越容易被识别，而对短外显子序列的识别有效性偏低。阈值的选取也对识别有很大的影响，由100个人和鼠类基因的识别率和信噪比关系，可以认为其基因识别的阈值较小，阈值在1左右才能维持较高的DNA序列中外显子编码区的识别率。

针对基因序列的识别预测，本文结合固定长度滑动窗口上频谱曲线和基于DNA序列上“移动序列”信噪比曲线的基因识别方法，参考相关文献，提出基于傅里叶变换和固定长度滑动窗口的信噪比识别法，对6个未识别的DNA序列的编码区域的预测，比较不同窗口长度选取较好的窗宽，观察不同阈值和制定最小外显子识别长度，从多方面多角度分析预测的结果，得到6个DNA序列最终的外显子个数和区间端点。

最后在参考一些学者研究成果的基础上，简单阐述了两个比较成熟的研究方向，多重特征指标在基因特征识别和结构性突变研究。前者引入多重指标，力求最大程度避免如信噪比等单个指标分析问题的缺陷，后者从统计、数据挖掘领域引入结构性突变研究方法，应用于DNA序列的预测识别和突变问题的研究，具有一定的借鉴意义。

三、展望

任何基因预测方法都有优势、不足及应用的有效性问题。仅仅依靠单一指标或者单一的基因预测方法很难取得理想的精度预测结果，在以后的工作中，应尝试各类混合基因预测方法，提高现有预测方法的预测精度，并使预测结果更加符合生物学意义。本文在提高基于功率谱和信噪比分析的基因预测方法的预测精度方面，取得了一些成果。但对这些问题的研究还是初步的，还有很多问题值得深入探讨。为此，提出下边可进一步的研究的几个方面：

1.根据频谱3-周期性，计算功率谱与信噪比的算法及阈值的确定，一个方面是功率谱与信噪比计算公式的选择（如均值、加权均值），另一个方面是，当基因序列不是3的倍数情况下，信噪比的简便算法和更精确的取值。

2.多数评价预测基因序列的方法，都提到外显子和内含子的对立关系，从对比分析的角度去考察预测结果的精确度，但是从预测的角度同时考察外显子和内含子的特征性质的方法并不多，探求多指标或者混合的基因预测方法成为一个研究的方向。

3.基因序列（外显子、内含子）所体现的特征的充分挖掘和应用，以及基本序列本身长度等特征与编码区、非编码区的内在关系，物种类型、基因序列的特征和各种指标、预测方法之间的最优分析方法。

参考文献

[1] 蒋长锦，等.快速傅里叶变换与其C程序[M].北京：中国科学技术大学出版社，2004.

[2] 邵建峰，等.DNA序列信号的3-周期特性[J].南京工业大学学报，201234（04）.

[3] Ahamd Rushdi and Jamal Tuqan，Gene identification using the Z-curve representation[J].Acoustics， Speech and Signal Processing，1024-1027，2006 .

[4] 刘俊峰.基于功率谱分析的蛋白编码区预测方法的研究[D].成都：电子科技大学，2007.

[5] 易丽君.基于基因表达谱的数据挖掘方法研究[D].长沙：中南大学，2007.

[6] Jushan Bai and Pierre Perron.Estimating and testing linear models with multiple structural changes[J].Econometrica，Vol.66，No.1，47-78，1998.

[7] Jushan Bai and Pierre putation and analysis of multiple structural change models[J].Journal of applied Econometrics，1-22，2003

[8] 张爽等.Z曲线在基因组中应用的研究与展望[J].生物信息学，2009，7（03）：212-214.

作者简介：甄自兴（1987-），河北邢台人，研究生，首都经济贸易大学统计学院，研究方向：统计学。