DNA序列问题模型

1、2015年芜湖三校数学建模竞赛题 目 DNA序列问题模型摘 要 DNA 序列是由A，T，G，C 四个表示4 种碱基的字符组成的序列。本文研究DNA 序列的结构找出序列间的差异和对八个物种的DNA序列进行分类。对于问题一首先对数据运用数理统计方法对数据进行计算，得到八个物种的DNA序列的碱基的丰度、碱基的重复出现情况、碱基之间的相邻情况、不同碱基的丰度之比的四个特征，通过对这四个特征作出相应的散点图比较得出八个物种的DNA序列间的差异：Human 、Opossum、 Lemur 、Rat等4种DNA序列的长度相同，其他四种DNA序列的长度各不相同，每种DNA序列四种碱基的的重复情况也各不相同；G

2、碱基的丰度相对于本序列的其他碱基的丰度都要高，碱基A的丰度在各个序列中丰度差不多，其他三种碱基在序列中波动性较大，差异性较大；8种DNA序列中GG、GT的相邻的状况比较明显；各个DNA序列中碱基丰度比、含量差不多且都含量比较高；其中，DNA序列中、含量差不多且都含量比较低。 对于问题二我们首先通过对问题一散点图的分析选取以碱基的丰度和碱基间的丰度之比为分类的指标，构建为分类的特征向量，但这些特征向量之间存在着一定的相关性，我们运用R型聚类选择出相关性程度差的特征向量为Q型聚类的指标。通过Q型聚类我们将这8种DNA序列分为3种分类方式，通过利用means方法，检验各类别在所有变量上的差异，再利用

3、单因素方差分析最终确定将8种DNA序列分为四类。分类结果如下： 第一类：Human 、Mouse； 第二类：Goat、Rabbit； 第三类：Opossum、 Lemur 、Rat ； 第四类：Gallus。关键词：数理统计；R型聚类；Q型聚类；means法；单因素方差分析法1 问题重述DNA 序列是由A，T，G，C 四个表示4 种碱基的字符组成的序列。研究DNA 序列的结构及序列中隐藏的规律，成为生物信息学的重要研究课题。根据表1 中八个物种的-球蛋白基因的第一个外显子序列，请解决以下问题： 1.建立数学模型刻画序列间的差异； 2. 对表1中八个物种的DNA序列进行分类。2 基本假设1）假设

4、所给的DNA序列片段中没有断句和标点符号；2）假设具有特殊碱基的DNA序列中，特殊碱基可以剔除，其影响可以忽略；3）8个物种DNA序列具有共同的特征；4） 假设给定的DNA序列均是从全序列中随机截取出来的,无法确定序列的起始位, 无法从序列中辨认出氨基酸，所以，在对DNA 序列分类时,从碱基层次上进行分类, 而不是从氨基酸层次上分类；5） 不考虑碱基序列的编码区和非编码区的区别；6） 题目中所给的样本信息量足够大；7） 题目附录中所给的数据真实可靠。3 符号说明：各个DNA序列中碱基出现的数量，i为A、T、C或G：第i个DNA序列的总碱基数目：各个DNA序列中碱基的丰度，i为A、T、C或G：各

5、个DNA序列中碱基i和碱基j的比值，i，j为A、T、C或G ：DNA序列中A、C、G、T的重复次数矩阵：DNA序列中A、C、G、T的所占百分量矩阵：第i个DNA序列相邻碱基占序列相邻情况的百分比，为A、C、T或G：R型聚类的特征向量：DNA序列中四个碱基之间丰度比矩阵 4 模型的建立及求解（1） 问题一模型的建立及求解1） 问题分析首先对数据运用数理统计方法对数据进行计算，得到八个物种的DNA序列的碱基的丰度、碱基的重复出现情况、碱基之间的相邻情况、不同碱基的丰度之比（如碱基A与碱基T的丰度之比）的四个特征，通过对这四个特征作出相应的散点图比较得出八个物种的DNA序列间的差异。 2）模型建立及

