1、DNA 序列分类模型DNA 序列分类模型摘要本文分析了已知类别的人工 DNA 序列的特征，建立了聚类分析延拓模型和马尔可夫模型，分别对未知类别的人工 DNA 序列和自然序列进行分类，根据分类效果选出了较优模型。首先对数据进行预处理，得到人工 DNA 序列的单个碱基丰度和不同碱基丰度之比等特征量，进而分析 A、B 两类的差异，得到合适的特征判定条件对未知类别的 DNA 序列进行分类。计算人工 DNA 序列的特征量，给出各序列的统计数据。其次用聚类分析延拓模型进行分类。用 A、B 两类具有明显差异的特征作为样品特征变量，得到欧式空间中表征编号 1-20 人工 DNA 序列的特征向量，计算两两之间的
2、 Lance 和 Williams 距离进行相似性度量，逐步选择相似性较大的归为一类，同时不断更新类内的标准比较特征向量，对聚类方法进行延拓，最终得到类内差异小、类间差异大的 A、B 两类，建立了聚类分析延拓模型。再对选取的特征变量进行改进，提高模型的分类效果。最后，借助均值、方差和相关系数等参数对改进模型的分类效果进行分析。再次用马尔可夫模型进行分类。将 DNA 序列看成是马尔可夫链，求出编号1-10 和 11-20 人工 DNA 序列在已知当前碱基种类的条件下，下一个碱基出现任一种的概率，结果存入概率转移矩阵 1 和 2，再利用矩阵 1 和 2 分别求出编号1-20 中任一条 DNA 序列
3、出现的概率，选择较大的一个作为该 DNA 序列的分类，建立马尔可夫模型。再进行与聚类分析延拓模型类似的改进和检验工作，然后对编号 21-40 人工 DNA 序列和 182 条自然序列进行分类，得到最终结果。最后，用层次分析法综合评价模型一与模型二，选择聚类分析延拓模型作为最终模型，其分类结果作为最终结果，具体如下：编号 21-40 人工 DNA 序列中属于 A 类的样品编号为：22，23，25 ， 27，29， 30，34，35，36，37，39 ；属于 B 类的样品编号为：21，24，26 ， 28，31，32，33，38，40。182 条自然序列中，属于 B 类的样品编号为：7，10， 1
4、2，22 ，23，24， 26，28，30，34，43，48，50，54，57，65，75，76，80， 84，85 ，86，92， 98，103，107，110，114，116，119，121，122，123，127，128，129，130，131，137，138，140， 142，143，144，146，151，156， 159， 161，162， 163，166，168，170，173，174，175，179，180，181，182；其余为 A 类。关键词 DNA 序列分类 聚类分析延拓法 Lance 和 Williams 距离 马尔可夫法 一、问题重述1.1 题目背景（1）2000 年
5、 6 月，人类基因组计划中 DNA 全序列草图完成，预计 2001 年可以完成精确的全序列图，此后人类将拥有一本记录着自身生老病死及遗传进化的全部信息的“天书” 。 （2）这本 “天书”是由 4 个字符 A，T，C，G 按一定顺序排成的无间隔的长约 30 亿的序列，除了这 4 个字符表示 4 种碱基以外，人们对它包含的“内容”知之甚少。因此，破译这部世界上最巨量信息的“天书”是二十一世纪最重要的任务之一。（3）为解读这部“天书” ，首先要研究 DNA 全序列具有什么结构，以及由这4 个字符排成的看似随机的序列中隐藏着什么规律，这也是生物信息学最重要的课题。1.2 题目信息（1）DNA 序列分为
6、编码区与非编码区。编码区是用于编码蛋白质的序列片段，即由这 4 个字符组成的 64 种不同的 3 字符串，其中大多数用于编码构成蛋白质的 20 种氨基酸。（2）在不用于编码蛋白质的序列片段中，A 和 T 的含量特别多些，于是以某些碱基特别丰富作为特征去研究 DNA 序列的结构也取得了一些结果。（3）利用统计的方法还发现序列的某些片段之间具有相关性。这些发现说明 DNA 序列中存在着局部的和全局性的结构，充分发掘序列的结构对理解 DNA 全序列有十分重要的意义。目前在这项研究中最普通的思想是省略序列的某些细节，突出特征，然后将其表示成适当的数学对象。 1.3 题目要求（1）有 20 个已知类别的
7、人工制造的 DNA 序列（见附件 1） ，其中序列标号110 为 A 类，11-20 为 B 类。从中提取特征，构造分类方法，并用这些已知类别的序列，衡量所选分类方法是否足够好。（2）用（1）中的分类方法对另外 20 个未标明类别的人工序列（见附件1，标号 2140）进行分类，根据分类效果对方法不断完善，将得到的最终结果用序号（按从小到大的顺序）标明它们的类别（A 类或 B 类，无法分类的不写入） 。 要求详细描述所选的分类方法，给出计算程序。若论文中部分地使用了现成的分类方法，应将方法名称准确注明。 （3）已知 182 个自然 DNA 序列（见附件 2） ，它们都较长。同样用以上所选的分类方
8、法对它们进行分类，并根据分类效果对方法不断完善，像（2）中一样给出最终的分类结果。 二、 名词解释1.编码区与非编码区：编码区是指 DNA 上编码蛋白质的序列片段，而非编码区不用于编码蛋白质。2.聚类分析：由已知数据，计算各个观察个体或变量之间亲疏关系的统计量。再根据某种准则（最短距离法、最长距离法、中间距离法、重心法等） ，使同一类内的差别较小，而类与类之间的差别较大，最终将观察个体或变量分为若干类的分类方法。其中，对样品所作的分类为Q-型聚类，对变量所作的分类为R-型聚类。3.相似性度量：对数值型数据而言，两个个体的相似度是指它们在欧氏空间中互相邻近的程度；而对分类型数据而言，两个个体的相
9、似度与它们取值相同的属性的个数有关。4.样品：每个观察个体即每条DNA序列为一个样品。5.样品变量：每个样品所具有的不同特征用不同的变量来表示，变量数等于特征数。6.碱基丰度：每条 DNA 序列中碱基 A、G、C 或 T 出现的频率。三、 问题分析DNA 序列分类问题要求在对 DNA 序列的一些规律和结构有所了解的基础上，从 20 个已知类别的人工制造的 DNA 序列中提取特征，构造分类方法，并用所选择的分类方法对其余未知类别的 20 个人工制造的 DNA 序列以及 182 个自然 DNA序列进行分类。3.1 建模目标的分析DNA 序列分类是一个复杂的统计分析问题，数据量大，影响因素多，无法直
10、接从 20 条已知类别的人工制造的 DNA 序列中提取出所有的有效特征，因此有必要对这 20 条 DNA 序列进行预处理。观察并分析数据预处理结果，归纳总结出 A 类和 B 类的有效特征，将其表示成适当的数学对象，并选择适当的分类方法，建立普遍意义下数学模型，再用得到的模型对其余未知类别的 20 个人工制造的 DNA 序列以及 182 个自然 DNA序列进行分类。由题意，建立的数学模型应该保证分类结果具有以下特点：（1）类别间差异尽量大；（2）类别内差异尽量小；（3）样品能够尽可能的落入 A、B 范围，且只能落入其中的一个。3.2 建模及求解方向1.分析已知类别的 DNA 序列 1-20 的结
11、构，提取出相应的特征。主要的特征有：碱基的丰度、碱基或碱基序列的重复出现情况、碱基或碱基序列之间的相邻情况、不同碱基的丰度之比（如碱基 A 与碱基 T 的丰度之比）等。2. 根据提取出的特征，选用合适的分类方法。对数据进行预处理后，尝试以下方法建立模型：（1）根据聚类分析法，建立模型一。由题意，DNA 序列分类属于对样品所做的分类，为 Q-型聚类。首先引入样品变量，例如可选择碱基 T 的丰度、碱基 G 的丰度、碱基 T 与碱基 G 的丰度之比、碱基 A 与碱基 T 的丰度之比等。由已知数据，计算出每条已知类别的人工制造的 DNA 序列的各个样品变量值，存入向量中。根据相似性度量原理，计算 20
12、 个样品两两之间的 Lance 和 Williams 距离，选择相距最远的两个样品（假设为样品 3 和样品 16）分别作为 A 类和 B 类，再分别以样品 3 和样品 16 为标准点，通过分别计算样品 3 和样品 16 与其余 18 个样品之间的 Lance 和 Williams 距离，找出与其相距最近的一个样品（假设为样品 1 和样品 18）归为一类。此时，新的标准点变为样品 1 与样品 3 的中点、样品 16 与样品 18 的中点。然后再以新的标准点为基准，分别找出与其相距最近的一个样品归为一类。逐步进行下去，直至 20 个样品被明显分成 A、B 两类。（2）根据马尔可夫法，建立模型二。以
13、单个碱基为单位，分别统计编号 1-10 和编号 11-20 人工制造的已知类别的 DNA 序列中 4 种碱基出现的次数，再以相邻的两个碱基为单位（共 16 种组合情况） ，分别统计编号 1-10 和编号 11-20 的 DNA 序列中 16 种碱基对出现的次数。为满足大样本需求，将 A 类和 B 类中的 10 条 DNA 序列组合起来看作两个大样品，单个碱基或碱基对出现（不包括上一条链的末尾碱基与下一条链的初始碱基组合的情况）的次数为 10 条序列之和。由条件概率的思想，分别求出 A 类和 B 类大样品中在已知当前碱基种类（可以为 A、G、C、T 中任何一个）的条件下，下一个碱基分别为 A、G
14、、C、T的概率，存入两个矩阵 1 和 2 中。对于任何一条给定的 DNA 序列，可将其看作一个已经发生的事件，说明该事件发生的概率比较大。用矩阵 1 和矩阵 2 分别求出这一事件发生（即形成当前 DNA 序列）的概率，若用矩阵 1 算出该编号的DNA 序列出现的概率较大，则该编号的 DNA 序列属于 A 类，否则属于 B 类。3.模型的初步检验与改进。用编号 1-20 已知类别的序列，分别衡量模型一与模型二中所选方法是否足够好，不断改进，尽可能使 1-20 号 DNA 序列在所选分类方法下，所得结果与已知分类完全一致。改进时，对于聚类分析法，可以尝试改变样品变量的个数或者改变样品变量的组合方式
15、；对于马尔可夫法，可以尝试引进中间变量，运用隐马尔可夫法求解。4.模型的进一步检验与完善。（1）用以上的得到的两种分类方法对编号 20-40 未知类别的人工序列、182个自然序列进行分类。（2）通过计算样品方差、均值等比较两种分类方法得到的分类结果与建模目标类别间差异尽量大、类别内差异尽量小、样品能够尽可能的落入 A、B范围，且只能落入其中的一个的接近程度。（3）选择更接近建模目标的一种分类方法作为最终的分类方法，其分类结果即为最终结果。四、基本假设1.假设所给的 DNA 序列片段中没有断句和标点符号。2.假设具有特殊碱基的 DNA 序列中，特殊碱基可以剔除，其影响可以忽略。3.较长的 182
16、 个自然序列与已知类别的 20 个样本序列具有共同的特征。4.假设给定的DNA序列均是从全序列中随机截取出来的,无法确定序列的起始位, 无法从序列中辨认出氨基酸，所以，在对DNA 序列分类时,从碱基层次上进行分类, 而不是从氨基酸层次上分类。五、定义与符号说明：各个样品中碱基 出现的数量，i 为 A、T 、C 或 Gin：第 i 个样品的总碱基数目iN：各个样品中碱基 的丰度，i 为 A、T、C 或 GiF：各个样品的第 i 个特征变量ix：各个样品中碱基 i 和碱基 j 的比值，i，j 为 A、T、C 或 Gijf：第 i 个样品的特征向量iY：向量 和向量 间的 Lance 和 Willi
17、ams 距离ijdijY：特征向量的分量个数，即向量的维数p：特征向量的第 k 个分量k：样品的个数n：特征向量 i 的第 k 个分量ikx：不同向量代表的 维空间中任意两点间 Lance 和 Williams 距离的最大值madp：不同向量代表的 维空间中任意两点间 Lance 和 Williams 距离的最小值in：聚类分析中 i 类的标准向量，i 为 A 或 Biy六、数据预处理1.A 类和 B 类样品单个碱基丰度的计算用 maTlab 编写程序（见附件 3） ，分别求出 20 条已知类别的人工制造的DNA 序列中，4 种碱基的丰度，绘出散点图如下：图 6.1.1 单个碱基丰度比较图分析
18、上图可得， A 类和 B 类 DNA 序列中碱基 T 和碱基 G 的丰度有明显差异，而碱基 A 和碱基 C 的丰度则比较接近。2. A 类和 B 类样品不同碱基丰度之比的计算用 matlab 编写程序（见附件 4） ，分别求出 20 条已知类别的人工制造的DNA 序列中，不同碱基的丰度之比，包括 、 、 、 、 、 ，绘出TAfCGAfCTGfC散点图如下：图 6.1.2 不同碱基丰度之比的比较图分析上图可得， A 类和 B 类 DNA 序列中，碱基 T 与碱基 A 的丰度之比、碱基 G 与碱基 A 的丰度之比、碱基 C 与碱基 T 的丰度之比、碱基 G 与碱基 T 的丰度之比有明显差异，而碱
19、基 C 与碱基 A 的丰度之比、碱基 G 与碱基 C 的丰度之比则比较接近。3.将编号 1-40 人工制造的 DNA 序列的中，碱基 T 的丰度、碱基 G 的丰度、碱基T 与碱基 A 的丰度之比、碱基 G 与碱基 A 的丰度之比、碱基 C 与碱基 T 的丰度之比、碱基 G 与碱基 T 的丰度之比，用表格的形式加以表达（见附件 5，表 1） 。4.统计所有 DNA 序列中碱基 A、T、C、G 的比例，发现在未知类别的人工制造的DNA 序列以及自然序列中并非只存在 A、T、C、G 四种碱基，还存在n、s、w、y 等特殊碱基，这可能和生物自身需要完成的特定功能有关，具体列表如下：表 2 特殊的 DN
20、A 序列及特殊碱基种类DNA 序列 特殊碱基 DNA 序列 特殊碱基人工37 号 s 自然131 n自然71 n 自然147 n自然101 n、s 自然169 n自然105 r、s、w、y由上表可知，编号 1-20 的人工制造的 DNA 序列中并未出现特殊碱基，所以在提取特征时不需要考虑特殊碱基的影响，同样，在处理编号 21-40 的人工制造的 DNA 序列以及 182 条自然序列时，也不必考虑特殊碱基的影响，使用数据时，可将特殊碱基直接剔除。七、模型的建立与求解7.1 模型一：聚类分析延拓模型要使 DNA 序列的分类能够尽量科学合理，集中要解决的问题是让分类后的样品满足：同类样品间的差异性尽
21、可能小，不同类样品间的差异性尽可能大。为达到上述目的，引入聚类分析模型对不同的 DNA 序列进行分类。7.1.1 模型一的建立聚类分析方法根据分类对象的不同可以分为两类：1.对样品所作的分类，即 Q-型聚类，2.对变量所作的分类，即 R-型聚类。此问题将给出的不同 DNA 序列看成是不同的样品，选用 Q-型聚类进行具体求解。（1）样品特征变量的引入为了刻画不同样品的性质，需要对样品引入统一的特征作为样品特征变量，特征变量的确定来源于聚类分析前对数据进行预处理得到的分析结果。1）样品中 A，C ，T ，G 的碱基丰度样品 i 中 A 碱基丰度的计算：/AiFnN（1）其他碱基丰度的计算方法同上。
22、绘出编号 1-20 的人工制造的已知类别的 DNA 序列中 4 种碱基丰度的离散统计图（图 6.1.1） 。观察该散点图，进行数据分析可得：DNA 序列中碱基 A 和碱基 C 在分类 A和 B 中的区分不大，均大致在相同的频率区间内波动，故不选用碱基 A 和碱基C 的丰度作为特征区分；而 DNA 序列中碱基 T 和碱基 G 在分类 A 和 B 中的区分较大，A 类和 B 类相应的碱基丰度分别集中在不同的频率区间范围内，故选用碱基 T 和碱基 G 的丰度作为特征区分。将 T 的碱基丰度作为样品的第 1 个特征变量，记为 。1x将 G 的碱基丰度作为样品的第 2 个特征变量，记为 。22）样品不同
23、碱基间的比例样品 i 中碱基 T 和碱基 A 的比值计算：/GTTfn（2）其他碱基比例的计算方法同上。绘出编号 1-20 的人工制造的已知类别的 DNA 序列中不同碱基的丰度之比的离散统计图（图 6.1.2） 。观察该散点图，进行数据分析可得：DNA 序列中碱基 T 和碱基 A 的丰度之比以及碱基 G 和碱基 T 的丰度之比在分类 A 和 B 中的区分较大，A 类和 B 类相应的碱基丰度之比分别集中在不同的频率区间范围内，故选用碱基 T 和碱基 A的丰度之比以及碱基 G 和碱基 T 的丰度之比作为特征区分。将碱基 T 和碱基 A 的比值作为样品的第 3 个特征变量，记为 。3x将碱基 G 和
24、碱基 T 的比值作为样品的第 4 个特征变量，记为 。4（2）样品特征数据的向量转化把上述得到的 4 种特征变量分别作为一个向量的四个分量，用该向量作为样品特征向量来描述不同样品。由附件 5 表 1，编号 1-40 样品的 、 、 和 的值分别为表中的第1x234x1、2 、 3、6 列。于是得到编号 1-20 的样品的 20 个特征向量如下：； ；1(0., .964, 0.5, .93)Y2(0.5, .1, 0.567, 2.9)Y； ；3271448380； ；5(.8, .3, ., .6)(.6, .94, ., 3.14)； ；7019640589Y8012067952Y； ；9
25、(.2, ., .71, 2.3)(.3, ., ., .0)；15125418；13(0.82, .173, .05, .46)Y14(0., .2, .67, 0.34)Y；15(.64, ., .8, .9)16(.3, .9, 1.5, .91)；17(0.23, .1, 0.7436, .52)Y18(0.5, ., .78, 0.)Y；19(.56, ., 2.8, .9)20(.63, ., 2.1, .9)。（3）不同样品的相似性度量（分析编号 1-20 的样品）因为 20 个已知类别的 DNA 序列的样品变量均属于数值型数据，所以两个个体的相似度是指它们在欧氏空间中互相邻近的
26、程度。据此，引用距离测度来描述不同样品的相似性。距离测度小的两个样品，相似性较高；反之，距离测度大的两个样品，相似性较低。为了排除不同变量之间的相互影响，以及减弱较大数据出现时对结果的不良影响，即减弱较大值（包括异常值）的敏感度。选用 Lance 和 Williams 距离来描述距离测度，进而衡量不同样品间的相似性。此外，Lance 和 Williams 距离还与样品变量的单位无关，使结果无量纲化。向量 和向量 间的 Lance 和 Williams 距离为：iYj1|()pikjijxd（3）用公式（3）计算所有向量所代表的 维空间中所有样品点之间的两两距离。由排列组合知识，所有向量（n 个
27、）进行两两组合的个数为： ，分别计2nC算出每个组合的 Lance 和 Williams 距离。本次聚类中选用的向量个数为 n=20，一共有 种组合，用 matlab2019C编程（见附件 6 ）求解出所有组合的 Lance 和 Williams 距离，并对数据进行比较得出 。max3.71d（4）根据距离测度进行分类1）样品数据分成两类由上述得到的 ，查找 所对应的向量组合，假定该向量组合是向量maxmaxd和向量 ，则将第 i 个样品和第 j 个样品分为 A，B 两类，可以令 i 样品为 AiYj类，令 j 样品为 B 类。分别将 和 作为 A，B 两类的标准向量 ， 对剩余iYj AyB
28、样品进行分类。2）剩余样品分类样品 i 和样品 j 分完类后，还剩余（n-2）个样品未进行分类，将这（n-2）个样品数据分别和 A 类的标准向量 进行组合，计算出每个组合的 Lance 和AyWilliams 距离，将所得的距离进行比较，得出最小的 ，查找 所对应的向mindmin量，假定该向量是 ，则将该向量和样品 i 分为一类，同属于 A 类。用同样的aY方法把这（n-2 ）个样品数据分别和 B 类的标准向量 进行组合，得出最小的By，假定该组合所对应的向量是 ，则将该向量和样品 j 分为一类，同属于mindbYB 类。此时得到 A 组为 ， 。B 组为 ， 。aibjA，B 两类标准的重
29、新计算：将此时 A，B 组中的所有向量分别求出平均值得到 A，B 类的新的标准向量。A 类的标准向量：()/2AaiyY（4）B 类的标准向量：()/Bbjy（5）3）上述步骤后还剩余（n-4）个样品未进行分类，依照 2）剩余样品分类给出的方法不断重复进行计算，对所有的剩余样品均实现分类。7.1.2 模型一的求解按照上述方法首先计算得到这些样品中向量 和向量 间的 Lance 和3Y20Williams 距离最大，则将第 3 个样品和第 20 个样品分为 A，B 两类。令第 3 个样品为 A 类，第 20 个样品为 B 类。按照 7.1.1 中的步骤依次进行分类，用matlab 编程（见附件
30、7）求解得到分类结果如下：A 类的样品编号为：1，2，3，5，6，7，8，9，10，17；B 类的样品编号为：4，11，12，13，14，15，16，18，19，20。7.1.3 模型一的检验与改进（1）模型一的改进与可行性分析由以上分类结果可知，用聚类分析延拓法对编号 1-20 人工制造的 DNA 序列进行分类的结果与已知分类结果并非完全一致。在此分类方法下，第 4 条DNA 序列不再属于 A 类，而属于 B 类；第 17 条 DNA 序列不再属于 B 类，而属于 A 类。因此，有必要对模型进行改进。可以改变样品变量的组合方式，选择碱基 T 的丰度、碱基 T 与碱基 A 的丰度之比、碱基 C
31、 与碱基 T 的丰度之比、碱基 G 与碱基 T 的丰度之比作为四个样品变量，分别设为 、 、 和 。1x234x由附件 5 表 1，编号 1-40 样品的 、 、 和 的值分别为表中的第123x41、3 、 5、6 列。得到编号 1-20 的样品的 20 个特征向量如下：； ；1(0., .4, .267, .93)Y2(0.5, .67, 1.058, 2.9)Y； ；36381448396.； ；17(0.263, .74, 0.931,.572)Y18(0.5,.718, 0.2364, .18)Y；958 63 59。用公式（3）计算 20 个向量所代表的 4 维空间中所有样品点两两之
32、间的Lance 和 Williams 距离，并按照 7.1.1 中的距离测度法对编号 1-20 人工制造的DNA 序列进行分类得到的分类结果如下：A 类的样品编号为：1，2，3，4，5，6，7，8，9，10；B 类的样品编号为：11， 12，13，14，15，16，17，18，19，20。由以上分类结果可知，改变样品变量的组合方式，选择碱基 T 的丰度、碱基 T 与碱基 A 的丰度之比、碱基 C 与碱基 T 的丰度之比、碱基 G 与碱基 T 的丰度之比作为四个样品变量后，用聚类分析延拓法对编号 1-20 人工制造的 DNA序列进行分类的结果与已知分类结果完全一致。所以，该分类方法可行。（2）模
33、型一的进一步检验与实践1）用模型一中改进后的聚类分析延拓法，对编号 21-40 人工制造的 DNA序列进行分类，对附件 7 中的程序稍作修改，求解得到分类结果如下：A 类的样品编号为：22，23，25，27，29，30，34，35，36，37，39；B 类的样品编号为：21， 24，26，28，31，32，33，38，40。2）用模型一中改进后的聚类分析延拓法，对 182 个自然 DNA 序列进行分类，同样对附件 7 中的程序稍作修改，求解得到分类结果如下：B 类的样品编号为：7，10， 12， 22，23，24， 26，28，30，34，43， 48，50，54，57，65，75，76，80
34、， 84， 85，86，92， 98，103，107，110，114 ，116，119，121，122，123，127，128，129，130 ，131，137，138，140， 142，143，144，146，151，156， 159， 161，162， 163，166，168，170，173 ，174，175，179，180，181，182；其余的自然 DNA 序列为 A 类。7.1.4 模型一改进后分类效果的评价（1）求出 A 类中 10 条 DNA 序列 4 个样品变量（碱基 T 的丰度、碱基 T 与碱基A 的丰度之比、碱基 C 与碱基 T 的丰度之比、碱基 G 与碱基 T 的丰度之比
35、）的平均值，作为 A 类的标准点 a；求出 B 类中 10 条 DNA 序列 4 个样品变量的平均值，作为 B 类的标准点 b：a=（0.1393，0.5311，1.5172 ，3.2803） ；b=（0.5020 ，1.8300，0.2618 ，0.2131） 。（2）计算 A 类中 10 个样品点与标准点 a 之间的 Lance 和 Williams 距离，并求出距离的平均值和标准差：平均值 ，标准差 ；计算 B 类10.72410.5298a中 10 个样品点与标准点 a 之间的 Lance 和 Williams 距离，并求出距离的平均值和标准差：平均值 ，标准差 。2.6932.6a（
36、3）计算 A 类中 10 个样品点与标准点 b 之间的 Lance 和 Williams 距离，并求出距离的平均值和标准差：平均值 ，标准差 ；计算 B 类1.510.593b中 10 个样品点与标准点 b 之间的 Lance 和 Williams 距离，并求出距离的平均值和标准差：平均值 ，标准差 。20.720.3b（4）对以上数据进行分析。若分类方法合理，那么不同类别之间的差别应尽可能大，即 与 的差别、1a2与 的差别应尽可能大；同类之间的差别应尽可能小，即 、 、 和1b2 1b应尽可能小。此外，定义相关系数：，12|XrX 为 a 时，表示选择标准点 a 进行评价时的相关系数，；X
37、 为 b 时，表示选择标准点 b 进行评价时的相关系120.548|r数， 。12.79|b由均值和标准差的含义，为使 A 类与 B 类之间的差别尽可能大，那么相关系数 r 应该尽可能小，由以上结果 和 的大小均为 0.5 左右，可知该分类方arb法合理，且能够达到较好的分类效果。7.2 模型二：马尔可夫模型7.2.1 模型的建立与求解（1）DNA 序列的马尔可夫链转换把 DNA 的一个样品序列看成是一个系统，组成该 DNA 序列的不同位置的碱基看成是这个系统中的相应的不同状态。DNA 的长度为 N，则该系统有 N 个状态，分别记为 ,每个状态对应一个碱基。这样给定的一条长度为 N1,2.NS
38、的 DNA 序列转化成有 N 个状态组成的系统，即为 。随着时间的推移，1,2.NS系统从某一状态转移到另一状态，设 为时间 t 的状态（ ） ，系统在 t 时tqtjq间的状态 只与其在时间 t-1 的状态相关（ ） ，其概率为 。这1tq 1tiS 1(|)tPq样将该系统转换成一个离散的一阶马尔可夫链。（2）不同碱基的组合情况将 4 个碱基进行两两组合，用表格的形式进行考虑。表 3 两个碱基的组合情况碱基 A 碱基 T 碱基 C 碱基 G碱基 A AA AT AC AG碱基 T TA TT TC TG碱基 C CA CT CC CG碱基 G GA GT GC GG两个碱基组合排列一共有
39、16 种情况。 （ 表示在前一状态为碱基 A 的PAT情况下后一状态出现碱基 T 的概率，其他字母表示意义和上述同）（3）中间状态发生的概率中间状态：系统中除第一个状态和最后一个状态外均称为中间状态。以 A 类情况下 AT 的情况为例进行计算：给定的 A 类样品为编号 1-10 的 DNA 序列，将这 10 条 DNA 序列组合作为一个大样品。表示在该样品中出现碱基 A 的个数， 表示一条 DNA 链中碱基 A 后出n ATn现碱基 T 的情况组合起来看成一个新的碱基（AT），该样品中碱基（AT）的个数。的概率： （6）PAPT(|)/AT根据上述的计算公式，算出 A 类其余 15 种组合的概
40、率。B 类的计算情况和 A 类的情况相同。根据 A 类，B 类给出的 DNA 序列具体进行计算，得到以下表格：表 4 A 类不同中间状态发生的概率表碱基 A 碱基 T 碱基 C 碱基 G碱基 A P(AA)=0.3662 P(AT)= 0.1911 P(AC)= 0.1911 P(AG)= 0.2516碱基 T P(TA)= 0.2456 P(TT)= 0.3275 P(TC)= 0.1696 P(TG)= 0.2573碱基 C P(CA)= 0.2359 P(CT)= 0.1385 P(CC)= 0.0974 P(CG)= 0.5282碱基 G P(GA)=0.2601 P(GT)= 0.0
41、644 P(GC)= 0.2029 P(GG)= 0.4726表 5 B 类不同中间状态发生的概率表碱基 A 碱基 T 碱基 C 碱基 G碱基 A P(AA)=0.3707 P(AT)= 0.4081 P(AC)= 0.1184 P(AG)= 0.1028碱基 T P(TA)= 0.2701 P(TT)= 0.5931 P(TC)= 0.0858 P(TG)= 0.0511碱基 C P(CA)= 0.2182 P(CT)= 0.4727 P(CC)= 0.1636 P(CG)= 0.1455碱基 G P(GA)=0.3063 P(GT)= 0.3964 P(GC)= 0.0811 P(GG)=
42、 0.2162（4）待判定分类 DNA 序列的概率计算给定一条长度为 N 的 DNA，将其转换为系统状态序列 ，每一个1,2.NS系统状态对应同一位置 DNA 序列给出的一个碱基，计算该 DNA 序列产生的概率。该 DNA 序列系统产生的概率计算公式：123121321Pn n（ S,.,） =P(S)|(|).P(|)（7 ）第一个状态的出现概率均设为 1,即 。1()=分别根据 A 类，B 类给出的中间状态出现的概率，得到该 DNA 序列产生概率。（5）DNA 分类的判定将上面得到的两个 DNA 序列产生概率经行比较，如果通过 A 类中间状态的概率计算值远远大于 B 类中间状态的概率计算值
43、 ,则将该状态归为 A 类；同样，若通过 B 类中间状态的概率计算值远远大于 A 类中间状态的概率计算值,则将该状态归为 B 类。（6）实际数据的代入计算对已知类别的 20 个样品依照上述方法进行分类。表 6 编号为 1-10 个样本产生概率统计表样品1样品2样品3样品4样品5样品6样品7样品8样品9样品10A 类数据计算0.10e-60 0.14e-590.14e-560.21e-650.14e-590.35e-560.63e-590.31e-610.74e-610.79e-62B 类数据计算0.14e-760.13e-750.86e-800.89e-650.60e-800.12e-730.
44、22e-710.24e-740.74e-820.52e-82所属类别 A A A B A A A A A A表 7 编号为 11-20 个样本产生概率统计表样品11样品12样品13样品14样品15样品16样品17样品18样品19样品20A 类数据计算0.25e-63 0.44e-640.28e-650.64e-670.28e-600.27e-650.51e-710.28e-650.20e-680.71e-70B 类数据计算0.20e-490.89e-510.97e-540.97e-550.26e-410.61e-540.29e-680.11e-520.19e-530.12e-54所属类别B B
45、 B B B B B B B B根据上面的判断结果，只有 4 号样品的类别出现了偏差，说明马尔可夫模型进行判定具有一定的合理性。可以进一步推广，对其他的 DNA 序列进行判别7.2.2 分类结果统计用 matlab 编写程序（见附件 8）对编号 21-40 人工 DNA 序列以及 182 个自然序列进行分类，结果如下。（1）编号 21-40 人工 DNA 序列分类结果：A 类:22，23，25，26，27，29，30，32，33，34，35，36,37,39；B 类: 21,24,28,31,38,40。（2）182 个自然序列分类结果：A 类：1，2，3，4，5，6，8，9，10，11，12
46、，13，14，15，16，17，18，19，20，21，22，24，25，26，27，28，29，31，32，33，34，35，36，37，38，39，40，41，42，44，45，46，47，48，49，50，51，52，53，54，55，56，57，58，59，60，61，62，63，64，66，67，68，69，70，71，72，73，74，77，78，79，80，81，82，83，84，86，87，88，89，90，91，93，94，95，96，97，98，99，100，101，102，104，105，106，108，109，111，112，113，115，116，117，118，12
47、0，121，123，124，125，126，127，129，130，132，133，134，135，136，137，139，140，141，142，143，145，146，147，148，149，150，152，153，154，155，157，158，160，164，165，167，168，169，171，172，173，174，175，176，177，178，179，180，181B 类：7， 23，30，43，65 ，75，76，85，92，103，107 ，110， 114，119，122，128，131，138，144，151 ，156，159，161，162 ，163，166，170
48、，182总计，182 个自然序列分入 A 类的个数：154 个；182 个自然序列分入 B 类的个数：28 个7.2.3 模型的评价对模型的分类效果进行评价，评价标准：同类样品间的差异性较小，不同类样品间的差异性较大，则这样的分类效果较好；反之同类样品间的差异性较大，而不同类样品间的差异性较小，这样的分类效果就不够理想。（1）检验样本的数据处理选取前 20 号 DNA 样品作为检验样本。依照上述马尔可夫模型的分类结果，A 类为：1，2，3,5,6,7,8,9,10。B 类为：4,11,12,13,14,15,16,17,18,19,20。将每一个样品的最大产生概率（每个样品有两个生成概率，分别
49、用 A 类和B 类的相应数据进行求解，其中较大的一个即为该样品的最大产生概率）作为该样品的特征变量。由于最大产生概率的数量级很小，为了后期数据处理的方便和准确性，将各样品最大产生概率通过分别取对数的方法进行数据处理。A 类样品对应的计算结果为：-61.0000，-59.8539，-56.8539，-59.8539，-56.4559，-59.2007，-61.5086，-61.1308，-62.1024B 类样品对应的计算结果为：-65.0506，-49.6990，-51.0506，-54.0132，-55.0132，-41.5850，-54.2147，-68.5376，-52.9586，-53.7212，-54.9208。（2）同类样品间的差异性判定对上面的数据，按照类为单位分别计算相应的平均值，得到 A 类的平均值为： ，B 类的平均值为： 。_60.38A_53.714B用公式 分别求 A 类和 B 类的无偏方差。 （ 为方差，2_21()niiS2S为第 i 个最大概率取对数的值， 为平均概率取对数的值， 为标准差）i_ 这样进一步得到 A 类和 B 类的标
