DNA序列问题模型详解

2015年芜湖三校数学建模竞赛

题目 DNA序列问题模型

摘要

DNA 序列是由A，T，G，C 四个表示4 种碱基的字符组成的序列。本文研究DNA 序列的结构找出序列间的差异和对八个物种的DNA序列进行分类。

对于问题一首先对数据运用数理统计方法对数据进行计算，得到八个物种的DNA序列的碱基的丰度、碱基的重复出现情况、碱基之间的相邻情况、不同碱

基的丰度之比的四个特征，通过对这四个特征作出相应的散点图比较得出八个物种的DNA序列间的差异：Human 、Opossum、 Lemur 、Rat等4种DNA序列的长度相同，其他四种DNA序列的长度各不相同，每种DNA序列四种碱基的的

重复情况也各不相同；G碱基的丰度相对于本序列的其他碱基的丰度都要高，碱基A的丰度在各个序列中丰度差不多，其他三种碱基在序列中波动性较大，差异性较大；8种DNA序列中GG、GT的相邻的状况比较明显；各个DNA序列中碱基

丰度比

f、GT f、CT f含量差不多且都含量比较高；其中，DNA序列中TA f、CA f、

GC

f含量差不多且都含量比较低。

GA

对于问题二我们首先通过对问题一散点图的分析选取以碱基的丰度和碱基间的丰度之比为分类的指标，构建为分类的特征向量，但这些特征向量之间存在着一定的相关性，我们运用R型聚类选择出相关性程度差的特征向量为Q型聚类的指标。通过Q型聚类我们将这8种DNA序列分为3种分类方式，通过利用means方法，检验各类别在所有变量上的差异，再利用单因素方差分析最终确定将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）假设所给的DNA序列片段中没有断句和标点符号；

2）假设具有特殊碱基的DNA序列中，特殊碱基可以剔除，其影响可以忽略；

3）8个物种DNA序列具有共同的特征；

4）假设给定的DNA序列均是从全序列中随机截取出来的,无法确定序列的起始位, 无法从序列中辨认出氨基酸，所以，在对DNA 序列分类时,从碱基层次上进行分类, 而不是从氨基酸层次上分类；

5）不考虑碱基序列的编码区和非编码区的区别；

6）题目中所给的样本信息量足够大；

7）题目附录中所给的数据真实可靠。

3 符号说明

n：各个DNA序列中碱基i出现的数量，i为A、T、C或G

i

N：第i个DNA序列的总碱基数目

i

F：各个DNA序列中碱基i的丰度，i为A、T、C或G

i

f：各个DNA序列中碱基i和碱基j的比值，i，j为A、T、C或G

ij

X：DNA序列中A、C、G、T的重复次数矩阵

1

X：DNA序列中A、C、G、T的所占百分量矩阵

2

XX：第i个DNA序列相邻碱基占序列相邻情况的百分比，X为A、C、T或G i

Y：R型聚类的特征向量

Z：DNA序列中四个碱基之间丰度比矩阵

4 模型的建立及求解

（一）问题一模型的建立及求解

1）问题分析

首先对数据运用数理统计方法对数据进行计算，得到八个物种的DNA序列的碱基的丰度、碱基的重复出现情况、碱基之间的相邻情况、不同碱基的丰度之比（如碱基A与碱基T的丰度之比）的四个特征，通过对这四个特征作出相应的散点图比较得出八个物种的DNA序列间的差异。

2）模型建立及求解

（1）碱基重复出现的情况

运用matlab 求出8种物种DNA 序列各自的碱基的重复出现的结果（即每种DNA 序列中碱基的个数）和每种DNA 序列的碱基数目（即序列的长度）。（matlab 运算的程序代码见附录一）其运算的结果如下： =1X [

17 21 19 35

17 17 17 35 21 22 20 29 19 15 23 34 19 23 15 35 17 23 20 34 17 20 16 37

20 21 18 33 ]

921=N 862=N 923=N 914=N 925=N 946=N 907=N

928=N

Human 、Opossum 、 Lemur 、Rat 等4种DNA 序列的长度相同，其他四种DNA 序列的长度各不相同；同时每种DNA 序列四种碱基的的重复情况也各不相同，其中，Human 、Goat 、Mouse 、Rabbit 碱基A 的重复情况一样；Gallus 、Lemur 碱基A 重复情况一致；Lemur 、Mouse 碱基T 的重复情况一致；Opossum 、Mouse 碱基C 的重复情况一致；Human 、Goat 、Lemur 碱基G 的重复情况一致；Gallus 、Mouse 碱基的重复情况一样；其他物种碱基重复情况各不相同。 （2）碱基的丰度