6、求解（1） 碱基重复出现的情况 运用matlab求出8种物种DNA序列各自的碱基的重复出现的结果（即每种DNA序列中碱基的个数）和每种DNA序列的碱基数目（即序列的长度）。（matlab运算的程序代码见附录一）其运算的结果如下： 17 21 19 3517 17 17 3521 22 20 2919 15 23 3419 23 15 3517 23 20 3417 20 16 37 20 21 18 33 Human 、Opossum、 Lemur 、Rat等4种DNA序列的长度相同，其他四种DNA序列的长度各不相同；同时每种DNA序列四种碱基的的重复情况也各不相同，其中，Human、Goat

7、、Mouse、Rabbit碱基A的重复情况一样；Gallus、Lemur碱基A重复情况一致；Lemur、Mouse碱基T的重复情况一致；Opossum、Mouse碱基C的重复情况一致；Human、Goat、Lemur碱基G的重复情况一致；Gallus、Mouse碱基的重复情况一样；其他物种碱基重复情况各不相同。 （2）碱基的丰度 对8种DNA序列碱基丰度的分析，i中A碱基丰度的计算： （4-1） 其他碱基T、C、G运算方式一样。通过matlab计算出8种序列的中A、T、C、G四种碱基的丰度结果如下（matlab运算的程序代码见附录一）： 0.1848 0.2283 0.2065 0.3804

8、0.1977 0.1977 0.1977 0.4070 0.2283 0.2391 0.2174 0.3152 0.2088 0.1648 0.2527 0.3736 0.2065 0.2500 0.1630 0.3804 0.1809 0.2447 0.2128 0.3617 0.1889 0.2222 0.1778 0.4111 0.2174 0.2283 0.1957 0.3587 并运用matlab作出8种DNA序列四种碱基丰度的散点图（matlab运算程序代码见附录二）如图4-1所示。图4-1 4种碱基的丰度散点图 通过上述散点图可知每种序列的碱基丰度各有不同，G碱基的丰度相对于本序

9、列的其他碱基的丰度都要高，碱基A的丰度在各个序列中丰度差不多，其他三种碱基在序列中波动性较大，差异性较大。（3） 碱基之间的相邻情况 运用matlab计算出DNA序列相邻碱基的情况，分别为各个序列的AA、AC、AG、AT、CA、CC、CG、CT、GA、GC、GG、GT、TA、TC、TG、TT的相邻次数占各条序列相邻情况的百分比，即如表4-1格式，运用matlab计算DNA序列相邻碱基占序列相邻情况的百分比结果如下（matlab运算程序代码见附录二）：表4-1相邻碱基在序列的排列情况碱基A碱基C碱基G碱基T碱基AAA ACAGAT碱基CCACCCGCT碱基GGAGCGGGT碱基TTATCTGTT

10、 DNA序列相邻碱基占序列相邻情况的百分比： 0.0471 0.0471 0.0824 0.0235 0.0353 0.0824 0.0235 0.0824 0.0824 0.0706 0.1412 0.1059 0.0235 0.0235 0.1647 0.0353 0.0588 0.0235 0.0941 0.0235 0.0353 0.0471 0.0235 0.0941 0.0941 0.1059 0.1412 0.0588 0 0.0235 0.1529 0.0235 0.0353 0.0824 0.0824 0.0353 0.0706 0.0471 0 0.1059 0.0941

11、0.0471 0.0941 0.0706 0.0235 0.0471 0.1294 0.0353 0.0588 0.0353 0.0706 0.0471 0.0588 0.0706 0.0353 0.0824 0.0706 0.1059 0.1412 0.0471 0.0118 0.0353 0.1294 0 0.0471 0.0235 0.0941 0.0471 0.0353 0.0235 0.0118 0.0941 0.1059 0.0706 0.1059 0.0824 0.0118 0.0471 0.1529 0.0471 0.0588 0.0353 0.0588 0.0353 0.02

