研究生数学建模A题-基因问题.pdf

参赛密码 由组委会填写 由组委会填写 由组委会填写 由组委会填写 第九届第九届第九届第九届第九届第九届第九届第九届 华为杯华为杯华为杯华为杯华为杯华为杯华为杯华为杯 全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛 学学校校北京化工大学北京化工大学 参赛队号参赛队号10010002100100021001000210010002 队员姓名队员姓名 1 1 1 1 韩萱韩萱 2 2 2 2 王天宇王天宇 3 3 3 3 彭越彭越 参赛密码 由组委会填写 由组委会填写 由组委会填写 由组委会填写 第九届第九届第九届第九届第九届第九届第九届第九届 华为杯华为杯华为杯华为杯华为杯华为杯华为杯华为杯 全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛全国研究生数学建模竞赛 题 目A基因识别问题及其算法研究 摘要 本文提出了 DNA 序列频谱与信噪比的快速计算方法 给出了不同物种类型 的信噪比阈值的判定方法 提出了固定长度滑动窗口上频谱曲线结合一维离散小 波去噪的方法 进行基因识别 最后 还对以干涉指标作为较短编码序列的判别 特征 以及利用频谱发现可能的基因突变做了研究 本文的主要工作包括以下几 个方面 第一问 本文在现有 DFT 以及 FFT 算法的基础上提出了基于频数方法的快 速计算方法 在 Voss 映射的基础上 引入 Z curve 映射以及实数映射 实验结果 表明 基于 Z curve 映射的快速算法计算速度最快 且降低了随机噪声的影响 第二问 常用的信噪比 0 2R 在实际分析中效果并不理想 本文在此基础上 还引入了均值法 距离判别法以及最优化方法 利用基于 Z curve 映射的快速算 法确定了哺乳动物类等 4 种物种类型的阈值 并且定义了正确率 c A作为评价标 准 实验结果表明 最优化方法确定的阈值正确率最高 第三问 本文对固定长度滑动窗口上频谱曲线的基因识别方法得到的频谱峰 值做一维离散小波去噪 利用第二问所求阈值对 DNA 序列的编码区域进行识别 实验结果表明 一维离散小波去噪可以有效的降低 DNA 序列随机噪声的影响 提高基因识别的准确率 第四问 指出信噪比不是基因识别的唯一特征 提出了一种新的短 DNA 序 列识别的特征 干涉指标 随后本文对基因突变过程做了随机模拟 根据突变 序列的频谱数据反推出其指示序列 得到了突变后的基因 并对突变结果进行验 证 最终证明了基因突变识别过程的准确性 关键词 基因识别 小波去噪 快速算法 干涉指标 基因突变 1 目录 一 问题重述 3 二 模型假设 4 三 符号说明 4 四 问题 I 的模型建立与求解 4 4 1 1 功率谱与信噪比的快速算法的分析 4 4 1 2 功率谱与信噪比的快速算法的模型建立 5 4 1 3 功率谱与信噪比的快速算法的模型建立 5 4 1 4 功率谱与信噪比的快速算法的模型结论 7 4 2 1Z curve 映射的频谱与信噪比的模型建立 7 4 2 2Z curve 映射的频谱与信噪比的模型求解 9 4 2 3Z curve 映射的频谱与信噪比的结论 10 4 3 1 实数映射的功率谱与信噪比的快速计算公式 11 4 3 2 实数映射的功率谱与信噪比的模型求解 12 4 3 3 实数映射的功率谱与信噪比的模型结论 13 五 模型 II 的建立与求解 14 5 1 模型 II 的分析 14 5 2 模型 II 的建立 14 5 3 模型 II 的求解 15 5 4 模型 II 的结果分析及结论 19 六 模型 III 的建立与求解 20 6 1 模型 III 的分析 20 6 2 模型 III 的建立 20 6 3 模型 III 的模型求解 21 6 4 模型 III 的结果分析及结论 23 6 5 未被注释的 DNA 序列预测 23 七 问题 IV 的建立与求解 25 7 1 1 模型 IV 中问题 1 的分析 25 7 1 2 问题 1 的建立 26 2 7 1 3 问题 1 的求解 26 7 1 4 问题 1 的结果分析 27 7 2 1 第 2 问的分析 27 7 2 1 第 2 问模型的建立 27 7 2 3 第 2 问的求解 28 7 2 4 第 2 问的结果分析及结论 29 八 模型的优缺点 29 8 1 模型的优点 29 8 2 模型的缺点 29 九 参考文献 29 3 一 问题重述一 问题重述一 问题重述一 问题重述 DNA 是生物遗传信息的载体 其化学名称为脱氧核糖核酸 DNA 分子是一 种长链聚合物 DNA 序列由腺嘌呤 A 鸟嘌呤 G 胞嘧啶 C 胸腺嘧 啶 T 这四种核苷酸符号按一定的顺序连接而成 其中带有遗传讯息的 DNA 片段称为基因 在真核生物的 DNA 序列中 基因通常被划分为许多间隔的片段 其中编码蛋白质的部分 即编码序列片段 称为外显子 不编码的部分称为内含 子 DNA 序列通过遗传编码来储存信息 指导蛋白质的合成 把遗传信息准确 无误地传递到蛋白质上去并实现各种生命功能 对大量 复杂的基因序列的分析 传统生物学解决问题的方式是基于分子实验的方法 但其代价高昂 于是 人们 设法通过物理或数学的方法从大量的 DNA 序列中获取丰富的生物信息 方便学 者对 DNA 序列有更好的认识 对给定的 DNA 序列进行基因预测 是一个尚未 完全解决的问题 也是当前生物信息学中最重要的问题之一 下面我们参考题目中给出的基于频谱 3 周期性的基因预测中提出的算法 结 合附录中给出的数据以及参考文献研究如下问题 I 试设计出功率谱与信噪比的快速算法 对于很长的 DNA 序列 在计算 其功率谱或信噪比时 离散 Fourier 变换 DFT 的总体计算量仍然很大 会影响到 所设计的基因识别算法的效率 研究了针对 Voss 映射 探求功率谱与信噪比的 某种快速计算方法 并且探讨 Z curve 映射的频谱与信噪比和 Voss 映射下的频谱 与信噪比之间的关系 此外 对实数映射给出功率谱与信噪比的快速计算公式 II 确定对不同物种类型基因的阈值 对特定的基因类型的 DNA 序列 将 其信噪比R的判别阈值取为 0 2R 带有一定的主观性 经验性 对不同的基因 类型 所选取的判别阈值也许应该是不同的 针对附件和生物数据库中的数据 找出具有代表性的基因序列 并对每类基因研究其阈值确定方法和阈值结果 此 外 对按照频谱或信噪比特征将编码与非编码区间分类的有效性 以及分类识别 时所产生的分类错误作适当分析 III 基因识别算法的实现 我们的目的是要探测 预报尚未被注释的 完 整的 DNA 序列的所有基因编码序列 外显子 目前基因识别方面的多数算法 结果还不是很充分 例如前面所列举的某些基因识别算法 由于 DNA 序列随机 噪声的影响等原因 还很难 精确地 确定基因外显子区间的两个端点 对此 提出更好的解决方法并对设计出的基因识别算法的准确率做出适当评估 最后将 算法用于对附件中给出的 6 个未被注释的 DNA 序列 gene6 的编码区域的预测 IV 延展性研究 探讨以下问题 1 总结甚至独自提出一些识别基因编码 序列的其它特征指数 并对此做相关的分析 2 利用频谱或信噪比方法去发现 基因编码序列可能存在的突变 二 模型假设二 模型假设二 模型假设二 模型假设 1 给出的 DNA 序列真实可靠 能够反映出指标所代表的真实含义 4 2 假设模型中所用的所有 DNA 序列都是真核基因 3 假设 DNA 序列中的外显子不存在重叠 4 假设在问题 I II III 中所使用的 DNA 序列仅由 A G C T 这四种核苷 酸符号按一定的顺序连接而成 三 符号说明三 符号说明三 符号说明三 符号说明 A T G C 四种核苷酸的符号 S DNA序列 R 信噪比 0 R 信噪比阈值 P k 功率谱 E 总功率谱的平均值 N DNA 序列长度 四 问题四 问题四 问题四 问题 I I I I 的模型建立与求解的模型建立与求解的模型建立与求解的模型建立与求解 4 1 14 1 14 1 14 1 1 功率谱与信噪比的快速算法的分析功率谱与信噪比的快速算法的分析 当 DNA 序列较长时 DFT 运算需要耗费大量的计算时间和内存空间 很难实 时地处理问题 于是主要针对降低 DFT 运算的复杂度 加快运行时间 现有的 FFT 算法只是在 DFT 的基础上做了少许改进 仍然不能避免 Fourior 变换过程带 来的复杂的计算工作量 我们考虑如何采用一种新的快速算法使得 DNA 功率谱 与信噪比的计算大大加快 由于 DNA 序列在 N 3 处的功率谱值以及信噪比可以 由该序列上 4 种核苷酸的频数分布直接得到 那么采用核苷酸的频数分布来计算 功率谱和信噪比就可以跳过对指示序列进行 Fourier 变换 而只需统计 4 种核苷 酸在 3 种不同位置上出现的频数 并且这样的计算还有累加功能 使得信噪比计 算工作量锐减 4 1 24 1 24 1 24 1 2 功率谱与信噪比的快速算法的模型建立功率谱与信噪比的快速算法的模型建立 1 1 1 1 令 bbb zyx 表示各核苷酸 b 在 3 个密码子位置上的频数 即核苷酸 b 分别在 序列的 0 3 6 和 1 4 7 以及 2 5 8 位置上分别出现的频数 不难得到 5 222 2 2 3222 11 333 00 222 3 11 21 2 1 211 2 1 21 21 nNnjj NN jj N bbbbbb nn bbbbbbbbb b T bbbbbb b N Uu n eu n exy ez e xyzx yx zy z x xyzyX MX z 4 1 其中 T bbbb zyxX M 为上述二次型的系数矩阵 系数矩阵 M 为为半正定阵 其特征值为 1 5 1 5 0 且当 bbb zyx 式 1 的值为 0 假设 DNA 序列S表示为 0 1 2 1 SS nS nInN 其中 IA T G C A T G C分别为四种核苷酸的符号 N 为序列长度 对于任意确定的bI 令 1 0 b S nb n S nb u 0 1 2 1nN 则 DNA 编码序列S采用快速算法的功率谱峰值为 b Ib T bMX XNP 3 4 2 其中 M 为 1 式中的系数矩阵 信噪比表达式为 N MXX E NP R b Ib T b 3 4 3 其中 N N N N kP E N k 2 1 0 2 为 S的总功率谱的平均值 4 1 34 1 34 1 34 1 3 功率谱与信噪比的快速算法的模型求解功率谱与信噪比的快速算法的模型求解 采用 Matlab 2008a 软件直接编写模型来求解 数据选用的是附件中基因数据 里的 AB304259 1 中的第一段完整的外显子序列 基因位置区间为 2165 3802 其中快速傅立叶变换 FFT 算法的计算是直接调用 Matlab 中 fft 函数 将采用 快速算法的结果与原 DFT 算法和改进的 FFT 算法作比较 分析算法的优劣性 得到以下结果 6 图 4 1采用 DFT 算法 左图 和 FFT 算法 右图 得出的 DNA 序列功率谱图 方法DFTFFT快速算法 所需时间 秒 13 1430630 1068220 007586 信噪比 R16 919413919356616 919413919413916 9194139194139 峰值位置 bp 271127112710 N 3 的位置 bp 271027102710 P N 3 277142771427714 表 4 1采用 FFT 快速算法和 DCF 算法所得信噪比及运算时间的对比 4 1 44 1 44 1 44 1 4 功率谱与信噪比的快速算法的模型结论功率谱与信噪比的快速算法的模型结论 从图 4 1 中可知 DFT 算法与 FFT 算法得到的功率谱的图形几乎一致 功率 谱图形在 3 N 处附近有明显的峰值 符合碱基的 3 周期性 可以看出采用 FFT 算法并没有影响 DFT 的数据计算结果 本节中式 4 3 给出的信噪比定义 改进了基于 DFT 变换的信噪比定义式 中 序列长度 N 必须是 3 的倍数的限制 使新的快速算法适用性更广泛 使用 也更方便 实际计算中得出的信噪比的误差非常小 几乎可以认为是相同的 但 是采用快速算法的计算时间比 DFT 快约 1732 倍 大大缩短了计算时间 这与我 们理论分析中所得的结果相同 4 2 14 2 14 2 14 2 1Z curveZ curveZ curveZ curve 映射的频谱与信噪比的模型建立映射的频谱与信噪比的模型建立 我们基于上一问的快速算法来建立 Z curve 映射的模型 定义三个序列 x n y n z n 2 2 2 nn nn nn x nAGn y nACn z nATn 4 4 令 1 0 x 1 0y 和 1 0z 以 及 1 x nx nx n 1 y ny ny n 和 7 1 z nz nz n 得到 Z curve 映射 1 2 3 1111 1111 1111 AA CC GG TT unun x na unun y na unun z na u nu n 2 1 1 0 A nk N j C N n G T un un x naX kx n e un u n 2 2 22 33 1 2 22 33 1 1 111111 1 3 11 21 2 1 211 2 1 21 21 jj TTT ACGTACGTACGT jj T N Xaxxxxeyyyyezzzz axeyez r r s tsR MR t 4 5 同理 22 222222 2 11 21 2 1 211 2 3 1 21 21 T r N Yr s tsR MR t 32 333333 3 11 21 2 1 211 2 3 1 21 21 T r N Zr s tsR MR t 4 6 其中 TTT ACGTACGTACGT xxxxxyyyyyzzzzz iiiiii ra x sa y ta z 1 2 3 T iiii Rr s ti 采用快速算法的功率谱峰值为 3 1 3 T ii i P NR MR 4 7 信噪比表达式为 3 1 3 3 T ii i R MR P N R EN 4 8 8 其中 1 2 0 3 3 N k P k N EN NN 是 Z curve 映射的平均功率谱 采用 DFT 算法的 Z curve 映射的总功率谱为 222 Z P kX kY kZ k 4 9 其中 X k Y k 和 Z k 分别表示数字序列 x n y n 和 z n 的快速傅立叶变 换 采用 DFT 算法的 Z curve 映射的信噪比为 222 333 3 Z Z NNN N XYZ P R EE 4 10 其中 1 0 N k P k E N 是 Z curve 映射的平均功率谱 4 2 24 2 24 2 24 2 2Z curveZ curveZ curveZ curve 映射的频谱与信噪比的模型求解映射的频谱与信噪比的模型求解 使用 Matlab 编程对上述模型进行编程 数据的选取与 4 1 3 节相同 选取 AB304259 1 中的第一段完整的外显子序列 得到模型求解的图像与结果如下 图 4 2采用 Z curve 映射 图左 与 Voss 映射 图右 的功率谱图像对比 FFT 9 图 4 3 在 FFT 算法下的 Z curve 与 Voss 映射功率谱的散点图 其中 Z curve 蓝色 的点是 Voss 变换 红色的 是 Z curve 变换 方法Voss FFT Z curve DFT Z curve FFT Z curve 快速算法 所需时间0 10682210 6147950 0472760 003076 信噪比 R16 9194139122 5592185622 5592185622 55921856 P N 3 27714110856110856110856 峰值位置2711271127112710 N 3 的位置2710271027102710 表 4 2采用 Z curve 映射的三种算法的统计数据 4 2 34 2 34 2 34 2 3Z curveZ curveZ curveZ curve 映射的频谱与信噪比的结论映射的频谱与信噪比的结论 从图 4 2 4 3 中可知在同时使用 FFT 算法的前提下 通过 Z curve 映射得到 的图像的纵坐标轴取值要大于采用 Voss 映射 可知 Z curve 映射下的功率谱的值 要大 但是采用两种映射的图像的形状走势几乎一致 在 3 N 处有明显的峰值 符合碱基的 3 周期性 同时可以验证我们在 4 1 节中使用的快速算法同样适用于 10 不同的映射方法 证明了新的快速算法适用性广泛 实际计算中 采用 Z curve 的三种算法之间的的信噪比和功率谱的误差相当小 可近似为 0 并且使用 Z curve 映射的外显子序列的功率谱曲线在频率 3kN 处 具有较大的频谱峰 值 相对 Voss 映射而言 大大降低了随机噪声的影响 在处理时间上可明显看 出采用 Z curve 映射的计算时间远远少于 Voss 映射 在采用相同算法下的运行时 间要节省约 2 5 倍 又缩短了计算时间 提高了计算效率 4 3 14 3 14 3 14 3 1 实数映射的功率谱与信噪比的快速计算公式实数映射的功率谱与信噪比的快速计算公式 任意的DNA序列 0 1 2 1 SS nS nInN 其中 A T G C分别为四种核苷酸的符号 令 0 0 1 2 3 AA CC GG TT un S nAun un S nCun T n un S nGun un S nTun 0 1 2 1nN 由此得到的映射序列中酮基 G T 大于氨基 A C 2 2 22 33 0 3 11 21 2 1 211 2 1 21 21 jj TTT ACGTACGTACGT T N Txxxxeyyyyezzzz r r s tsR MR t 4 11 其中 其中 TTT ACGTACGTACGT xxxxxyyyyyzzzzz 000 ra x sa y ta z TRr s t 实数映射采用快速算法模型功率谱峰值 3 T P NR MR 4 12 信噪比表达式为 3 4 T P NR MR R EN 4 11 其中 1 0 4 N k P k E N 是实数映射的平均功率谱 实数映射采用 FFT 算法模型如下 11 1 1 0 2 1 0 NkenSkU N nk j N n S的功率谱序列 P k 2 kUkP 信噪比为 3 N P R E 其中 1 0 N k P k E N 为S的总功率谱的平均值 4 3 24 3 24 3 24 3 2 实数映射的功率谱与信噪比的模型求解实数映射的功率谱与信噪比的模型求解 选取与 4 1 3 节相同的数据 利用 Matlab 编程得到模型求解的图像与结果如 下 图 4 4 实数映射的功率谱图 FFT 12 信噪比 R2 4394 计算时间 秒 FFT 算法0 861886 快速算法0 059543 表 4 3 采用实数映射的统计数据 4 3 34 3 34 3 34 3 3 实数映射的功率谱与信噪比的模型结论实数映射的功率谱与信噪比的模型结论 从图 4 中可知采用实数映射得出的功率谱图噪声较大 但在 3 N 处仍然可以 看到比较明显的峰值 这是由于实数映射并非线性变换 导致随机噪声比较大 且计算量大大增加 正因为如此 较 Z curve 映射的模型 实数映射模型计算所 需的时间也大大增加 计算效率降低 五 模型五 模型五 模型五 模型 II II II II 的建立与求解的建立与求解的建立与求解的建立与求解 5 15 15 15 1 模型模型 II II II II 的分析的分析 现有的外显子的信噪比2R 是它一个普遍的特征 因此大部分文献以阈值 取 2 为指标来区别和区分外显子和内含子 但是在实际分析中 图 5 1 外显子与内含子的信噪比统计 可以看出不同基因的信噪比均值有较大差异 且外显子的信噪比的标准差远 大于内含子 但信噪比大于 2 的仅为 56 86 本文考虑对不同的基因分别确定 一个阈值 0 R 然后将现有的信噪比R与 0 R做比较 使 0 RR 作为外显子的判别 0 RR 作为内含子的判别 最后将判别的结果作正确率的分析 考虑 0 R是否是 13 一个稳健的估计值 5 25 25 25 2 模型模型 II II II II 的建立的建立 本文分别采用四种不同的方法来确定阈值 0 R 同时基因序列的处理方法采 取的是 4 2 节中提出的 Z curve 映射的快速算法 一 均值法 对于已有的基因数据 设所有外显子的信噪比均值为 1 m 所有内含子的信 噪比均值为 2 m 一种简单确定阈值的方法是将它们作算术平均 即令 01 0 2 mm R 5 1 二 直接定义阈值 0 2R 三 距离判别法 距离判别方法是模式识别理论中的一种重要分类方法 也是生物信息学中 广泛应用的分类方法 距离判别法的基本思想是根据已知分类的数据 分别计算 各类的中心 即分组 的均值 判别准则是对任给的一次观测 若它与第i类的 重心距离最近 就认为它来自第i类 由此达到分类的目的 3 四 最优化算法 1 设所有外显子的信噪比值组成的集合为 1 S 所有内含子信噪比值组成的集合 为 2 S 欲寻求的最优分类阈值为 0 R 设 1 2 12 ij RS RS 求解阈值 0 R的优化模 型为 1 2 00 maxsgnsgn ij ij RRRR 5 2 其中 a b 0 aRb 为信噪比值域区间 即在基因外显子 内含子信 噪比样本集上 优化模型求的是判别正确率达到最大的阈值解 0 R 为了评价判别的结果 我们设计了如下指标 nPPN STTF 5 3 PNNP STTF 5 4 其中 n S为敏感性指标 p S是专一性指标 P T表示被正确判为外显子的个数 N T 表示被正确判为内含子的个数 N F表示被错误递判为内含子的个数 P F表示错 误的判为外显子的个数 14 定义阈值判别的总正确率 c A为 2 np c SS A 5 5 5 35 35 35 3 模型模型 II II II II 的求解的求解 数据采用的是题目附件中所给出的 Genes100 mat 与 genes200 mat 中所含的 Homo sapiens Mus musculus Rattus norvegicus 和哺乳动物类的总计 1265 条外 显子基因序列和 1563 条内显子基因序列样本 其中由于 Mus musculus 和哺乳动 物类样本数据量较大 我们选取样本数据中的 70 作为训练样本 剩余的 30 作为测试样本来验证阈值的有效性和正确率 使用软件 Matlab2008a 来对以上模 型进行编程 得到以下结果 基因种类 均值平均法 0 R 2距离判别法最优化方法 1 1 1 1 0 0 0 0 R R R R c c c c A A A A 2 2 2 2 0 0 0 0 R R R R c c c c A A A A 3 3 3 3 0 0 0 0 R R R R c c c c A A A A 4 4 4 4 0 0 0 0 R R R R c c c c A A A A Homo sapiens 4 47750 625020 73742 550 71591 29000 7917 Rattus norvegicus 4 65030 568720 73432 6450 65051 16300 7788 Mus musculus 训练 3 73680 614820 75632 100 75731 35800 7905 Mus musculus 预测 3 73680 643520 76982 100 75421 35800 8737 哺乳动物 训练 3 73040 628820 77702 0850 76651 37800 8097 哺乳动物 预测 3 73040 643520 78012 0850 76361 37800 8006 表 5 1采取四种不同的方法预测的阈值和正确率 15 图 5 2人类基因序列判别正确率曲线 图 5 3 小鼠基因序列判别正确率曲线 16 图 5 4家鼠基因序列判别正确率曲线 图 5 5 哺乳动物基因序列判别正确率曲线 17 5 45 45 45 4 模型模型 II II II II 的结果分析及结论的结果分析及结论 从表5 1中可以得出 4种基因序列在采用最优化算法得出的阈值判别的正确 率最高 在采用均值算法的时候正确率最低 在阈值的选取过程中 可以看出最 优化方法对阈值的判别最为准确 也可以看出通过最优化方法得出的各种基因类 型的阈值都在2以下 也就是说如果采用简单地取 0 2R 将把部分外显子误判 为内含子 使基因数据发生误差 采取四种方法在进行DNA序列信噪比的统计时 正确率并不能达到90 以上 原因是部分外显子的信噪比并不显著 同样地也存在信噪比较大的内含子 产生 这一现象的原因至少有下面两点 1 不同基因类型对信噪比的影响 对于不同类 型的基因 其外显子的信噪比也有着显著的差别 2 同一类型中不同外显子信 噪比的独特性与差异性 由于每段外显子都可以计算得到自己独特的信噪比 采 用统一的阈值作为外显子与内含子的判断标准必然会造成一定的错误 本文所采 用的最优化方法已使得正确率达到最高 也就是最低的误判概率 通过我们最终得出的数据结果 用于识别基因外显子的信噪比 不能仅凭经 验或者简单统计就可以获得 事实上 基因序列统计样本的大小 基因外显子序 列的长短 甚至同种生物的不同基因类型 对信噪比阈值的确定都有很大的影响 六 模型六 模型六 模型六 模型 IIIIIIIIIIII 的建立与求解的建立与求解的建立与求解的建立与求解 6 16 16 16 1 模型模型 IIIIIIIIIIII 的分析的分析 如上节所述 由于在第 5 3 节中得出的阈值 0 R的正确率 c A是小于 1 的数 所以在基因识别过程中存在内含子和外显子误判的两种情况 存在外显子真实的 信噪比 0 RR 这就是说在基因序列上存在着信息 重叠的情况 而在第 5 3 节我们已经得出了不同物种基因序列的最优阈值 0 R 本 文基于模型 II 确定最优阈值的方法来实现基因识别算法 获取的基因频谱的信号是含噪声信号 由于受噪声影响 将对确定外显子区 间造成影响 故在对信号进一步分析之前 需要将真实信号提取出来 目前 人 们已根据噪声的统计特征和频谱分布规律 开发出了多种多样的信号去噪方法 其中最为直观的一种方法是 根据噪声能连一般集中于高频 而信号频谱分布于 一个有限区间的特点 用傅里叶变换将含噪信号变换到频域 然后采用低通滤波 器进行滤波 6 26 26 26 2 模型模型 IIIIIIIIIIII 的建立的建立 本文采用基于固定长度滑动窗口上频谱曲线的基因识别方法结合一维离散 18 小波去噪方法进行基因识别 步骤一 固定长度滑动窗口求 3 M p n值 对一个DNA序列S作 Z curve 映射 0 1 2 1nN 取长度M作为固定 窗口长度 对任意n 01nN 在以n为中心的长度为M的序列片段 n 1 2 M n 1 2 M 上 当n接近序列的两端时 窗口实际有效长度可能会小于M 采用 快速算法 求出它在 M 3 处的功率谱峰值 3 M p n 即 3 1 3 T ii i M PR MR 6 1 步骤二 一维离散小波去噪 5 使用 MATLAB 小波 GUI 利用步骤一得到的功率谱峰值 3 M p n 0 1 2 1nN 输出含噪声信号 并且进行一维离散小波去噪 得到去噪后 的功率谱峰值 3 M p n 以及降噪信号 步骤三 结合阈值 预测外显子区间 利用去噪后的功率谱峰值 3 M p n 结合模型 2 中阈值判定的最优化方法 得到最优阈值 并预测外显子区间 6 36 36 36 3 模型模型 IIIIIIIIIIII 的模型求解的模型求解 为了验证去噪过程的有效性 我们选择 gnens200 mat 中的第 199 组基因 RNGMTG 进行判断 并对采用了去噪方法和未采用去噪方法的基因识别结果进 行比对 为此我们定义准确率 为评价指标 正确判断为外显子的基因长度 正确判断为外显子的基因长度 误判为内含子的基因长度 误判为外显子的基因长度 取 5 3 节所求该组基因所对应的哺乳动物类最优化阈值 R0 1 378 窗宽 M 为 97 利用 Matlab 编程得到结果如下 19 图 6 1 未采用去噪方法的基因频谱及识别结果 图 6 2 采用去噪方法的基因频谱及识别结果 上述两图中 红色粗线为真实外显子区间 绿色细线为基因识别模型处理结果 20 外显子123456 区间1 94430 6281291 13871432 17561816 21102146 2503 外显子789101112 区间3088 32773730 42854327 44744498 47834843 49425170 5259 表 6 1 外显子区间判断结果 采用去噪方法 未去噪方法识别准确率23 1 去噪方法识别准确率32 2 表 6 2 两种识别方法的准确率 6 46 46 46 4 模型模型 IIIIIIIIIIII 的结果分析及结论的结果分析及结论 根据 6 3 节结果显示 采用了去噪方法后的基因识别结果有所改善 基因识 别准确率由 23 1 提高至 32 2 在图中无论是频谱图像还是识别结果也可以明 显的看出去噪效果 将基于固定长度滑动窗口上频谱曲线的基因识别方法与一维 离散小波去噪方法相结合 可以在一定程度上改善因识别的准确率问题 但是 这种结合方法的准确率仍然有限 究其原因 一方面是 5 4 节所述的阈值确定方 法的正确度 c A小于 1 每段外显子具有独特性与差异性 另一方面是有些外显子 如外显子长度较短 的 3 周期特性并不明显 在进一步的研究中 可以结合生 物学基因本身的特性加以考虑 6 56 56 56 5 未被注释的未被注释的 DNADNADNADNA 序列预测序列预测 下面按照题目要求 对 genes6 mat 文件中的未被注释的 DNA 序列进行编码 区域的预测 由于时间原因 选择了第一组 第二组和第四组基因进行了预测 结果如下 外显子1234 区间1160 159521691 18832420 25342672 2933 外显子5678 区间3188 33623635 38844460 47005060 5155 表 6 3 第一组外显子区间判断结果 采用去噪方法 21 图 6 3 第一组基因频谱与预测结果 左图未去噪方法 右图去噪方法 外显子123 区间1254 14223387 35613885 4038 外显子45 区间4116 43505574 5984 表 6 4 第二组外显子区间判断结果 采用去噪方法 图 6 4 第二组基因频谱与预测结果 左图未去噪方法 右图去噪方法 外显子123 区间1499 16912558 35034736 4952 外显子45 区间5030 53035465 5689 表 6 5 第四组外显子区间判断结果 采用去噪方法 22 图 6 5 第四组基因频谱与预测结果 左图未去噪方法 右图去噪方法 七 问题七 问题七 问题七 问题 IVIVIVIV 的建立与求解的建立与求解的建立与求解的建立与求解 7 1 17 1 17 1 17 1 1 模型模型 IVIVIVIV 中问题 中问题 1 1 1 1 的分析 的分析 采用频谱或信噪比这样单一的判别特征 也许是影响 限制基因识别正确率 的一个重要原因 在进行 DNA 序列信噪比的统计时 注意到部分外显子的信噪比并不显著 同样的也存在信噪比比较大的内含子 这是造成误判的根源 产生这一现象的原 因有以下两点 1 不同基因类型对信噪比的影响 对于不同类型的基因 其外显子的信噪比也有着显著的差别 2 序列长度对信噪比的影响 另外统计了 genes200 mat 中 290 个外显子和 112 内含子片段的信噪比 按长 度将外显子分成 3 类 长度 100bp 为短外显子 长度 100bp 为长外显子 内 含子也做同样处理 用确定阈值的优化方法得到判别正确率统计如表 3 所示 由 表 3 可见 哺乳动物类基因序列随着外显子长度的增加 外显子序列的信噪比和 判别正确率也随之增加 短外显子在哺乳动物类外显子中占的比率最大 判别正 确率也最低 在短外显子所占比率较大的生物基因中 就易于出现外显子判别正 确率偏低的现象 非均匀指数 HI 也是一个十分重要的判别特征 由偏好模分析法和关联函数 谱分析法 都表明外显子中存在三重读码框架 而内含子及其他非编码区中 一 般不存在读码框架 根据这一性质 我们引入了反映密码子 3 个位点上核苷酸分 布不均匀参数 HI 2 2 2 2 34343434 1 1 1 1 llllllll l l l l bbbbbbbb b b b b lblblblb NNNNNNNNNNNNNNNN HINHINHINHIN NNNNNNNN 7 1 这里 1 2 3 41 2 3 41 2 3 41 2 3 4 b b b b NbNbNbNb 表示序列中的 4 种碱基数 b b b b b b b b NNNNNNNN 为序列长度 23 1 2 31 2 31 2 31 2 3 l l l l NlNlNlNl 为 3 个子序列的长度 3 3 3 3 l l l l N N N N N N N N l l l l b b b b N N N N为第l l l l个子序列中第b b b b个 种碱基数 7 1 27 1 27 1 27 1 2 问题 问题 1 1 1 1 的建立 的建立 由于一些外显子和内含子序列较短 统计特征不明显 可将n n n n个序列 外显 子按框架结构 拼接 定义其非均匀指标的模型值和观测值分别是 mod12mod12mod12mod12 n n n n obsobsobsobs HIHIHIHIHIHIHIHIHIHIHIHIHIHIHIHI HIHIHIHIHIHIHIHI 上式中 1 2 1 2 1 2 1 2 i i i i HIinHIinHIinHIin 表示被拼接的外显子或内含子的 HI 值 而HIHIHIHI是指多个 外显子或内含子被拼接后 称为一段多外显子 Multi Exon ME 或多内含子 Multi Intron MI 序列的 HI 值 由此我们定义n n n n个片段的干涉指标T T T T为 1 1 1 1 n n n n i i i i modimodimodimodi obsobsobsobs HIHIHIHI HIHIHIHI T T T T HIHIHIHIHIHIHIHI 根据编码区和非编码区的统计特征 理论上内含子的干涉指标T T T T应该为n n n n 而外 显子的T T T T值为 1 7 1 37 1 37 1 37 1 3 问题 问题 1 1 1 1 的求解 的求解 采用 genes200 mat 文件中 bp 值小于 100 的外显子和内含子序列 建立模型 使用 Matlab 求解干涉指标T T T T 统计结果显示外显子的T T T T值为 2 16 内含子的T T T T值 为 48 23 7 1 47 1 47 1 47 1 4 问题 问题 1 1 1 1 的结果分析 的结果分析 当序列片段数n n n n足够大 ME 和 MI 的T T T T值分布的差别将会显示出来 可以 用来区别编码区和非编码区 7 2 17 2 17 2 17 2 1 第 第 2 2 2 2 问的分析 问的分析 基因突变具有随机性 低频性 可逆性 可遗传性等特点 一般发生的概率 在 1 2 左右 本文将采用随机模拟的方法模拟基因突变的发生过程 得到基 因突变后 DNA 序列的频谱图像 与未发生突变的 DNA 序列频谱进行对比 将 发生了基因突变的频谱数据反推出其指示序列 得到突变后的基因 验证突变结 果以表明该过程的准确性 由于时间原因 只考虑单核苷酸发生替换变异的情况 24 7 2 27 2 27 2 27 2 2 第 第 2 2 2 2 问模型的建立 问模型的建立 本节采用 Z curve 的 FFT 算法进行计算 得到频谱图像是建立在 Z curve 映 射的基础上 故由发生了突变后的 Z curve 映射进行傅立叶变换后的 X k Y k 和 Z k 作为频谱的替代值 将三序列采用 IFFT 算法进行逆 Fourior 变换 得到 x n y n 和 z n 回忆 Z curve 映射 1111 1111 1111 A C G T un x n un y n un z n un 由于上式的系数矩阵不是方阵 该式不能反解 但考虑由于是核苷酸替换过程 A un C un G un 与 T un 之和不变 故得到反解过程如下 1 1 11 1 11 1 1 1 1 11 1111 A C G T unx