对8种DNA 序列碱基丰度的分析，i 中A 碱基丰度的计算： /A A i F n N = （4-1）

其他碱基T 、C 、G 运算方式一样。通过matlab 计算出8种序列的中A 、T 、C 、G 四种碱基的丰度结果如下（matlab 运算的程序代码见附录一）：

=2X [

0.1848 0.2283 0.2065 0.3804 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碱基的丰度相对于本序列的其他碱基的丰度都要高，碱基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相邻碱基在序列的排列情况

DNA序列相邻碱基占序列相邻情况的百分比：

XX[

1

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 ]

XX[

=

2

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 ]

XX[

=

3

0.0353 0.0824 0.0824 0.0353

0.0706 0.0471 0 0.1059

0.0941 0.0471 0.0941 0.0706

0.0235 0.0471 0.1294 0.0353 ]

XX[

=

4

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 ]

XX[

=

5

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 ]

XX[

=

6

0.0588 0.0353 0.0588 0.0353

0.0235 0.0824 0.0118 0.0941

0.0941 0.0706 0.0824 0.0941

0 0.0353 0.1882 0.0353 ]

XX[

=

7

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 ]

XX[

=

8

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 序列相邻碱基占序列相邻情况可知，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 的比值计算：

A T TA n n f /= （4-2）

碱基C 与A 、碱基G 与A 、碱基C 与T 、碱基G 与T 、碱基G 与C 运算方式一样。通过matlab 计算出DNA 序列中四个碱基之间的丰度之比结果如下(matlab 运算程序见附录三)：

=Z [

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 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和碱基A的丰度之比含量差不多且都含量比较低。

综上八个物种的DNA序列的碱基的丰度、碱基的重复出现情况、碱基之间的相邻情况、不同碱基的丰度之比的四个特征，通过对这四个特征作出相应的散点图比较得出八个物种的DNA序列间的差异。

（二）问题二模型的建立及求解

1、问题分析

为了使DNA序列的分类能够尽量科学合理，集中要解决的问题是让分类后的样品满足：同类样品间的差异性尽可能小，不同类样品间的差异性尽可能大。为达到上述目的，引入聚类分析模型对不同的DNA序列进行分类。首先我们分析DNA的序列结构，提取出相应的特征。我们分析DNA特征主要从碱基的丰度、不同碱基的丰度之比方面进行入手。我们把8种序列的DNA的碱基丰度和不同碱基丰度作为分类的特征。但这些特征之间存在着一定的相关性，采用R型聚类从中选取代表性的指标。再以这些代表型的指标，采用Q型聚类，对8种DNA序列进行分类。

2、模型的建立与求解

（1）DNA序列特征变量的引入

从问题一我们可知影响DNA序列的差异有DNA序列的碱基的丰度、碱基的重复出现情况、碱基之间的相邻情况、不同碱基的丰度之比，经过比较数据的差异我们选择序列中A，C，T，G的碱基丰度和序列中A，C，T，G的碱基间的

丰度比作为对8种DNA 序列分类的指标。我们以A F 、T F 、C F 、G F 、TA f 、CA f 、

GA f 、CT f 、GT f 、GC f 为指标建立相应的特征向量Y 。

Y [

0.1848 0.2283 0.2065 0.3804 1.2353 1.1176 2.0588 0.9048 1.6667 1.8421

0.1977 0.1977 0.1977 0.4070 1.0000 1.0000 2.0588 1.0000 2.0588 2.0588

0.2283 0.2391 0.2174 0.3152 1.0476 0.9524 1.3810 0.9091 1.3182 1.4500

0.2088 0.1648 0.2527 0.3736 0.7895 1.2105 1.7895 1.5333 2.2667 1.4783

0.2065 0.2500 0.1630 0.3804 1.2105 0.7895 1.8421 0.6522 1.5217 2.3333

0.1809 0.2447 0.2128 0.3617 1.3529 1.1765 2.0000 0.8696 1.4783 1.7000

0.1889 0.2222 0.1778 0.4111 1.1765 0.9412 2.1765 0.8000 1.8500 2.3125

0.2174 0.2283 0.1957 0.3587 1.0500 0.9000 1.6500 0.8571 1.5714 1.8333 ]

(2) R 型聚类分析

可以看出某些特征之间确实存在着一定的相关性，因此可以考虑从这些指标中选取几个有代表性的指标进行聚类分析。为此，把10个指标根据其相关性进行R 型聚类，再从每个类中选取代表性的指标。运用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个指标对8个物种DNA序列进行聚类分析。

（3）Q型聚类分析

根据这6个指标对8个物种DNA的序列进行聚类分析。运用matlab对这6个指标数据分别进行标准化处理，样本间相似性采用欧氏距离度量，类间距离的计算选用类平均法。运用matlab处理数据得到聚类树型图如图4-4所示。

图4-4物种DNA聚类树型图

（4）模型结果分析

由matlab软件运行的结果可知：8种物种的DNA序列存在着较大的差异。如果根据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序列分为五类，结果为：

第一类： 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序列分三类时，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；

第三类：Opossum、 Lemur 、Rat ；

第四类：Gallus。

5 模型的评价

1）模型的优点

（1）、问题一的求解简单容易，易于理解，运用数理统计的方法转化原问题并且求出各个序列所具有的统计特征。

（2）、在对于问题二我们首先在对问题一分析的基础上对碱基的丰度，和碱基之间的碱基比散点图，根据其波动的大小初步判断出那些变量作为聚类分析的特征向量，为后续的分析缩短了数据处理的过程，简化了模型。

（3）、对于问题二我们采用R型聚类分析我们对分类的指标进行了相关性分析，选择了特征之间一些关联性差的特征作为Q型聚类分析的指标，方法简单易行，有较好的普遍性。

（4）、对于问题二，在R型聚类类的基础上我们建立了Q型聚类方法得到了三种分类方式，考虑到了一些因素，最终结果更加真实可靠。

（5）、在问题二时，使用了SPSS对三种分类的结果进行了单因素方差分析，从而得到分类的结果更合理，更可靠。

（6）、对于问题二，从R聚类分析、Q型聚类分析、means方法、到单因素方差分析，从而使得分类的结果比较切合实际的。

2）模型的缺点

（1）、在进行模型一和模型二的最终选择时，综合考虑的因素还有欠缺。（2）、在问题二的基础上建立的聚类分类模型，考虑的指标不够全面。（3）、在建立模型的时候对数据处理的不合理，不是很到位，对模型的结果有一定的影响，模型的结果与实际结果之间存在着一定的差距。

（4）、在用MATLAB编程时所编写的代码复杂程度不够好，致使在考虑DNA序列样本更多的时候，代码不太具有参考性，带有局限性。

（5）、对于问题一分析DNA序列差异的指标不够全面，导致对DNA序列差异分析的不到位，同时问题二选择的指标也存在着比较片面。

（6）、问题一采用数理分析过于简单未能将DNA序列之间的差异很好的表现出来，处理数据方面处理的不是很理想。

6 模型的改进

DNA序列问题我们考虑了碱基的丰度特征值等，而一个序列所含的信息远不止每个碱基的丰度等特征,还有基于碱基所在位置的有关特征，即碱基在序列中出现的规律性、碱基和它前后若干个碱基的相关性等等，我们可将DNA序列中碱

基的排列看成是一个随机过程,如果是着重研究随机过程发生的规律,并设碱基分别为：A=0,T=1,C=2,G=3,那么这些序列就可以看作一组离散的数字信号,则可以用数字信号处理的理论来进行研究与分析即将离散的数字信号转化为一组波形图，用周期性、数学期望、相关函数以及频率谱等方法描述其规律,这当然就增加了问题的复杂性,但就一个庞大的DNA 序列规律破解工程来说,不失为一种值得考虑的方法，这里,可以仅以周期性为例讨论如下：

对某个碱基,以a 为例,假设它在序列中的t1,t2,…,tk+1个位置出现,可试图找出这些碱基之间的关联。首先,可以认识到考察ti 的分布及绝对值意义不大,因为序列是一大段DNA 中的一个片段,片段的起始段不同会导致ti 的不同,于是考虑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 的二阶中心矩,即方差:

∑-=-

-=n

i ara

s s i n v 12

)

(11 （6-1）

∑=-

=n

i i

s n s 1

1

（6-2）

同理可求Varg,Vart,Varc1由于a,t,c,g 成对出现,作判别函数F=Varg/Vart

可较好的对序列进行分类。

7 参考文献

[1] 汤诗杰，周亮，王晓玲．DNA 序列分类模型[J]．预防医学情报杂志，2005，(6): 83-85． [2] 陈合格．三种鳖线粒体DNA 部分基因序列的比较分析和分子鉴定标记[J]．工业工程，2006，(5):23-27． [3] 罗贤晖，江从喜，洪翔. 基于神经网络集成的DNA 序列分类方法研究[J]. 中国市场，2012，(49):76-77． [4] 孙晓敏，张厚粲．聚类分析方法在DNA 序列分类中的应用[J]．心理科学．2005，（3）：646-649.．

[5] 姜启源．数学模型[M]．北京：高等教育出版社，1993．

[6] 姜启源.邢文训.谢金星．杨顶辉.大学数学实验.清华大学出版社．

[7] 陈东彦，李冬梅，王树忠．《数学建模》．北京.2009．

[8] 张德丰. MATLAB概率数理统计分析[M].北京：机械工业出版社，2010．

附录：

（一）碱基重复出现的情况及各碱基丰度matlab程序

%输入数据

A1='ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGG TGAACGTGGATTAAGTTGGTGGTGAGGCCCTGGGCAG';

A2='ATGCTGACTGCTGAGGAGAAGGCTGCCGTCACCGGCTTCTGGGGCAAGGTGAAA GTGGATGAAGTTGGTGCTGAGGCCCTGGGCAG';

A3='ATGGTGCACTTGACTTCTGAGGAGAAGAACTGCATCACTACCA TCTGGTCTAAGG TGCAGGTTGACCAGACTGGTGGTGAGGCCCTTGGCAG';

A4='ATGGTGCACTGGACTGCTGAGGAGAAGCAGCTCATCACCGGCCTCTGGGGCAAG GTAATGTGGCCGAATGTGGGGCCGAAGCCCTGGCCAG';

A5='ATGACTTTGCTGAGTGCTGAGGAGAATGCTCATGTCACCTCTCTGTGGGGCAAGG TGGATGTAGAGAAAGTTGGTGGCGAGGCCTTGGGCAG';

A6='ATGGTTGCACCTGACTGATGCTGAGAAGTCTGCTGTCTCTTGCCTGTGGGCAAAG GTGAACCCCGATGAAGTTGGTGGTGAGGCCCTGGGCAGG';

A7='ATGGTGCATCTGTCCAGTGAGGAGAAGTCTGCGGTCACTGCCCTGTGGGGCAAGG TGAATGTGGAAGAAGTTGGTGGTGAGGCCCTGGGC';

A8='ATGGTGCACCTAACTGATGCTGAGAAGGCTACTGTTAGTGGCCTGTGGGGAAAGG TGAACCCTGATAATGTTGGCGCTGAGGCCCTGGGCAG';

%求单个碱基丰度及丰度之比

x1为A、T、C、G的重复次数

x2为A、T、C、G的丰度

x1=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;

end

end

x2(1,:)=x1(1,:)/y;

[u,y]=size(A2);

for j=1:y

m=A2(j);

switch m

case 'A';

x1(2,1)=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;

end

end

x2(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 'G';

x1(3,4)=x1(3,4)+1;

end

end

x2(3,:)=x1(3,:)/y;

[u,y]=size(A4);

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;

end

end

x2(4,:)=x1(4,:)/y;

[u,y]=size(A5);

for j=1:y

m=A5(j);

switch 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;

end

end

x2(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';

x1(6,3)=x1(6,3)+1;

case 'G';

x1(6,4)=x1(6,4)+1;

end

x2(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;

end

x2(7,:)=x1(7,:)/y;

end

[u,y]=size(A8);

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;

end

end

x2(8,:)=x1(8,:)/y;

%画散点图

subplot(2,2,1)

plot(1:8,x2(1:8,1),'*') xlabel('八种DNA序列') ylabel('A碱基丰度')

title('单个碱基丰度')

axis([1,8,0.15,0.45])

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 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;

Yca=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;

end

case'C'

switch n;

case'A'

Yca=Yca+1;

case'C'

Ycc=Ycc+1;

case'G'

Ycg=Ycg+1;

case'T'

Yct=Yct+1;

end

case'G'

switch n;

case'A'

Yga=Yga+1;

case'C'

Ygc=Ygc+1;

case'G'

Ygg=Ygg+1;

case'T'

Ygt=Ygt+1;

end

case'T'

switch n

case'A'

Yta=Yta+1;

case'C'

Ytc=Ytc+1;

case'G'

Ytg=Ytg+1;

case'T'

Ytt=Ytt+1;

end

end

end

A1_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;

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

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),'+')