12、35 0.0824 0.0118 0.0941 0.0941 0.0706 0.0824 0.0941 0 0.0353 0.1882 0.0353 0.0588 0.0118 0.0941 0.0353 0.0471 0.0588 0.0118 0.0588 0.0824 0.0588 0.1294 0.1294 0 0.0471 0.1647 0.0118 0.0706 0.0471 0.0588 0.0471 0.0118 0.0706 0.0118 0.1059 0.0824 0.0824 0.1059 0.0706 0.0471 0 0.1647 0.0235 分析DNA序列相邻碱基

13、占序列相邻情况可知，Human序列中GG、GT、TG的含量较多；Goat序列中出现GG、TG、GC的含量较多，其中TA为0；Opossum序列中出现TG、CT的含量较多，其中CG为0；Gallus序列中出现GC、GG、TG的含量较多，其中TT为0；Lemur序列中出现GA、GG、TG的含量较多；Rabbit序列中出现GT、GG、TG的含量较多，其中TA为0;Rat序列中出现TG、GG、CT的含量较多，其中TC为0。其中，在8种DNA序列中GG、GT的相邻的状况比较明显。（4）不同碱基的丰度之比 为了比较DNA序列之间的差异，不同碱基的丰度之比，是影响其差异的重要原因之一。样品i中碱基T和碱基A

14、的比值计算： （4-2）碱基C与A、碱基G与A、碱基C与T、碱基G与T、碱基G与C运算方式一样。通过matlab计算出DNA序列中四个碱基之间的丰度之比结果如下(matlab运算程序见附录三)： 1.2353 1.1176 2.0588 0.9048 1.6667 1.8421 1.0000 1.0000 2.0588 1.0000 2.0588 2.0588 1.0476 0.9524 1.3810 0.9091 1.3182 1.4500 0.7895 1.2105 1.7895 1.5333 2.2667 1.4783 1.2105 0.7895 1.8421 0.6522 1.5217

15、 2.3333 1.3529 1.1765 2.0000 0.8696 1.4783 1.7000 1.1765 0.9412 2.1765 0.8000 1.8500 2.3125 1.0500 0.9000 1.6500 0.8571 1.5714 1.8333 运用matlab绘制出DNA序列中四种碱基之间丰度之比的散点图如图4-2所示。图4-2碱基丰度比散点图通过图4-2，进行数据分析可得：各个DNA序列中碱基G和碱基C的丰度之比，碱基G和碱基T的丰度之比，碱基C和碱基T的丰度之比含量差不多且都含量比较高；其中，DNA序列中碱基T和碱基A的丰度之比，碱基C和碱基A的丰度之比，碱基G和碱

16、基A的丰度之比含量差不多且都含量比较低。综上八个物种的DNA序列的碱基的丰度、碱基的重复出现情况、碱基之间的相邻情况、不同碱基的丰度之比的四个特征，通过对这四个特征作出相应的散点图比较得出八个物种的DNA序列间的差异。（二）问题二模型的建立及求解1、问题分析 为了使DNA序列的分类能够尽量科学合理，集中要解决的问题是让分类后的样品满足：同类样品间的差异性尽可能小，不同类样品间的差异性尽可能大。为达到上述目的，引入聚类分析模型对不同的DNA序列进行分类。首先我们分析DNA的序列结构，提取出相应的特征。我们分析DNA特征主要从碱基的丰度、不同碱基的丰度之比方面进行入手。我们把8种序列的DNA的碱基

17、丰度和不同碱基丰度作为分类的特征。但这些特征之间存在着一定的相关性，采用R型聚类从中选取代表性的指标。再以这些代表型的指标，采用Q型聚类，对8种DNA序列进行分类。2、 模型的建立与求解（1）DNA序列特征变量的引入 从问题一我们可知影响DNA序列的差异有DNA序列的碱基的丰度、碱基的重复出现情况、碱基之间的相邻情况、不同碱基的丰度之比，经过比较数据的差异我们选择序列中A，C，T，G的碱基丰度和序列中A，C，T，G的碱基间的丰度比作为对8种DNA序列分类的指标。我们以、为指标建立相应的特征向量。0.1848 0.2283 0.2065 0.3804 1.2353 1.1176 2.0588 0

