数学建模_DNA序列分类模型论文终稿.doc

dna序列分类模型 dna序列分类模型毕业设计（论文）原创性声明和使用授权说明原创性声明本人郑重承诺：所呈交的毕业设计（论文），是我个人在指导教师的指导下进行的研究工作及取得的成果。尽我所知，除文中特别加以标注和致谢的地方外，不包含其他人或组织已经发表或公布过的研究成果，也不包含我为获得 及其它教育机构的学位或学历而使用过的材料。对本研究提供过帮助和做出过贡献的个人或集体，均已在文中作了明确的说明并表示了谢意。作 者 签 名： 日 期： 指导教师签名： 日期： 使用授权说明本人完全了解 大学关于收集、保存、使用毕业设计（论文）的规定，即：按照学校要求提交毕业设计（论文）的印刷本和电子版本；学校有权保存毕业设计（论文）的印刷本和电子版，并提供目录检索与阅览服务；学校可以采用影印、缩印、数字化或其它复制手段保存论文；在不以赢利为目的前提下，学校可以公布论文的部分或全部内容。作者签名： 日 期： 学位论文原创性声明本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任何其他个人或集体已经发表或撰写的成果作品。对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本人完全意识到本声明的法律后果由本人承担。作者签名： 日期： 年 月 日学位论文版权使用授权书本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅。本人授权 大学可以将本学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。涉密论文按学校规定处理。作者签名：日期： 年 月 日导师签名： 日期： 年 月 日注 意 事 项1.设计（论文）的内容包括：1）封面（按教务处制定的标准封面格式制作）2）原创性声明3）中文摘要（300字左右）、关键词4）外文摘要、关键词 5）目次页（附件不统一编入）6）论文主体部分：引言（或绪论）、正文、结论7）参考文献8）致谢9）附录（对论文支持必要时）2.论文字数要求：理工类设计（论文）正文字数不少于1万字（不包括图纸、程序清单等），文科类论文正文字数不少于1.2万字。3.附件包括：任务书、开题报告、外文译文、译文原文（复印件）。4.文字、图表要求：1）文字通顺，语言流畅，书写字迹工整，打印字体及大小符合要求，无错别字，不准请他人代写2）工程设计类题目的图纸，要求部分用尺规绘制，部分用计算机绘制，所有图纸应符合国家技术标准规范。图表整洁，布局合理，文字注释必须使用工程字书写，不准用徒手画3）毕业论文须用a4单面打印，论文50页以上的双面打印4）图表应绘制于无格子的页面上5）软件工程类课题应有程序清单，并提供电子文档5.装订顺序1）设计（论文）2）附件：按照任务书、开题报告、外文译文、译文原文（复印件）次序装订3）其它摘要本文分析了已知类别的人工dna序列的特征，建立了聚类分析延拓模型和马尔可夫模型，分别对未知类别的人工dna序列和自然序列进行分类，根据分类效果选出了较优模型。首先对数据进行预处理，得到人工dna序列的单个碱基丰度和不同碱基丰度之比等特征量，进而分析a、b两类的差异，得到合适的特征判定条件对未知类别的dna序列进行分类。计算人工dna序列的特征量，给出各序列的统计数据。其次用聚类分析延拓模型进行分类。用a、b两类具有明显差异的特征作为样品特征变量，得到欧式空间中表征编号1-20人工dna序列的特征向量，计算两两之间的lance和williams距离进行相似性度量，逐步选择相似性较大的归为一类，同时不断更新类内的标准比较特征向量，对聚类方法进行延拓，最终得到类内差异小、类间差异大的a、b两类，建立了聚类分析延拓模型。再对选取的特征变量进行改进，提高模型的分类效果。最后，借助均值、方差和相关系数等参数对改进模型的分类效果进行分析。 再次用马尔可夫模型进行分类。将dna序列看成是马尔可夫链，求出编号1-10和11-20人工dna序列在已知当前碱基种类的条件下，下一个碱基出现任一种的概率，结果存入概率转移矩阵1和2，再利用矩阵1和2分别求出编号1-20中任一条dna序列出现的概率，选择较大的一个作为该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，12，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年6月，人类基因组计划中dna全序列草图完成，预计2001年可以完成精确的全序列图，此后人类将拥有一本记录着自身生老病死及遗传进化的全部信息的“天书”。 （2）这本 “天书”是由4个字符a，t，c，g按一定顺序排成的无间隔的长约30亿的序列，除了这4个字符表示4种碱基以外，人们对它包含的“内容”知之甚少。因此，破译这部世界上最巨量信息的“天书”是二十一世纪最重要的任务之一。（3）为解读这部“天书”，首先要研究dna全序列具有什么结构，以及由这4个字符排成的看似随机的序列中隐藏着什么规律，这也是生物信息学最重要的课题。1.2题目信息（1）dna序列分为编码区与非编码区。编码区是用于编码蛋白质的序列片段，即由这4个字符组成的64种不同的3字符串，其中大多数用于编码构成蛋白质的20种氨基酸。（2）在不用于编码蛋白质的序列片段中，a和t的含量特别多些，于是以某些碱基特别丰富作为特征去研究dna序列的结构也取得了一些结果。（3）利用统计的方法还发现序列的某些片段之间具有相关性。这些发现说明dna序列中存在着局部的和全局性的结构，充分发掘序列的结构对理解dna全序列有十分重要的意义。目前在这项研究中最普通的思想是省略序列的某些细节，突出特征，然后将其表示成适当的数学对象。 1.3题目要求 （1）有20个已知类别的人工制造的dna序列（见附件1），其中序列标号110 为a类，11-20为b类。从中提取特征，构造分类方法，并用这些已知类别的序列，衡量所选分类方法是否足够好。（2）用（1）中的分类方法对另外20个未标明类别的人工序列（见附件1，标号2140）进行分类，根据分类效果对方法不断完善，将得到的最终结果用序号（按从小到大的顺序）标明它们的类别（a类或b类，无法分类的不写入）。 要求详细描述所选的分类方法，给出计算程序。若论文中部分地使用了现成的分类方法，应将方法名称准确注明。 （3）已知182个自然dna序列（见附件2），它们都较长。同样用以上所选的分类方法对它们进行分类，并根据分类效果对方法不断完善，像（2）中一样给出最终的分类结果。 二、 名词解释1.编码区与非编码区：编码区是指dna上编码蛋白质的序列片段，而非编码区不用于编码蛋白质。2.聚类分析：由已知数据，计算各个观察个体或变量之间亲疏关系的统计量。再根据某种准则（最短距离法、最长距离法、中间距离法、重心法等），使同一类内的差别较小，而类与类之间的差别较大，最终将观察个体或变量分为若干类的分类方法。其中，对样品所作的分类为q-型聚类，对变量所作的分类为r-型聚类。3.相似性度量：对数值型数据而言，两个个体的相似度是指它们在欧氏空间中互相邻近的程度；而对分类型数据而言，两个个体的相似度与它们取值相同的属性的个数有关。4.样品：每个观察个体即每条dna序列为一个样品。5.样品变量：每个样品所具有的不同特征用不同的变量来表示，变量数等于特征数。6.碱基丰度：每条dna序列中碱基a、g、c或t出现的频率。三、 问题分析dna序列分类问题要求在对dna序列的一些规律和结构有所了解的基础上，从20个已知类别的人工制造的dna序列中提取特征，构造分类方法，并用所选择的分类方法对其余未知类别的20个人工制造的dna序列以及182个自然dna序列进行分类。3.1建模目标的分析dna序列分类是一个复杂的统计分析问题，数据量大，影响因素多，无法直接从20条已知类别的人工制造的dna序列中提取出所有的有效特征，因此有必要对这20条dna序列进行预处理。观察并分析数据预处理结果，归纳总结出a类和b类的有效特征，将其表示成适当的数学对象，并选择适当的分类方法，建立普遍意义下数学模型，再用得到的模型对其余未知类别的20个人工制造的dna序列以及182个自然dna序列进行分类。由题意，建立的数学模型应该保证分类结果具有以下特点：（1）类别间差异尽量大；（2）类别内差异尽量小；（3）样品能够尽可能的落入a、b范围，且只能落入其中的一个。3.2建模及求解方向1.分析已知类别的dna序列1-20的结构，提取出相应的特征。主要的特征有：碱基的丰度、碱基或碱基序列的重复出现情况、碱基或碱基序列之间的相邻情况、不同碱基的丰度之比（如碱基a与碱基t的丰度之比）等。2. 根据提取出的特征，选用合适的分类方法。 对数据进行预处理后，尝试以下方法建立模型：（1）根据聚类分析法，建立模型一。由题意，dna序列分类属于对样品所做的分类，为q-型聚类。首先引入样品变量，例如可选择碱基t的丰度、碱基g的丰度、碱基t与碱基g的丰度之比、碱基a与碱基t的丰度之比等。由已知数据，计算出每条已知类别的人工制造的dna序列的各个样品变量值，存入向量中。根据相似性度量原理，计算20个样品两两之间的lance和williams距离，选择相距最远的两个样品（假设为样品3和样品16）分别作为a类和b类，再分别以样品3和样品16为标准点，通过分别计算样品3和样品16与其余18个样品之间的lance和williams距离，找出与其相距最近的一个样品（假设为样品1和样品18）归为一类。此时，新的标准点变为样品1与样品3的中点、样品16与样品18的中点。然后再以新的标准点为基准，分别找出与其相距最近的一个样品归为一类。逐步进行下去，直至20个样品被明显分成a、b两类。（2）根据马尔可夫法，建立模型二。以单个碱基为单位，分别统计编号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、c、t的概率，存入两个矩阵1和2中。对于任何一条给定的dna序列，可将其看作一个已经发生的事件，说明该事件发生的概率比较大。用矩阵1和矩阵2分别求出这一事件发生（即形成当前dna序列）的概率，若用矩阵1算出该编号的dna序列出现的概率较大，则该编号的dna序列属于a类，否则属于b类。3.模型的初步检验与改进。用编号1-20已知类别的序列，分别衡量模型一与模型二中所选方法是否足够好，不断改进，尽可能使1-20号dna序列在所选分类方法下，所得结果与已知分类完全一致。改进时，对于聚类分析法，可以尝试改变样品变量的个数或者改变样品变量的组合方式；对于马尔可夫法，可以尝试引进中间变量，运用隐马尔可夫法求解。4.模型的进一步检验与完善。（1）用以上的得到的两种分类方法对编号20-40未知类别的人工序列、182个自然序列进行分类。（2）通过计算样品方差、均值等比较两种分类方法得到的分类结果与建模目标类别间差异尽量大、类别内差异尽量小、样品能够尽可能的落入a、b范围，且只能落入其中的一个的接近程度。（3）选择更接近建模目标的一种分类方法作为最终的分类方法，其分类结果即为最终结果。四、基本假设1.假设所给的dna序列片段中没有断句和标点符号。2.假设具有特殊碱基的dna序列中，特殊碱基可以剔除，其影响可以忽略。3.较长的182个自然序列与已知类别的20个样本序列具有共同的特征。4.假设给定的dna序列均是从全序列中随机截取出来的,无法确定序列的起始位, 无法从序列中辨认出氨基酸，所以，在对dna 序列分类时,从碱基层次上进行分类, 而不是从氨基酸层次上分类。五、定义与符号说明：各个样品中碱基出现的数量，i为a、t、c或g：第i个样品的总碱基数目：各个样品中碱基的丰度，i为a、t、c或g：各个样品的第i个特征变量：各个样品中碱基i和碱基j的比值，i，j为a、t、c或g：第i个样品的特征向量：向量和向量间的lance和williams距离：特征向量的分量个数，即向量的维数：特征向量的第k个分量：样品的个数：特征向量i的第k个分量：不同向量代表的维空间中任意两点间lance和williams距离的最大值：不同向量代表的维空间中任意两点间lance和williams距离的最小值：聚类分析中i类的标准向量，i为a或b六、数据预处理1.a类和b类样品单个碱基丰度的计算用matlab编写程序（见附件3），分别求出20条已知类别的人工制造的dna序列中，4种碱基的丰度，绘出散点图如下：图6.1.1 单个碱基丰度比较图分析上图可得， a类和b类dna序列中碱基t和碱基g的丰度有明显差异，而碱基a和碱基c的丰度则比较接近。2. a类和b类样品不同碱基丰度之比的计算用matlab编写程序（见附件4），分别求出20条已知类别的人工制造的dna序列中，不同碱基的丰度之比，包括、，绘出散点图如下：图6.1.2 不同碱基丰度之比的比较图分析上图可得， a类和b类dna序列中，碱基t与碱基a的丰度之比、碱基g与碱基a的丰度之比、碱基c与碱基t的丰度之比、碱基g与碱基t的丰度之比有明显差异，而碱基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 特殊的dna序列及特殊碱基种类dna序列特殊碱基dna序列特殊碱基人工37号s自然131n自然71n自然147n自然101n、s自然169n自然105r、s、w、y 由上表可知，编号1-20的人工制造的dna序列中并未出现特殊碱基，所以在提取特征时不需要考虑特殊碱基的影响，同样，在处理编号21-40的人工制造的dna序列以及182条自然序列时，也不必考虑特殊碱基的影响，使用数据时，可将特殊碱基直接剔除。七、模型的建立与求解7.1模型一：聚类分析延拓模型要使dna序列的分类能够尽量科学合理，集中要解决的问题是让分类后的样品满足：同类样品间的差异性尽可能小，不同类样品间的差异性尽可能大。为达到上述目的，引入聚类分析模型对不同的dna序列进行分类。7.1.1模型一的建立聚类分析方法根据分类对象的不同可以分为两类：1.对样品所作的分类，即q-型聚类，2.对变量所作的分类，即r-型聚类。此问题将给出的不同dna序列看成是不同的样品，选用q-型聚类进行具体求解。（1）样品特征变量的引入为了刻画不同样品的性质，需要对样品引入统一的特征作为样品特征变量，特征变量的确定来源于聚类分析前对数据进行预处理得到的分析结果。1）样品中a，c，t，g的碱基丰度样品i中a碱基丰度的计算： （1）其他碱基丰度的计算方法同上。绘出编号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个特征变量，记为。将g的碱基丰度作为样品的第2个特征变量，记为。2）样品不同碱基间的比例样品i中碱基t和碱基a的比值计算： （2）其他碱基比例的计算方法同上。绘出编号1-20的人工制造的已知类别的dna序列中不同碱基的丰度之比的离散统计图（图6.1.2）。观察该散点图，进行数据分析可得：dna序列中碱基t和碱基a的丰度之比以及碱基g和碱基t的丰度之比在分类a和b中的区分较大，a类和b类相应的碱基丰度之比分别集中在不同的频率区间范围内，故选用碱基t和碱基a的丰度之比以及碱基g和碱基t的丰度之比作为特征区分。将碱基t和碱基a的比值作为样品的第3个特征变量，记为。将碱基g和碱基t的比值作为样品的第4个特征变量，记为。（2）样品特征数据的向量转化把上述得到的4种特征变量分别作为一个向量的四个分量，用该向量作为样品特征向量来描述不同样品。由附件5表1，编号1-40样品的、和的值分别为表中的第1、2、3、6列。于是得到编号1-20的样品的20个特征向量如下：；。（3）不同样品的相似性度量（分析编号1-20的样品）因为20个已知类别的dna序列的样品变量均属于数值型数据，所以两个个体的相似度是指它们在欧氏空间中互相邻近的程度。据此，引用距离测度来描述不同样品的相似性。距离测度小的两个样品，相似性较高；反之，距离测度大的两个样品，相似性较低。为了排除不同变量之间的相互影响，以及减弱较大数据出现时对结果的不良影响，即减弱较大值（包括异常值）的敏感度。选用lance和williams距离来描述距离测度，进而衡量不同样品间的相似性。此外，lance和williams距离还与样品变量的单位无关，使结果无量纲化。向量和向量间的lance和williams距离为： （3）用公式（3）计算所有向量所代表的维空间中所有样品点之间的两两距离。由排列组合知识，所有向量（n个）进行两两组合的个数为：，分别计算出每个组合的lance和williams距离。本次聚类中选用的向量个数为n=20，一共有种组合，用matlab编程（见附件6 ）求解出所有组合的lance和williams距离，并对数据进行比较得出。（4）根据距离测度进行分类1）样品数据分成两类由上述得到的，查找所对应的向量组合，假定该向量组合是向量和向量，则将第i个样品和第j个样品分为a，b两类，可以令i样品为a类，令j样品为b类。分别将和作为a，b两类的标准向量，对剩余样品进行分类。2）剩余样品分类样品i和样品j分完类后，还剩余（n-2）个样品未进行分类，将这（n-2）个样品数据分别和a类的标准向量进行组合，计算出每个组合的lance和williams距离，将所得的距离进行比较，得出最小的，查找所对应的向量，假定该向量是，则将该向量和样品i分为一类，同属于a类。用同样的方法把这（n-2）个样品数据分别和b类的标准向量进行组合，得出最小的，假定该组合所对应的向量是，则将该向量和样品j分为一类，同属于b类。此时得到a组为，。b组为，。a，b两类标准的重新计算：将此时a，b组中的所有向量分别求出平均值得到a，b类的新的标准向量。a类的标准向量： （4）b类的标准向量： （5）3）上述步骤后还剩余（n-4）个样品未进行分类，依照2）剩余样品分类给出的方法不断重复进行计算，对所有的剩余样品均实现分类。7.1.2模型一的求解按照上述方法首先计算得到这些样品中向量和向量间的lance和williams距离最大，则将第3个样品和第20个样品分为a，b两类。令第3个样品为a类，第20个样品为b类。按照7.1.1中的步骤依次进行分类，用matlab编程（见附件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与碱基t的丰度之比、碱基g与碱基t的丰度之比作为四个样品变量，分别设为、和。由附件5表1，编号1-40样品的、和的值分别为表中的第1、3、5、6列。得到编号1-20的样品的20个特征向量如下：；.；。用公式（3）计算20个向量所代表的4维空间中所有样品点两两之间的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）模型一的进一步检验与实践 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，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的丰度之比）的平均值，作为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个样品点与标准点a之间的lance和williams距离，并求出距离的平均值和标准差：平均值，标准差。（3）计算a类中10个样品点与标准点b之间的lance和williams距离，并求出距离的平均值和标准差：平均值，标准差；计算b类中10个样品点与标准点b之间的lance和williams距离，并求出距离的平均值和标准差：平均值，标准差。（4）对以上数据进行分析。若分类方法合理，那么不同类别之间的差别应尽可能大，即与的差别、与的差别应尽可能大；同类之间的差别应尽可能小，即、和应尽可能小。此外，定义相关系数：， x为a时，表示选择标准点a进行评价时的相关系数，；x为b时，表示选择标准点b进行评价时的相关系数，。由均值和标准差的含义，为使a类与b类之间的差别尽可能大，那么相关系数r应该尽可能小，由以上结果 和的大小均为0.5左右，可知该分类方法合理，且能够达到较好的分类效果。7.2模型二：马尔可夫模型7.2.1模型的建立与求解（1）dna序列的马尔可夫链转换把dna的一个样品序列看成是一个系统，组成该dna序列的不同位置的碱基看成是这个系统中的相应的不同状态。dna的长度为n，则该系统有n个状态，分别记为,每个状态对应一个碱基。这样给定的一条长度为n的dna序列转化成有n个状态组成的系统，即为。随着时间的推移，系统从某一状态转移到另一状态，设为时间t的状态（），系统在t时间的状态只与其在时间t-1的状态相关（），其概率为。这样将该系统转换成一个离散的一阶马尔可夫链。（2）不同碱基的组合情况将4个碱基进行两两组合，用表格的形式进行考虑。表3 两个碱基的组合情况碱基a碱基t碱基c碱基g碱基aaaatacag碱基ttatttctg碱基ccactcccg碱基ggagtgcgg两个碱基组合排列一共有16种情况。（表示在前一状态为碱基a的情况下后一状态出现碱基t的概率，其他字母表示意义和上述同）（3）中间状态发生的概率中间状态：系统中除第一个状态和最后一个状态外均称为中间状态。以a类情况下at的情况为例进行计算：给定的a类样品为编号1-10的dna序列，将这10条dna序列组合作为一个大样品。表示在该样品中出现碱基a的个数，表示一条dna链中碱基a后出现碱基t的情况组合起来看成一个新的碱基（at），该样品中碱基（at）的个数。的概率： （6）根据上述的计算公式，算出a类其余15种组合的概率。b类的计算情况和a类的情况相同。根据a类，b类给出的dna序列具体进行计算，得到以下表格：表4 a类不同中间状态发生的概率表碱基a碱基t碱基c碱基g碱基ap(aa)=0.3662p(at)= 0.1911p(ac)= 0.1911p(ag)= 0.2516碱基tp(ta)= 0.2456p(tt)= 0.3275p(tc)= 0.1696p(tg)= 0.2573碱基cp(ca)= 0.2359p(ct)= 0.1385p(cc)= 0.0974p(cg)= 0.5282碱基gp(ga)=0.2601p(gt)= 0.0644p(gc)= 0.2029p(gg)= 0.4726表5 b类不同中间状态发生的概率表碱基a碱基t碱基c碱基g碱基ap(aa)=0.3707p(at)= 0.4081p(ac)= 0.1184p(ag)= 0.1028碱基tp(ta)= 0.2701p(tt)= 0.5931p(tc)= 0.0858p(tg)= 0.0511碱基cp(ca)= 0.2182p(ct)= 0.4727p(cc)= 0.1636p(cg)= 0.1455碱基gp(ga)=0.3063p(gt)= 0.3964p(gc)= 0.0811p(gg)= 0.2162（4）待判定分类dna序列的概率计算给定一条长度为n的dna，将其转换为系统状态序列，每一个系统状态对应同一位置dna序列给出的一个碱基，计算该dna序列产生的概率。该dna序列系统产生的概率计算公式： （7）第一个状态的出现概率均设为1,即。分别根据a类，b类给出的中间状态出现的概率，得到该dna序列产生概率。（5）dna分类的判定将上面得到的两个dna序列产生概率经行比较，如果通过a类中间状态的概率计算值远远大于b类中间状态的概率计算值,则将该状态归为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.22e-710.24e-740.74e-820.52e-82所属类别aaabaaaaaa表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所属类别bbbbbbbbbb根据上面的判断结果，只有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，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，120，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，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。将每一个样品的最大产生概率（每个样品有两个生成概率，分别用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类的平均值为：。用公式分别求a类和b类的无偏方差。（为方差，为第i个最大概率取对数的值，为平均概率取对数的值，为标准差）这样进一步得到a类和b类的标准差分别为，该标准差的值较小说明同类样品间的差异性较小。（3）不同类样品间的差异性判定将a类的9个样品分别用a类的中间状态概率和b类的中间状态概率进行计算得到两组平均值：，。（表示a类的9个样品用a类的中间状态概率得到平均值，表示a类的9个样品用b类的中间状态概率得到平均值）同理将b类的11个样品分别用a类的中间状态概率和b类的中间状态概率进行计算得到两组平均值：，。（表示b类的11个样品用a类的中间状态概率得到平均值，表示b类的11个样品用b类的中间状态概率得到平均值）a类样品的两个标准差:,(表示a类样品通过a类的中间状态概率计算得到的方差，表示a类样品通过a类的中间状态概率计算得到的标准差，表示a类样品通过b类的中间状态概率计算得到的方差，表示a类样品通过b类的中间状态概率计算得到的标准差)a类样品中分别用a类的中间状态概率和b类的中间状态概率进行计算，所得结果的差异性：（表示a类样品中两种中间状态概率计算结果的相关系数）把b类的11个样品用上述的方法做同样的处理，得到b类样品中两种中间状态概率计算结果的相关系数：上述相关系数间接刻画了两类别间的差异性大小，值越小，则说明不同类样品间的差异性越大；反之,值越大，则说明不同类样品间的差异性越小。选取两个系数中的较差值作为最后的相关性系数（依照最差值选取原则，若最差值也能取得较好的结果，则得到整个结果的合理性），该值较小，从而得到不同类样品间的差异性较大。依照上述给定的判断标准，得到马尔可夫模型的分类效果较好的结论。7.3模型一与模型二选择的层次分析模型本文给出了两种不同的模型对dna序列进行了分类，现用层次分析法综合对两者进行比较，选出其中更优的一个模型。该层次分析中，目标层为选取两个模型中更优的一个；准则层为模型好坏的评价标准，确定为两点：1.模型分类情况和给定前20个样品的分类情况的相似程度，2.模型自身的分类效果（同类样品间的差异性较小，不同类样品间的差异性较大，则自身的分类效果较好）；方案层为聚类分析延拓模型和马尔科夫模型。根据上述的层次，给出层次分析的结构图：图7.3.1 层次分析法的结构图1）构造准则层对于目标层的成对比较矩阵：，其中，a是表示模型与给定分类相似程度与模型自身分类效果重要性之比。求得权向量为：；参数确定：a取3，表示模型与给定分类相似程度比自身分类效果稍微重要。2） 构造方案层对准则层各因素的成对比较矩阵并求出权向量相对于指标层模型与给定分类相似程度： ，其中，b1表示聚类分析延拓模型与马尔科夫模型对于指标层模型与给定分类相似程度。 求得权向量为：；参数确定：由于聚类分析延拓模型与模型给定人工序列的分类结果完全吻合，而马尔科夫模型有一条链出现偏差，取b1=5;求得相对于指标层模型自身分类效果矩阵：，其中，b2表示聚类分析延拓模型与马尔科夫模型对于指标层模型自身分类效果；求得权向量为：，参数确定：由于聚类分析延拓模型和马尔科夫模型相比，模型自身分类效果较佳，取b2=3；由此求得方案最终的权重矩阵为： ，代入数据解得： .由结果可得聚类分析延拓模型明显优于马尔科夫模型，所以采用聚类分析延拓模型作为最终的dna序列分类模型。八、模型的评价与改进方向8.1模型的优点（1）对整个问题的解决，基于聚类分析法和马尔可夫法分别建立了模型一和模型二，从中选择较优的模型，避免单一模型解决问题，引起较大误差。（2） 建立模型一时，对传统的聚类分析法进行改进，建立了聚类分析延拓模型。（3）进行模型二的求解时，将结果存入矩阵，简化了运算，增加了效率。（4）在进行模型一和模型二的选择时，采用层次分析法，综合考虑评判分类效果的多个因素。8.2模型的缺点（1）在用马尔可夫模型进行分类时，编号1-20人工dna序列的分类结果没有达到与已知分类结果完全一致。（2）在进行模型一和模型二的最终选择时，综合考虑的因素还有欠缺。8.3模型的改进方向（1）对编号1-20人工dna序列特征的分析还有改进之处，例如可以考虑以三联子密码子为单位、以相邻四个碱基为单位，分析a类和b类的差异。（2）在模型二中，采用马尔可夫法对后20条人工dna序列和自然dna序列进行分类时，效果不是很理想，可以考虑引进中间变量，采用隐马尔可夫法，对模型进行改进。九、参考文献1徐全智、杨晋浩，数学建模（第二版），北京：高等教育出版社，2008.62姜启源、谢金星、叶俊，数学模型（第三版），北京：高等教育出版社，2003.83梁之舜、邓集贤、杨维权等，概率论与数理统计，中山大学数学系，1987.54彭放、杨瑞琰、罗文强等，数学建模方法，北京：科学出版社，2007.55刘祖洞，遗传学（第二版），北京：高等教育出版社，19916刘秉权，隐马尔可夫模型，2006.117吕金翅、马小龙、曹芳，dna 序列分类的数学模型，数学的实践与认识，第31卷第1期，2001.18樊龙江，序列分析与联配，生物信息学札记9李银山、杨春燕、张伟，dna 序列分类的神经网络方法，计算机仿真，第20卷第2期，2003.210孟大志，dna 序列中的结构与简化模型，数学的实践与认识，第31卷第1 期，2001 .1 11电子科技大学应用数学