1、建模实例：DNA 序列分类2000 年 6 月, 人类基因组计划中 DNA 全序列草图完成, 预计 2001 年可以完成精确的全序列图, 此后人类将拥有一本记录着自身生老病死及遗传进化的全部信息的“天书”. 这本大自然写成的“天书”是由 4 个字符 A, T, C, G 按一定顺序排成的长约 30 亿的序列, 其中没有“断句”也没有标点符号, 除了这 4 个字符表示 4 种碱基以外, 人们对它包含的“内容”知之甚少, 难以读懂. 破译这部世界上最巨量信息的 “天书”是二十一世纪最重要的任务之一. 在这个目标中, 研究 DNA 全序列具有什么结构, 由这 4 个字符排成的看似随机的序列中隐藏着什
2、么规律, 又是解读这部天书的基础, 是生物信息学（Bioinformatics）最重要的课题之一.虽然人类对这部“天书”知之甚少, 但也发现了 DNA 序列中的一些规律性和结构. 例如, 在全序列中有一些是用于编码蛋白质的序列片段, 即由这 4 个字符组成的 64 种不同的3 字符串, 其中大多数用于编码构成蛋白质的 20 种氨基酸 . 又例如, 在不用于编码蛋白质的序列片段中, A 和 T 的含量特别多些, 于是以某些碱基特别丰富作为特征去研究 DNA 序列的结构也取得了一些结果. 此外 , 利用统计的方法还发现序列的某些片段之间具有相关性, 等等. 这些发现让人们相信, DNA 序列中存在
3、着局部的和全局性的结构, 充分发掘序列的结构对理解 DNA 全序列是十分有意义的. 目前在这项研究中最普通的思想是省略序列的某些细节, 突出特征, 然后将其表示成适当的数学对象. 这种被称为粗粒化和模型化的方法往往有助于研究规律性和结构.作为研究 DNA 序列的结构的尝试, 提出以下对序列集合进行分类的问题：1）下面有 20 个已知类别的人工制造的序列（见后面）, 其中序列标号 110 为 A 类, 11-20 为 B 类. 请从中提取特征, 构造分类方法, 并用这些已知类别的序列, 衡量你的方法是否足够好. 然后用你认为满意的方法, 对另外 20 个未标明类别的人工序列（标号 2140）进行
4、分类, 把结果用序号（按从小到大的顺序）标明它们的类别（无法分类的不写入）：A 类 ； B 类 .请详细描述你的方法, 给出计算程序. 如果你部分地使用了现成的分类方法, 也要将方法名称准确注明. 这 40 个序列也放在如下地址的网页上, 用数据文件 Art-model-data 标识, 供下载：网易网址： 教育频道 在线试题；教育网： News mcm2000教育网： Nat-model-data 中给出了 182 个自然 DNA 序列, 它们都较长. 用你的分类方法对它们进行分类, 像 1）一样地给出分类结果.提示：衡量分类方法优劣的标准是分类的正确率, 构造分类方法有许多途径, 例如
5、提取序列的某些特征, 给出它们的数学表示：几何空间或向量空间的元素等, 然后再选择或构造适合这种数学表示的分类方法；又例如构造概率统计模型, 然后用统计方法分类等. Art-model-data1.aggcacggaaaaacgggaataacggaggaggacttggcacggcattacacggaggacgaggtaaaggaggcttgtctacggccggaagtgaagggggatatgaccgcttgg2.cggaggacaaacgggatggcggtattggaggtggcggactgttcggggaattattcggtttaaacgggacaaggaaggcggctggaac
6、aaccggacggtggcagcaaagga3.gggacggatacggattctggccacggacggaaaggaggacacggcggacatacacggcggcaacggacggaacggaggaaggagggcggcaatcggtacggaggcggcgga4.atggataacggaaacaaaccagacaaacttcggtagaaatacagaagcttagatgcatatgttttttaaataaaatttgtattattatggtatcataaaaaaaggttgcga5.cggctggcggacaacggactggcggattccaaaaacggaggaggcggac
7、ggaggctacaccaccgtttcggcggaaaggcggagggctggcaggaggctcattacggggag6.atggaaaattttcggaaaggcggcaggcaggaggcaaaggcggaaaggaaggaaacggcggatatttcggaagtggatattaggagggcggaataaaggaacggcggcaca7.atgggattattgaatggcggaggaagatccggaataaaatatggcggaaagaacttgttttcggaaatggaaaaaggactaggaatcggcggcaggaaggatatggaggcg8.atggccgatc
8、ggcttaggctggaaggaacaaataggcggaattaaggaaggcgttctcgcttttcgacaaggaggcggaccataggaggcggattaggaacggttatgagg9.atggcggaaaaaggaaatgtttggcatcggcgggctccggcaactggaggttcggccatggaggcgaaaatcgtgggcggcggcagcgctggccggagtttgaggagcgcg10.tggccgcggaggggcccgtcgggcgcggatttctacaagggcttcctgttaaggaggtggcatccaggcgtcgcacgctcggc
9、gcggcaggaggcacgcgggaaaaaacg11.gttagatttaacgttttttatggaatttatggaattataaatttaaaaatttatattttttaggtaagtaatccaacgtttttattactttttaaaattaaatatttatt12.gtttaattactttatcatttaatttaggttttaattttaaatttaatttaggtaagatgaatttggttttttttaaggtagttatttaattatcgttaaggaaagttaaa13.gtattacaggcagaccttatttaggttattattattatttggat
10、tttttttttttttttttttaagttaaccgaattattttctttaaagacgttacttaatgtcaatgc14.gttagtcttttttagattaaattattagattatgcagtttttttacataagaaaatttttttttcggagttcatattctaatctgtctttattaaatcttagagatatta15.gtattatatttttttatttttattattttagaatataatttgaggtatgtgtttaaaaaaaatttttttttttttttttttttttttttttttaaaatttataaatttaa16.gttat
11、ttttaaatttaattttaattttaaaatacaaaatttttactttctaaaattggtctctggatcgataatgtaaacttattgaatctatagaattacattattgat17.gtatgtctatttcacggaagaatgcaccactatatgatttgaaattatctatggctaaaaaccctcagtaaaatcaatccctaaacccttaaaaaacggcggcctatccc18.gttaattatttattccttacgggcaattaattatttattacggttttatttacaattttttttttttgtcctatagaga
12、aattacttacaaaacgttattttacatactt19.gttacattatttattattatccgttatcgataattttttacctcttttttcgctgagtttttattcttactttttttcttctttatataggatctcatttaatatcttaa20.gtatttaactctctttactttttttttcactctctacattttcatcttctaaaactgtttgatttaaacttttgtttctttaaggattttttttacttatcctctgttat21.tttagctcagtccagctagctagtttacaatttcgacacc
13、agtttcgcaccatcttaaatttcgatccgtaccgtaatttagcttagatttggatttaaaggatttagattga22.tttagtacagtagctcagtccaagaacgatgtttaccgtaacgtqacgtaccgtacgctaccgttaccggattccggaaagccgattaaggaccgatcgaaaggg 23.cgggcggatttaggccgacggggacccgggattcgggacccgaggaaattcccggattaaggtttagcttcccgggatttagggcccggatggctgggaccc24.tttagctagc
14、tactttagctatttttagtagctagccagcctttaaggctagctttagctagcattgttctttattgggacccaagttcgacttttacgatttagttttgaccgt25.gaccaaaggtgggctttagggacccgatgctttagtcgcagctggaccagttccccagggtattaggcaaaagctgacgggcaattgcaatttaggcttaggcca26.gatttactttagcatttttagctgacgttagcaagcattagctttagccaatttcgcatttgccagtttcgcagctcagtttta
15、acgcgggatctttagcttcaagctttttac 27.ggattcggatttacccggggattggcggaacgggacctttaggtcgggacccattaggagtaaatgccaaaggacgctggtttagccagtccgttaaggcttag28.tccttagatttcagttactatatttgacttacagtctttgagatttcccttacgattttgacttaaaatttagacgttagggcttatcagttatggattaatttagcttattttcga29.ggccaattccggtaggaaggtgatggcccgggggttccc
16、gggaggatttaggctgacgggccggccatttcggtttagggagggccgggacgcgttagggc30.cgctaagcagctcaagctcagtcagtcacgtttgccaagtcagtaatttgccaaagttaaccgttagctgacgctgaacgctaaacagtattagctgatgactcgta31.ttaaggacttaggctttagcagttactttagtttagttccaagctacgtttacgggaccagatgctagctagcaatttattatccgtattaggcttaccgtaggtttagcgt32.gctaccgggc
17、agtctttaacgtagctaccgtttagtttgggcccagccttgcggtgtttcggattaaattcgttgtcagtcgctctrtgggtttagtcattcccaaaagg33.cagttagctgaatcgtttagccatttgacgtaaacatgattttacgtacgtaaattttagccctgacgtttagctaggaatttatgctgacgtagcgatcgactttagcac34.cggttagggcaaaggttggatttcgacccagggggaaagcccgggacccgaacccagggctttagcgtaggctgacgctaggc
18、ttaggttggaacccggaaa35.gcggaagggcgtaggtttgggatgcttagccgtaggctagctttcgacacgatcgattcgcaccacaggataaaagttaagggaccggtaagtcgcggtagcc36.ctagctacgaacgctttaggcgcccccgggagtagtcgttaccgttagtatagcagtcgcagtcgcaattcgcaaaagtccccagctttagccccagagtcgacg37.gggatgctgacgctggttagctttaggcttagcgtagctttagggccccagtctgcaggaaatg
19、cccaaaggaggcccaccgggtagatgccasagtgcaccgt38.aacttttagggcatttccagttttacgggttattttcccagttaaactttgcaccattttacgtgttacgatttacgtataatttgaccttattttggacactttagtttgggttac39.ttagggccaagtcccgaggcaaggaattctgatccaagtccaatcacgtacagtccaagtcaccgtttgcagctaccgtttaccgtacgttgcaagtcaaatccat40.ccattagggtttatttacctgtttattt
20、tttcccgagaccttaggtttaccgtactttttaacggtttacctttgaaatttttggactagcttaccctggatttaacggccagttt本题是 2000 年网易杯全国大学生数学建模竞赛题目 A 题.作者：李清亮 王晓波 杜皓 华中农业大学 问题的简述生物学家发现 DNA 序列是由四种碱基 A, T, C, G 按一定顺序排列而成, 其中既没有“断句”, 也没有标点符号, 同时也发现 DNA 序列的某些片段具有一定的规律性和结构. 例如, 在全序列中有一些是用于编码蛋白质的序列片段, 即由这 4 个字符组成的 64 种不同的 3 字符串, 其中大多数用于编
21、码构成蛋白质的 20 种氨基酸. 而在不用于编码蛋白质的序列片段中, A 和 T 的含量特别多些. 由此人工制造两类序列 (A 类编号为 110, B 类编号为1120 ), 现在的问题是如何找出比较满意的方法来识别未知的序列 (编号为 2140 ), 并判断它们各属于哪一类. 问题的分析由于 DNA 序列全是 A, T, C, G 组成, 且长短不一, 所以我们采用提取 DNA 序列中 A, T, C, G 的百分含量这一特征来对已知 DNA 序列进行模糊分类和对未知 DNA 序列进行模糊识别.表 5-7 和表 5-8 中分别列出了已知 DNA 序列和未知 DNA 序列中含 A, T, C,
22、 G 的个数 ( A 类编号为 110, B 类编号为 1120 ; 未知 DNA 序列编号为 2140).表 5-7 已知 DNA 序列含碱基(A,T,C,G)的个数No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20A 33 30 30 47 26 39 39 31 23 20 39 36 28 33 32 40 39 32 24 22T 15 17 7 32 12 14 21 21 17 15 55 55 57 55 71 51 29 55 62 62C 19 18 24 12 26 14 11 18 23 30 5 3 11 9
23、0 9 27 13 16 19G 44 46 50 20 47 44 40 41 48 45 11 16 14 13 7 10 15 10 8 7表 5-8 未知 DNA 序列含碱基(A,T,C,G)的个数No. 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40A 31 30 18 24 26 25 24 30 15 31 27 19 30 24 25 24 22 26 29 23T 41 23 19 47 23 44 24 52 19 27 40 36 37 17 21 22 21 51 25 50C 22 25 26 2
24、2 24 24 21 17 22 26 20 25 21 24 22 32 26 20 30 23G 19 26 39 22 32 21 35 18 45 23 25 29 23 37 35 27 34 20 22 20 已知类别 DNA 序列的模糊分类 模糊聚类提取已知类别的 20 个 DNA 序列的 A, T, C, G 的百分含量构成如下矩阵：X = (xij)204，其中 xi1, xi2, xi3, xi4 分别表示第 个 DNA 系列中的 A, T, C, G 的百分含量. 采用切比雪i夫距离法建立模糊相似矩阵，然后用传递闭包法进行聚类，动态聚类图如图 5-12. 确定最佳分类设
25、, 称为总体样本的中心向量, 对应于nikk mx1),2( ),(1mx值的分类数为 r, 第 j 类的样本数为 nj , 样本记为：x 1(j), x2(j), , xn(j), 聚类中心为向量, 其中 为第 k 个特征的平均值：)()(21) mjjjx )(k，),(1xnxjikjjk ）（）（作 F 统计量：rjnijjirjj rnx12)()()( )/| 1/|上式中 为 与 的距离, 为第 j 类中样本mkkjjx12)()(| )(jx|)()(jjixxi(j)与中心 的距离.由于 F 服从自由度为 r 1, n r 的 F 分布，其分子表示类与类间的)(j距离，分母表
26、示类本身的距离，那么 F 的值越大，则说明类与类间的距离越大，即分类的结果越好. 如果 FF ( r 1, n r) ( = 0.05 )，则根据数理统计方差分析理论知道类与类之间差异是显著的，说明分类比较合理，如果满足不等式 FF ( r 1, n r) 的 F 值不止一个，则可进一步考查差 (F F) / F的大小，从较大者中找一个满意的 F 值即可.图 5-12 已知 DNA 系列的分类图 - 特征值 r - 分类数根据上述的原则，分别计算出各种分类的 F 值 ( 表 5 - 9 ) 得，当 = 0.971, F = 43.13, (FF) / F= 12.74 时，将 20 个已知 D
27、NA 序列分成如下 11 类为最佳：A1 = 1, 2, 8, A2 = 6, A3 = 3, 5, 9, A4 = 10, A5 = 7, A6 = 4, A7 = 17, A8 = 11, 12, 16, A9 = 13, 14, 18, A10 = 19, 20, A11 = 15.表 5-9 已知 DNA 序列分类的 F 值分类数 19 18 15 12 11 10 9 7 6 4 3 2F 值 54.08 55.85 44.18 44.97 43.13 36.05 33.22 15.94 13.76 18.46 25.82 33.48临界值 247.0 19.44 4.63 3.32
28、 3.14 3.02 2.95 2.52 2.96 3.24 3.59 4.41比 值 -0.78 1.87 8.54 12.55 12.74 10.94 10.26 5.33 3.65 4.70 6.19 6.59 未知 DNA 序列的模糊识别根据中的结论, 可建立标准模型库：A 1, A2, , A11. 自然我们将 A1, A2, , A6 归为 A类，将 A7, A8, , A11 归为 B 类. 但从动态聚类图 5-12 上看 , 以及模糊相似矩阵的实际意义,可将 A6 = 4, A7 = 17归为非 A, B 类.采用格贴近度公式：0(A, B)= A B + (1 AB),21计
29、算待检测 DNA 序列与标准模型库中各元素的贴近程度, 并运用“最大隶属原则”对未知类别的 DNA 序列进行模糊识别. 计算结果如下：属于 A 类的序号有：22，23，25，27，29，34，35，36，37 ;属于 B 类的序号有：21，24，26，28，31，32，33，38，40;既不属于 A 类又不属于 B 类的序号有：30，39. 模型的评价与分析 该模型成功地对 DNA 序列进行了分类，突出特点是引入模糊数学的思想，使分类效果上更具灵活性. 该模型虽然只考虑了 A, T, C, G 的百分含量四个指标的情况，但是就此模型来说它具有广泛的开放性，即该模型可以对更多的性状指标进行模糊聚类，模糊识别，特别是一些生物学指标，如(A+T) 或(C+G)的值和(A+T)/(C+G)的值，将更有利于对 DNA 序列进行分类. 该模型对于一段 DNA 片段上的碱基的相关性，即携带的遗传信息的差异等等指标考虑不足.