18、.9048 1.6667 1.84210.1977 0.1977 0.1977 0.4070 1.0000 1.0000 2.0588 1.0000 2.0588 2.05880.2283 0.2391 0.2174 0.3152 1.0476 0.9524 1.3810 0.9091 1.3182 1.45000.2088 0.1648 0.2527 0.3736 0.7895 1.2105 1.7895 1.5333 2.2667 1.47830.2065 0.2500 0.1630 0.3804 1.2105 0.7895 1.8421 0.6522 1.5217 2.33330.180

19、9 0.2447 0.2128 0.3617 1.3529 1.1765 2.0000 0.8696 1.4783 1.70000.1889 0.2222 0.1778 0.4111 1.1765 0.9412 2.1765 0.8000 1.8500 2.31250.2174 0.2283 0.1957 0.3587 1.0500 0.9000 1.6500 0.8571 1.5714 1.8333 (2) R型聚类分析 可以看出某些特征之间确实存在着一定的相关性，因此可以考虑从这些指标中选取几个有代表性的指标进行聚类分析。为此，把10个指标根据其相关性进行R型聚类，再从每个类中选取代表性的

20、指标。运用matlab对每个指标的数据分别进行标准化处理，指标间相近性度量采用相关系数，类间相似性度量的计算选用类平均法，用matlab处理数据得到聚类树形图如图4-3所示。图4-3指标聚类树形图 从聚类图4-3中可以看出，碱基A的丰度、碱基C的丰度、碱基G的丰度、碱基T与碱基A的丰度比、碱基G与碱基A的风度比等5个指标之间有较大的相关性，最先聚到一起。如果把10个指标分为6个类，其他5个指标各自为一类。这样就从10个指标中选定了6个分析指标。 碱基T的丰度、碱基C与碱基A的丰度比、碱基G与碱基A的丰度比、碱基C与碱基T的丰度比、碱基G与碱基T的风度比、碱基G与碱基C的风度比 可根据这6个指标

21、对8个物种DNA序列进行聚类分析。（3）Q型聚类分析 根据这6个指标对8个物种DNA的序列进行聚类分析。运用matlab对这6个指标数据分别进行标准化处理，样本间相似性采用欧氏距离度量，类间距离的计算选用类平均法。运用matlab处理数据得到聚类树型图如图4-4所示。表4-2标准化后数据碱基T的丰度碱基C与A丰度比碱基G与A丰度比碱基C与T丰度比碱基G与T丰度比碱基G与C丰度比Human0.22810.73220.9982-0.1396-0.1139-0.0546Goat-0.8603-0.0747-0.00700.22811.06590.5322Opossum0.6122-0.4025-0.

22、9091-0.1215-1.2530-1.2621Gallus-2.03051.37000.74432.28171.7182-1.1781Lemur0.9999-1.5228-1.2636-1.1118-0.61631.3427Mouse0.81141.13701.4619-0.2740-0.7525-0.5259Rabbit0.0111-0.4788-0.0838-0.54120.41251.2790Rat0.2281-0.7606-0.9400-0.3217-0.4609-0.1333图4-4物种DNA聚类树型图（4）模型结果分析 由matlab软件运行的结果可知：8种物种的DNA序列存在

23、着较大的差异。如果根据DNA序列的碱基T的丰度、碱基C与碱基A的丰度比、碱基G与碱基A的丰度比、碱基C与碱基T的丰度比、碱基G与碱基T的丰度比、碱基G与碱基C的丰度比等6个指标可以把8种物种的DNA序列分为三类，结果为：第一类：Opossum、 Lemur、Rat；第二类：Human、Goat、Mouse、Rabbit；第三类：Gallus。如果根据这6个指标把8种物种的DNA序列分为四类，结果为：第一类：Human 、Mouse；第二类：Goat、Rabbit；第三类：Opossum、 Lemur 、Rat ；第四类：Gallus。如果根据这6个指标把8种物种的DNA序列分为五类，结果为：第

