题目DNA序列摘要本文主要研究DNA序列的结构问题,通过建立相应的数学模型,对DNA序列中所隐藏的规律进行研究和分析,给出了解决问题的最优方案,并且对模型进行了评价和推广。对于问题一,为了挖掘DNA序列的特征将其分为A类和B类,以20种基本氨基酸为目标,利用Matlab软件编程得出每一行每一种氨基酸出现的概率;再运用主成分分析法进行降维,利用SPSS软件进行数据处理得到矩阵;然后再将模糊聚类问题转化为如下优化问题:2111min(,)(())..1(1,2,6)01ncqikikkicikiikJUVudstuku用模糊聚类分析方法来获取样本与聚类中心的加权距离最小的最佳分类,使其分类结果与问题所给的答案一致,如果如下表:A类1235678910B类411121314151617181920由此可见,误差较小可忽略不计,证明此种方法可行。再以此种方法将21到40行DNA序列分为两类,其中A类,14681011131820;B类,235791214151617对于问题二:在问题一的基础上,对182个氨基酸序列进行分类,采取与问题一相同的方法进行分类,分类结果见问题二的求解。总的来说,本模型在未知数据特征的情况下很好的将数据进行分类,成功地解决了此次数学建模的DNA序列问题,是聚类分析问题的一个有效而且具有较强实用性的方法。关键词:主成分分析模糊聚类分析Matlab软件Spss软件1一、问题重述1.1背景分析随着DNA测序时代的到来,越来越多生物的全基因组序列正逐渐展现于人们的眼前。如何从中挖掘有用的信息成为对当今生物学乃至整个科学领域的一个挑战。本文主要致力于对DNA序列结构以及序列中所隐藏规律的研究。1.2问题重述2000年6月,人类基因组计划中DNA全序列草图完成,预计2001年可以完成精确的全序列图,此后人类将拥有一本记录着自身生老病死及遗传进化的全部信息的“天书”。这本大自然写成的“天书”是由4个字符A,T,C,G按一定顺序排成的长约30亿的序列,其中没有“断句”也没有标点符号,除了这4个字符表示4种碱基以外,人们对它包含的“内容”知之甚少,难以读懂。破译这部世界上最巨量信息的“天书”是二十一世纪最重要的任务之一。在这个目标中,研究DNA全序列具有什么结构,由这4个字符排成的看似随机的序列中隐藏着什么规律,又是解读这部天书的基础,是生物信息学(Bioinformatics)最重要的课题之一。虽然人类对这部“天书”知之甚少,但也发现了DNA序列中的一些规律性和结构。例如,在全序列中有一些是用于编码蛋白质的序列片段,即由这4个字符组成的64种不同的3字符串,其中大多数用于编码构成蛋白质的20种氨基酸。又例如,在不用于编码蛋白质的序列片段中,A和T的含量特别多些,于是以某些碱基特别丰富作为特征去研究DNA序列的结构也取得了一些结果。此外,利用统计的方法还发现序列的某些片段之间具有相关性,等等。这些发现让人们相信,DNA序列中存在着局部的和全局性的结构,充分发掘序列的结构对理解DNA全序列是十分有意义的。目前在这项研究中最普通的思想是省略序列的某些细节,突出特征,然后将其表示成适当的数学对象。这种被称为粗粒化和模型化的方法往往有助于研究规律性和结构。作为研究DNA序列的结构的尝试,提出以下对序列集合进行分类的问题:问题一:下面有20个已知类别的人工制造的序列(见附录),其中序列标号1—10为A类,11-20为B类。请从中提取特征,构造分类方法,并用这些已知类别的序列,衡量你的方法是否足够好。然后用你认为满意的方法,对另外20个未标明类别的人工序列(标号21—40)进行分类,把结果用序号(按从小到大的顺序)标明它们的类别(无法分类的不写入):A类;B类。请详细描述你的方法,给出计算程序。如果你部分地使用了现成的分类方法,也要将方法名称准确注明。2问题二:在同样网址的数据文件Nat-model-data中给出了182个自然DNA序列,它们都较长。用你的分类方法对它们进行分类,像1)一样地给出分类结果。二、模型假设结合本题实际,为确保模型求解的准确性和合理性,我们排除了一些因素的干扰,提出以下几点假设:1、DNA序列具有连续有序性;2、较长的182个自然序列与已知类别的20个样本序列具有共同的特征;3、64种3字符串压缩为20组后不影响分类的结果;4、各序列中DNA碱基三联组(即3字符串)的起始位置和基因表达不影响分类的结果。5、分类的结果在一定误差范围内认为合理;三、符号说明3.1为了便于问题的求解,我们给出以下符号说明:ix表示第i种氨基酸x表示20个变量构成的20维随机变量x原始数据的平均值a随机变量x的均值R协方差矩阵i矩阵的特征值iY主成分的贡献率ijl主成分与原始指标的相关系数mF氨基酸综合评分函数iku样本隶属于第i类的隶属度U模糊划分矩阵ikd样本kw与第i类聚类中心之间的距离iw分类标准原始数据标准差3.2名词解释:1、氨基酸:含有氨基和羧基的一类有机化合物的通称。生物功能大分子蛋白质的基本组成单位,是构成动物营养所需蛋白质的基本物质。2、模糊聚类分析法:聚类分析是对事物按一定要求进行分类的数学方法。实际的分类问题常伴有模糊性,因此,聚类问题用模糊数学的方法解决更确切。3、隶属度:若对论域(研究的范围)A中的任一元素y,都有一个数(0,1)iku与之对应,则称iku为A上的模糊集,iku称为y对A的隶属度。3四、问题分析4.1对问题一的分析问题一中需要20行DNA序列的特征并对其进行分类,使其分为两类,并使前十行属于A类后十行属于B类中,因为DNA序列较长,且为字符,不能直接对其进行分类,这里我们将其转化为三字符的氨基酸形式,以查询资料得到的20种基本氨基酸为目标,利用Matlab软件编程求出每一行这20种氨基酸分别出现的概率,以此作为实验数据。然而20种氨基酸种类仍然较多,需要对其进行主成分分析将其降维以方便下面的运算。数据处理完之后,建立相应的20*6的矩阵,利用模糊聚类分析方法从事先给出的C个划分出发,通过不断地反复修改样本的类别、聚类中心以及样本隶属于各类别的隶属度,来获得样本与聚类中心的加权距离最小的最佳分类。再将分类结果与已分好的类别作比较,判断误差是否较小结果是否正确。在以相同的模糊聚类分析方法将后20行分为A、B两类,分别填入题中横线。4.2对问题二的分析因为前面的模糊聚类分析方法经第一次分类的检验,证明其误差较小,结果比较准确,可以继续用这种方法对接下来的182个DNA序列进行分类。因为序列较长且多,所以可以用Matlab编程直接导入数据,算出每行20种氨基酸出现的概率,然而数据太过庞大,根据第一题可以看出DNA序列中各种氨基酸出现的概率之间方差较大,可以先对其进行主成分分析,然后重新建立目标矩阵直接编程对其进行分类。五、模型的建立与求解5.1问题一的模型建立与求解5.11根据对问题一第一小问的分析可知,附件里的DNA序列较长较多,都为字符,而且总体上只有四个字符数,不能直接对这些字符进行分类,并且这些字符所代表的DNA在生物上会三三结合形成氨基酸以创造出蛋白质,所以其有一定的连续性,必须将这些单个DNA字符转化成三个字符相结合而成的氨基酸形式来进行研究以找出其相应特征进行分类。根据资料可知,,,,ATCG这四个字符组成了64种不同的3字符串,这64种3字符串构成生物蛋白质的20种氨基酸;通过对这20种氨基酸出现概率的特征研究结合Matlab软件编程(程序见附录),得出前20行DNA序列的20*20概率矩阵,见附录。然而数据还是太大,氨基酸种类太多,不利于特征的研究;所以根据这20*20概率矩阵利用降维的思想将多个变量转化为少数几个综合变量对其进行主成分分析。主成分分析模型建立如下:假设原来研究对象是20种氨基酸,分别用1220,xxx来表示,这20个变量构成的20维随机向量为1220(,)txxxx。设随机向量x的均值为a,为反应标准化的数据之间相关关系密切程度的统计指标建立协方差矩阵为R,值越大,说明有必要对数据进行主成分分析,其中,(,1,2,,20)ijRij为原始变量ix,jy的相关系数。R为对称矩阵(即ijjiRR),只需计算其上三角元素或下三角元素即4可,其计算公式为:1221()()()()nkjikjjkijnkjikjjkxxxxRxxxx然后对x进行线性变化,考虑原始变量的线性组合:再根据协方差矩阵R求出特征值、主成分贡献率和累计方差贡献率,确定主成分个数。特征值是各主成分的方差,它的大小反映了各个主成分的影响力。解特征方程0,ER求出(1,2,,20)ii。因为R是正定矩阵,所以其特征值i都为正数,将其按大小顺序排列,即120i。特征值是各主成分的方差,它的大小反映了各个主成分的影响力。主成分iy的贡献率1/pijjjY,累计贡献率为11/pmjjjj。根据选取主成分个数的原则,特征值要求大于1且累计贡献率达80%-95%的特征值12,,m所对应的1,2,···,()mmp,其中整数m即为主成分的个数。接下来建立初始因子载荷矩阵,解释主成分。因子载荷矩阵是主成分iy与原始指标ix的相关系数(,)(,1,2,,20)ijiiiijlRyxeij,揭示了主成分与氨基酸比率之间的相关程度。过程详细数据见附件SPSS软件运行结果图。最后根据特征值及累计贡献率的数值,分析得到3个主成分,根据筛选出来的3个主成分重新设立成分矩阵(详情见附录图1)。利用模糊聚类分析将前二十行DNA序列按照各个样本点之间、样本点与样本点子集合之间的关系体系确定氨基酸之间的“亲疏关系”从而对其进行正确而合理的分类。模糊聚类分析模型建立过程如下:设待分类的m个样本为12,6,,个特征值,即12{,,,}iiiin,根据题意可知需要将样本划分为2类。为了获得一个最佳的模糊分类,需要一个分类准则,定义目标函数为样本与聚类中心的加权距离,加权系数为样本的隶属度函数的q次方。于是,模糊聚类问题可以转化为如下优化问题如下:2111min(,)(())..1(1,2,6)01ncqikikkicikiikJUVudstuku其中iku为样本隶属于第i类的隶属度,模糊划分矩阵26[]ikUu;12{,}|Vvv为1111122120202211222220202020112022202020yaxaxaxyaxaxaxyaxaxax5两个聚类中心集合;[1,]q为加权指数,当1q时,模糊聚类就退化为硬2均值聚类;通常q的最佳选择范围为[1.5,2.5],一般2q是比较理想的取值。ikd表示样本kw与第i类聚类中心之间的距离,定义为2()()TikkikikidwvwuAwu式中A为pp的正定矩阵,当1A时,即为欧氏距离。样本kw隶属于i类的程度,即隶属度,iku可在0到1之间取值。这样,样本kw不在明确的属于某一类,而是对于每一个类别都有一个隶属度,隶属度的数值越大表明样本隶属于该类别的程度越大,反之则说明隶属于该样本的程度越小。模糊聚类的这种模糊划分描述了样本聚类过程中的模糊现象,从而可以获得更为合理的聚类结果。为了使目标函数(,)JUV到达极小点,UV,其应该满足:122/1111,,()nqikkkiiknqqikikkkjkUXVUdUd利用上述公式,运用迭代方法通过Matlab软件(编程见附录)可以获得聚类分析结果,如下:A1235678910B411121314151617181920根据模糊聚类分析得出的结果,1到10行除了第4行均属于一类,11到20行以及第4行属于另一类。与题中给出的序列标号1—10为A类,11-20为B类的分类情况相比,只有第4行出现误差。与数据庞大的题目相比,误差相对较小,可以忽略不计;所以用模糊聚类分析结合