24、一类： Lemur；第二类：Opossum、Rat；第三类：Human、Mouse；第四类：Goat、Rabbit；第五类：Gallus。（2） 、运用spss对分类的结果进行检验 运用spss以上述6个指标利用means方法，检验各类别在所有变量上的差异，如果差异显著，我们认为分类结果是可靠的，我们利用spss算出三种分类结果各自的平均数（三种分类结果各自的平均数见附表五）再对8种物种DNA序列进行分类的结果单因素方差分析，判断分类的结果如表4-3所示。表4-3（a)序列分三类方差分析表表4-3（b）序列分四类方差分析表表4-3（c）序列分五类方差分析表方差分析结果显示，8种DNA序列分三类

25、时，DNA序列的碱基T的丰度、碱基G与碱基A的丰度比、碱基G与碱基T的丰度比等三项指标达到显著水平，说明这种分类效果不是很好。8种DNA序列分四类时，DNA序列的碱基T的丰度、碱基C与碱基A的丰度比、碱基G与碱基A的丰度比、碱基C与碱基T的丰度比、碱基G与碱基T的丰度比等五项指标达到显著水平，说明这种分类效果是比较明显的。8种DNA序列分五类时，碱基C与碱基A的丰度比、碱基G与碱基A的丰度比、碱基C与碱基T的丰度比等三项指标达到显著水平，说明这种分类效果不是很好。综上分析可知我们最终确定将8种物种DNA序列分为四类，分别为： 第一类：Human 、Mouse； 第二类：Goat、Rabbit；

26、 第三类：Opossum、 Lemur 、Rat ； 第四类：Gallus。5 模型的评价1） 模型的优点（1）、问题一的求解简单容易，易于理解，运用数理统计的方法转化原问题并且求出各个序列所具有的统计特征。（2）、在对于问题二我们首先在对问题一分析的基础上对碱基的丰度，和碱基之间的碱基比散点图，根据其波动的大小初步判断出那些变量作为聚类分析的特征向量，为后续的分析缩短了数据处理的过程，简化了模型。（3）、对于问题二我们采用R型聚类分析我们对分类的指标进行了相关性分析，选择了特征之间一些关联性差的特征作为Q型聚类分析的指标，方法简单易行，有较好的普遍性。（4）、对于问题二，在R型聚类类的基础上

27、我们建立了Q型聚类方法得到了三种分类方式，考虑到了一些因素，最终结果更加真实可靠。（5）、在问题二时，使用了SPSS对三种分类的结果进行了单因素方差分析，从而得到分类的结果更合理，更可靠。（6）、对于问题二，从R聚类分析、Q型聚类分析、means方法、到单因素方差分析，从而使得分类的结果比较切合实际的。2） 模型的缺点（1）、在进行模型一和模型二的最终选择时，综合考虑的因素还有欠缺。（2）、在问题二的基础上建立的聚类分类模型，考虑的指标不够全面。（3）、在建立模型的时候对数据处理的不合理，不是很到位，对模型的结果有 一定的影响，模型的结果与实际结果之间存在着一定的差距。（4）、在用MATLAB

28、编程时所编写的代码复杂程度不够好，致使在考虑DNA序列样本更多的时候，代码不太具有参考性，带有局限性。（5）、对于问题一分析DNA序列差异的指标不够全面，导致对DNA序列差异分析的不到位，同时问题二选择的指标也存在着比较片面。（6）、问题一采用数理分析过于简单未能将DNA序列之间的差异很好的表现出来，处理数据方面处理的不是很理想。6 模型的改进 DNA序列问题我们考虑了碱基的丰度特征值等，而一个序列所含的信息远不止每个碱基的丰度等特征,还有基于碱基所在位置的有关特征，即碱基在序列中出现的规律性、碱基和它前后若干个碱基的相关性等等，我们可将DNA序列中碱基的排列看成是一个随机过程,如果是着重研究

29、随机过程发生的规律,并设碱基分别为：A=0,T=1,C=2,G=3,那么这些序列就可以看作一组离散的数字信号,则可以用数字信号处理的理论来进行研究与分析即将离散的数字信号转化为一组波形图，用周期性、数学期望、相关函数以及频率谱等方法描述其规律,这当然就增加了问题的复杂性,但就一个庞大的DNA序列规律破解工程来说,不失为一种值得考虑的方法，这里,可以仅以周期性为例讨论如下：对某个碱基,以a为例,假设它在序列中的t1,t2,tk+1个位置出现,可试图找出这些碱基之间的关联。首先,可以认识到考察ti的分布及绝对值意义不大,因为序列是一大段DNA中的一个片段,片段的起始段不同会导致ti的不同,于是考虑

30、a的间距：si=ti+1-ti (i=1,2,k)1可以看出,序列s1,s2,sn的大小包含的信息是a的“稠密度”,也可以看成一个与频率有关的量,这个在求解问题一时已讨论过，而在这我们考察序列s1,s2,sn的波动幅度,幅度越小,说明si(i=1,2,k)的值越趋于统一,即a的周期性越大1而表征波动幅度的量在统计中是中心矩，现求si的二阶中心矩,即方差: （6-1） （6-2）同理可求Varg,Vart,Varc1由于a,t,c,g成对出现,作判别函数F=Varg/Vart可较好的对序列进行分类。7 参考文献1 汤诗杰，周亮，王晓玲DNA序列分类模型J预防医学情报杂志，2005，(6

31、): 83-852 陈合格三种鳖线粒体DNA部分基因序列的比较分析和分子鉴定标记J工业工程，2006，(5):23-273 罗贤晖，江从喜，洪翔. 基于神经网络集成的DNA序列分类方法研究J. 中国市场，2012，(49):76-774 孙晓敏，张厚粲聚类分析方法在DNA序列分类中的应用J心理科学2005，（3）：646-649.5 姜启源数学模型M北京：高等教育出版社，19936 姜启源.邢文训.谢金星杨顶辉.大学数学实验.清华大学出版社7 陈东彦，李冬梅，王树忠数学建模北京.20098 张德丰. MATLAB概率数理统计分析M.北京：机械工业出版社，2010附录：（1） 碱基重复出现的情况

32、及各碱基丰度matlab程序 %输入数据A1='ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATTAAGTTGGTGGTGAGGCCCTGGGCAG'A2='ATGCTGACTGCTGAGGAGAAGGCTGCCGTCACCGGCTTCTGGGGCAAGGTGAAAGTGGATGAAGTTGGTGCTGAGGCCCTGGGCAG' A3='ATGGTGCACTTGACTTCTGAGGAGAAGAACTGCATCACTACCATCTGGTCTAAGGTGCAGGTTGACC

33、AGACTGGTGGTGAGGCCCTTGGCAG' A4='ATGGTGCACTGGACTGCTGAGGAGAAGCAGCTCATCACCGGCCTCTGGGGCAAGGTAATGTGGCCGAATGTGGGGCCGAAGCCCTGGCCAG' A5='ATGACTTTGCTGAGTGCTGAGGAGAATGCTCATGTCACCTCTCTGTGGGGCAAGGTGGATGTAGAGAAAGTTGGTGGCGAGGCCTTGGGCAG'A6='ATGGTTGCACCTGACTGATGCTGAGAAGTCTGCTGTCTCTTGCCTGTGGGC

34、AAAGGTGAACCCCGATGAAGTTGGTGGTGAGGCCCTGGGCAGG'A7='ATGGTGCATCTGTCCAGTGAGGAGAAGTCTGCGGTCACTGCCCTGTGGGGCAAGGTGAATGTGGAAGAAGTTGGTGGTGAGGCCCTGGGC'A8='ATGGTGCACCTAACTGATGCTGAGAAGGCTACTGTTAGTGGCCTGTGGGGAAAGGTGAACCCTGATAATGTTGGCGCTGAGGCCCTGGGCAG'%求单个碱基丰度及丰度之比x1为A、T、C、G的重复次数x2为A、T、C、G的丰度x1=

35、zeros(8,4);x2=zeros(8,4); u,y=size(A1);for j=1:y m=A1(j); switch m case 'A' x1(1,1)=x1(1,1)+1; case 'T' x1(1,2)=x1(1,2)+1; case 'C' x1(1,3)=x1(1,3)+1; case 'G' x1(1,4)=x1(1,4)+1; endendx2(1,:)=x1(1,:)/y;u,y=size(A2);for j=1:y m=A2(j); switch m case 'A' x1(2,1)

36、=x1(2,1)+1; case 'T' x1(2,2)=x1(2,2)+1; case 'C' x1(2,3)=x1(2,3)+1; case 'G' x1(2,4)=x1(2,4)+1; endendx2(2,:)=x1(2,:)/y;u,y=size(A3);for j=1:y m=A3(j); switch m case 'A' x1(3,1)=x1(3,1)+1; case 'T' x1(3,2)=x1(3,2)+1; case 'C' x1(3,3)=x1(3,3)+1; case &#

37、39;G' x1(3,4)=x1(3,4)+1; endendx2(3,:)=x1(3,:)/y;u,y=size(A4);for j=1:y m=A4(j); switch m case 'A' x1(4,1)=x1(4,1)+1; case 'T' x1(4,2)=x1(4,2)+1; case 'C' x1(4,3)=x1(4,3)+1; case 'G' x1(4,4)=x1(4,4)+1; endendx2(4,:)=x1(4,:)/y;u,y=size(A5);for j=1:y m=A5(j); switch

38、 m case 'A' x1(5,1)=x1(5,1)+1; case 'T' x1(5,2)=x1(5,2)+1; case 'C' x1(5,3)=x1(5,3)+1; case 'G' x1(5,4)=x1(5,4)+1; endendx2(5,:)=x1(5,:)/y;u,y=size(A6);for j=1:y m=A6(j); switch m case 'A' x1(6,1)=x1(6,1)+1; case 'T' x1(6,2)=x1(6,2)+1; case 'C'

39、 x1(6,3)=x1(6,3)+1; case 'G' x1(6,4)=x1(6,4)+1; endendx2(6,:)=x1(6,:)/y;u,y=size(A7);for j=1:y m=A7(j); switch m case 'A' x1(7,1)=x1(7,1)+1; case 'T' x1(7,2)=x1(7,2)+1; case 'C' x1(7,3)=x1(7,3)+1; case 'G' x1(7,4)=x1(7,4)+1; endx2(7,:)=x1(7,:)/y;endu,y=size(A8

40、);for j=1:y m=A8(j); switch m case 'A' x1(8,1)=x1(8,1)+1; case 'T' x1(8,2)=x1(8,2)+1; case 'C' x1(8,3)=x1(8,3)+1; case 'G' x1(8,4)=x1(8,4)+1; endendx2(8,:)=x1(8,:)/y;%画散点图 subplot(2,2,1) plot(1:8,x2(1:8,1),'*') xlabel('八种DNA序列') ylabel('A碱基丰度')

41、 title('单个碱基丰度') axis(1,8,0.15,0.45) hold on subplot(2,2,2) plot(1:8,x2(1:8,2),'*') xlabel('八种DNA序列') ylabel('T碱基丰度') axis(1,8,0.15,0.45) hold on subplot(2,2,3) plot(1:8,x2(1:8,3),'*') xlabel('八种DNA序列') ylabel('C碱基丰度') axis(1,8,0.15,0.45) hold

42、on subplot(2,2,4) plot(1:8,x2(1:8,4),'*') xlabel('八种DNA序列') ylabel('G碱基丰度') axis(1,8,0.15,0.45)（二）相邻碱基序列情况matlab程序%求相邻碱基序列情况A1_XX、A2_XX、A3_XX、A4_XX、A5_XX、A6_XX、A7_XX、A8_XX分别为十六种碱基相邻情况的分析结果%以求A1样本为例，其余求解代码只需将A1分别改为A2、A3、A4、A5、A6、A7、A8即可char m;char n; Yaa=0;Yag=0;Yat=0;Yac=0; Y

43、ca=0;Ycc=0;Ycg=0;Yct=0; Yga=0;Ygc=0;Ygg=0;Ygt=0; Yta=0;Ytc=0;Ytg=0;Ytt=0; u,y=size(A2);for j=1:y-1;m=A1(j);n=A1(j+1);switch m;case 'A' switch n; case 'A' Yaa=Yaa+1; case 'C' Yac=Yac+1; case 'G' Yag=Yag+1; case 'T' Yat=Yat+1; endcase 'C' switch n; case

44、'A' Yca=Yca+1; case 'C' Ycc=Ycc+1; case 'G' Ycg=Ycg+1; case 'T' Yct=Yct+1; endcase 'G' switch n; case 'A' Yga=Yga+1; case 'C' Ygc=Ygc+1; case 'G' Ygg=Ygg+1; case 'T' Ygt=Ygt+1; endcase 'T' switch n case 'A' Yta=Y

45、ta+1; case'C' Ytc=Ytc+1; case 'G' Ytg=Ytg+1; case 'T' Ytt=Ytt+1; endendendA1_r=Yaa,Yac,Yag,Yat; Yca,Ycc,Ycg,Yct; Yga,Ygc,Ygg,Ygt; Yta,Ytc,Ytg,Ytt Y=Yaa+Yac+Yag+Yat+Yca+Ycc+Ycg+Yct+Yga+Ygc+Ygg+Ygt+Yta+Ytc+Ytg+Ytt;A1_XX=A1_X/Y（三）求两个碱基间丰度之比（，）表示/、/、/、/、/、/p=zeros(8,6);for i=1:8;

46、 p(i,1)=x1(i,2)/x1(i,1) p(i,2)=x1(i,3)/x1(i,1) p(i,3)=x1(i,4)/x1(i,1) p(i,4)=x1(i,3)/x1(i,2) p(i,5)=x1(i,4)/x1(i,2) p(i,6)=x1(i,4)/x1(i,3)End%画出相应的散点图 subplot(2,3,1) plot(1:8,x2(1:8,1),'+') xlabel('八种DNA序列') ylabel('T/A碱基丰度之比') title('碱基对丰度之比') axis(1,8,0,0.7) hold on

47、 subplot(2,3,2) plot(1:8,x2(1:8,2),'+') xlabel('八种DNA序列') ylabel('C/A碱基丰度之比') axis(1,8,0,0.7) hold on subplot(2,3,3) plot(1:8,x2(1:8,3),'+') xlabel('八种DNA序列') ylabel('G/A碱基丰度之比') axis(1,8,0,0.7) hold on subplot(2,3,4) plot(1:8,x2(1:8,4),'+') xl

48、abel('八种DNA序列') ylabel('C/T碱基丰度之比') axis(1,8,0,0.7) subplot(2,3,5) plot(1:8,x2(1:8,4),'+') xlabel('八种DNA序列') ylabel('G/T碱基丰度之比') axis(1,8,0,0.7) subplot(2,3,6) plot(1:8,x2(1:8,4),'+') xlabel('八种DNA序列') ylabel('G/C碱基丰度之比') axis(1,8,0,0.7

49、)（四）R型聚类分析matlab程序 clear>> X=0.1848 0.2283 0.2065 0.3840 1.2354 1.1174 2.0588 0.9045 1.6820 1.8596 0.1977 0.1977 0.1977 0.4070 1.0000 1.0000 2.0588 1.0000 2.0587 2.0587 0.2283 0.2391 0.2174 0.3152 1.0473 0.9523 1.3810 0.9092 1.3183 1.4499 0.2088 0.1648 0.2527 0.3736 0.7893 1.2102 1.7895 1.5334

50、 2.2670 1.4784 0.2065 0.2500 0.1630 0.3804 1.2107 0.7893 1.8421 0.6520 1.5216 2.3337 0.1809 0.2447 0.2128 0.3617 1.3527 1.1763 2.0000 0.8696 1.4781 1.6997 0.1889 0.2222 0.1778 0.4111 1.1763 0.9412 2.1765 0.8002 1.8501 2.3121 0.2174 0.2283 0.1957 0.3587 1.0501 0.9002 1.6500 0.8572 1.5712 1.8329 ;>> d=pdist(X');z=linkage(d,'average'