07DNA序列分析.docx
《07DNA序列分析.docx》由会员分享，可在线阅读，更多相关《07DNA序列分析.docx（33页珍藏版）》请在冰豆网上搜索。
07DNA序列分析
2014高教社杯全国大学生数学建模竞赛
承诺书
我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。
如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）：
D
我们的参赛报名号为（如果赛区设置报名号的话）：
201507
所属学校（请填写完整的全名）：
南京工程学院
参赛队员（打印并签名）：
1.王星云
2.顾玖
3.周洁
指导教师或指导教师组负责人（打印并签名）：
（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。
以上内容请仔细核对，提交后将不再允许做任何修改。
如填写错误，论文可能被取消评奖资格。
）
日期：
年月日
赛区评阅编号（由赛区组委会评阅前进行编号）：
2014高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
评
阅
人
评
分
备
注
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
DNA序列分类
摘要
随着生物信息的迅猛发展，海量的DNA序列数据被收集汇编在各种数据库中。
将DNA序列进行分类，发现其规律和结构，对认识生命的起源有着重要意义。
本文首先将问题中的A、B两类DNA序列转换成氨基酸序列，然后，针对该序列进行特征提取和选择，从而把研究DNA序列结构转换为分析氨基酸的频率问题。
通过MATLAB统计出20种氨基酸的频率，并将其作为20个初始特征，通过Ward最小方差法，进行谱系聚类。
再采用SPSS软件作出树状图，对202个序列进行分类。
其中对21~40号人工序列号的分类如下：
A类：
2124262831333840
B类：
222325272930323435363739
对于1~182号自然序列，由于单个数据序列过长，我们选取其具有实际意义的前120组数据，计算每组数据的20种氨基酸频率。
因为精氨酸频率与合成防御素密切相关，所以筛选出精氨酸含量前40位的数据进行ward聚类。
对1~182号自然序列中精氨酸含量高的序列的分类如下表：
A类：
155171162711213645141324979159551041513290
1116910970169
B类：
14158117284210010611611351791767114747881
关键词：
谱系聚类；氨基酸的频率；DNA序列分类
一、问题重述
DNA全序列图是由4个字符A，T，C，G按一定顺序排成的长约30亿的序列，其中没有“断句”也没有标点符号，除了这4个字符表示4种碱基以外，人们对它包含的“内容”知之甚少，难以读懂。
研究DNA全序列具有什么结构，由这4个字符排成的看似随机的序列中隐藏着什么规律，是生物信息学最重要的课题之一。
人们也发现了DNA序列中的一些规律性和结构。
例如，在全序列中有一些是用于编码蛋白质的序列片段，即由这4个字符组成的64种不同的3字符串，其中大多数用于编码构成蛋白质的20种氨基酸。
又例如，在不用于编码蛋白质的序列片段中，A和T的含量特别多些，于是以某些碱基特别丰富作为特征去研究DNA序列的结构也取得了一些结果。
此外，利用统计的方法还发现序列的某些片段之间具有相关性，等等。
这些发现让人们相信，DNA序列中存在着局部的和全局性的结构，充分发掘序列的结构对理解DNA全序列是十分有意义的。
目前在这项研究中最普通的思想是省略序列的某些细节，突出特征，然后将其表示成适当的数学对象。
这种被称为粗粒化和模型化的方法往往有助于研究规律性和结构。
作为研究DNA序列的结构的尝试，提出以下对序列集合进行分类的问题：
（1）数据文件中有20个已知类别的人工制造的序列，其中序列标号1—10为A类，11-20为B类。
请从中提取特征，构造分类方法，并用这些已知类别的序列，衡量你的方法是否足够好。
用你认为满意的方法，对另外20个未标明类别的人工序列（标号21—40）进行分类，把结果用序号（按从小到大的顺序）标明它们的类别（无法分类的不写入）：
A类；B类。
请详细描述你的方法，给出计算程序。
如果你部分地使用了现成的分类方法，也要将方法名称准确注明。
（2）在数据文件中给出了182个自然DNA序列，它们都较长。
用你的分类方法对它们进行分类，一样地给出分类结果。
二、问题的分析
问题一：
针对问题一，主要是从已知类别的人工制造序列中提取特征，构造分类方法，并将21~40的人工序号其分为A、B两类。
因为A、T、C、G这四个字符组成64种不同的3字符串，而这64种密码子构成20种氨基酸（包括61种密码子和3种终止密码），所以可以将DNA序列转化为氨基酸序列。
然而，在对DNA序列根据遗传密码翻译为氨基酸序列时，由于起始点位置的不同，含有三种可能的氨基酸序列，因此我们提取出20种氨基酸出现的平均频率。
以每种氨基酸出现的平均频率为特征，利用谱系聚类方法对20个已知类别的DNA序列进行分类，并用这些已知类别检验分类是否合理。
问题二：
针对问题二，对于1~182号自然序列，由于单个数据序列过长，我们根据实际意义进行简化。
抗微生物多肽也称防御素，对细菌、病菌、真菌等病原微生物能产生较强的杀伤作用，它们的分子量小于4000，仅有30多个氨基酸。
因此我们选取每个序列的前120的数据作为一组，计算每组数据的20种氨基酸频率。
因为精氨酸频率与合成防御素密切相关，将氨基酸频率按精氨酸频率进行排序，筛选出前40组数据作为肽链中第一个氨基酸最有可能为防御素的DNA序列。
再对这40组序列进行类似第一问ward聚类。
三、模型的假设
1.假设DNA序列按照文本的从左向右的顺序翻译成氨基酸序列；
2.不考虑氨基酸排列顺序对DNA序列所表示的生物信息的影响；
3.假设精氨酸的含量越高，越容易翻译为防御素；
4.只考虑DNA的前120个序列是否可以翻译为防御素；
四、符号定义及说明
符号
含义
xi
密码子平均频率
N
一条DNA序列中提取出的密码子总数
Ni
一条DNA序列中提取的第i种氨基酸出的密码子总数
Xi
第i条DNA序列的特征向量
Xi（k）
第i个特征在第k列序列中出现的频率
T
总离差平方和
欧式长度
Wk
类内离差平方和
PG
各类的类内离差平方和的总和
D（k）
距离矩阵，表示两类之间的距离
Gi（0）
第i个样本的初始状态的值
五、模型的建立与求解
5.1问题1
5.1.1基于氨基酸含量的特征提取法
首先将DNA序列按氨基酸的密码子翻译成氨基酸序列。
由生物学的知识，应用碱基互补配对的原则，将DNA序列翻译为RNA序列，再根据RNA中的64个密码子对应的氨基酸信息进行分类。
64种密码子构成了20种氨基酸，和编号为21的终止信息码。
各密码子所对应的氨基酸如表1所示。
表5-1-1每种氨基酸所对应的密码子
由于给出的DNA人工序列只是一个长长的DNA序列中的一个片段，不能够确定某个氨基酸是从哪个碱基开始，所以起点未知。
因此每条DNA序列可以译成三条可能的氨基酸序列。
为了说明方便，我们定义密码子平均频率xi：
i种密码子在三种情况下出现的总次数占取出的密码子总数的比例，即
（1）
对于一个给定的DNA序列，例如：
tgacctc，统计密码子个数时从第一个字母开始读起，对应密码子为tga，接着从第二个字母g读起，对应密码子为gac，依次类推，依次得到密码子（tga）（gac）（acc）（cct）（ctc）。
不难看出，此法提取出的密码子总数N为DNA序列长度D减去2，（程序见附录1）即
（2）
由此可以得到每条给定序列中20种氨基酸及终止信息码的含量，使之与一个21维的向量相对应，其中向量的分量表示各类氨基酸中各密码子的频率之和。
氨基酸的出现频率，终止信息码出现的频率如表5-1所示。
表5-1-221个特征值x1,x2,…,x21
氨基酸
苯丙氨酸
亮氨酸
异亮氨酸
甲硫氨酸
缬氨酸
特征
X1
X2
X3
X4
X5
氨基酸
丝氨酸
脯氨酸
苏氨酸
丙氨酸
组氨酸
特征
X6
X7
X8
X9
X10
氨基酸
酪氨酸
半胱氨酸
色氨酸
精氨酸
甘氨酸
特征
X11
X12
X13
X14
X15
氨基酸
天门冬氨酸
谷氨酸
赖氨酸
谷氨酰胺
天门冬酰胺
特征
X16
X17
X18
X19
X20
氨基酸
终止
特征
X21
记Xi={x1（k）,x2（k）,x3（k）,…,x21（k）}，其中xi（K）为第i个特征在第k列序列中出现的频率。
由MATLAB计算得到20种氨基酸及终止码出现的频率，统计结果如表所示（MATLAB程序见附录1-1，频率统计结果见附录1-2）
特征
氨基酸号
X1
X2
X3
X4
X5
1
0
0.072072072
0.027027027
0.009009009
0.027027027
2
0.027027
0.045045045
0.027027027
0.009009009
0.045045045
3
0.009009
0.009009009
0.036036036
0
0.009009009
4
0.054054
0.072072072
0.081081081
0.036036036
0.045045045
5
0.027027
0.054054054
0.018018018
0
0.009009009
6
0.045045
0.009009009
0.054054054
0.009009009
0.009009009
7
0.027027
0.045045045
0.063063063
0.045045045
0.009009009
8
0.036036
0.072072072
0.045045045
0.018018018
0.018018018
9
0.027027
0.045045045
0.018018018
0.027027027
0.036036036
10
0.027027
0.045045045
0.018018018
0
0.036036036
11
0.171171
0.126126126
0.126126126
0.018018018
0.045045045
12
0.162162
0.153153153
0.099099099
0.009009009
0.072072072
13
0.216216
0.153153153
0.081081081
0.018018018
0.045045045
14
0.18018
0.144144144
0.117117117
0.009009009
0.045045045
15
0.387387
0.09009009
0.117117117
0.009009009
0.036036036
16
0.117117
0.162162162
0.135135135
0.009009009
0.027027027
17
0.027027
0.099099099
0.072072072
0.036036036
0.027027027
18
0.162162
0.171171171
0.117117117
0
0.036036036
19
0.207207
0.189189189
0.135135135
0
0.027027027
20
0.243243
0.216216216
0.054054054
0
0.036036036
5.1.2基于谱系聚类法确定代表其信息的氨基酸序列
谱系聚类是一种逐次合并类的方法，最后得到一个聚类的二叉树聚类图。
其想法是，对于n个样本，先计算两两距离得到一个距离矩阵，然后将找到距离最近的两个样本合并为一个类，于是剩下n-1类（每个单独的未合并的样本作为一个类）。
计算n-1个类两两间的距离，找到距离最近的两个类将其合并，只剩下n-2个类……直至剩下两类。
类间距离的计算方法选择采用ward最小方差法
设观测的个数为n，变量的个数为v，G为在某一聚类水平上的个数，xi为第i个观测，Ck是当前（水平G）的第k类，Nk为Ck中的观测个数，
为均值向量，
为类Ck中的均值向量（中心），
为欧氏长度，则有
总离差平方和
（3）
类Ck的类内离差平方和
（4）
聚类水平G对应的各类的类内离差平方和的总和
（5）
假设某一步聚类把类Ck和类CL合并为下一水平的类CM，则定义
BK=WM-WK-WL为合并导致的类内离差平方和的增量。
DK表示类Gx和类Gy之间的距离。
（6）
下面我们将样本按距离准则逐步分类，类别由多到少，直至获得合适的分类要求位置，算法如下：
步骤一，将20个样本（20条人造DNA序列号）各自成一类，即建立20类。
记为G1（0），G2（0），…Gn（0）。
步骤二，计算20个样品两两间距离，得到一个20*20维的距离矩阵，记作D0。
这里的标号（0）表示聚类开始运行前的状态。
步骤三，合并距离最近的两类为一新类。
若前一步聚类运算中所得到的距离矩阵是D（k）（k为逐次聚类合并的次数），求出距离矩阵D（k）中的最小元素。
如果它是Gx（k）和Gy（k）两类之间的距离，则将Gx（k）和Gy（k）两类合并为一类Gxy（k+1），即Gxy（k+1）=Gx（k）
Gy（k）并由此建立新的分类：
G1（k+1），G2（k+1），…，G20（k+1）。
步骤四，计算出各新类与当前各类的距离，记为D（k+1）。
重复步骤三，继续计算及合并，直到达到所需要的聚类数目。
借助SPSS聚类得到的案例处理摘要如表5-1-3所示
案例处理摘要a
案例
有效
缺失
合计
N
百分比
N
百分比
N
百分比
20
100.0%
0
0.0%
20
100.0%
a.平方Euclidean距离已使用
图1
借助SPSS聚类得到的二叉树聚类如图1所示：
图1
得到的冰挂图如图2所示：
图2
由此得到的前20条DNA的分类如下
A：
1,2,3,5,6,7,8,9,10；
B：
4,11,12,13,14,15,16,17,18,19,20
与给出序列类型（1~10为A类，11~20为B类）相比，只有第四条DNA不符合已知序列。
样本回代正确率达95%，是较好的分类方法。
5.1.3用上述的谱系聚类法对编号为21~40继续分类：
首先统计出每条DNA的20种氨基酸及终止码的平均频率（MATLAB程序及统计结果见附录1-3）
对这20组序列进行氨基酸频率聚类，使用系统聚类中ward聚类方法，度量标准中区间选为平方Euclidean距离，标准为Z值。
可得聚类结果如下，表5-1-4为编号21~40的案例处理摘要，表5-1-5为编号21~40的查类表,图3为编号21~40的冰挂图，图4为编号21~40的二叉树聚类图：
案例处理摘要a
案例
有效
缺失
合计
N
百分比
N
百分比
N
百分比
20
100.0%
0
0.0%
20
100.0%
a.平方Euclidean距离已使用
表5-1-4
聚类表
阶
群集组合
系数
首次出现阶群集
下一阶
群集1
群集2
群集1
群集2
1
7
15
4.838
0
0
6
2
18
20
9.747
0
0
10
3
4
11
15.337
0
0
12
4
12
16
22.747
0
0
11
5
3
9
30.264
0
0
9
6
7
14
38.461
1
0
9
7
1
6
46.808
0
0
12
8
5
10
55.808
0
0
13
9
3
7
68.155
5
6
15
10
8
18
82.121
0
2
14
11
12
19
98.917
4
0
13
12
1
4
116.679
7
3
16
13
5
12
134.827
8
11
17
14
8
13
154.640
10
0
16
15
2
3
176.775
0
9
18
16
1
8
201.175
12
14
19
17
5
17
237.047
13
0
18
18
2
5
294.519
15
17
19
19
1
2
380.000
16
18
0
表5-1-5
图3
图4
5.2问题2
先使用MATLAB读入1~182序列中每组序列的前120个数据（MATLAB程序见附录2-1），统计每组序列的氨基酸含量（见附录2-2）。
找出精氨酸含量前40的序列见下表5-2-1：
DNA序列号
精氨酸频率
聚类代号
第16
0.237704918
1
第155
0.219512195
2
第171
0.219512195
3
第5
0.181818182
4
第90
0.180327869
5
第27
0.172131148
6
第111
0.170731707
7
第112
0.170731707
8
第49
0.163934426
9
第13
0.155737705
10
第69
0.155737705
11
第79
0.155737705
12
第14
0.147540984
13
第104
0.146341463
14
第109
0.146341463
15
第141
0.146341463
16
第159
0.146341463
17
第55
0.139344262
18
第61
0.139344262
19
第64
0.139344262
20
第15
0.131147541
21
第32
0.131147541
22
第106
0.130081301
23
第135
0.130081301
24
第158
0.130081301
25
第11
0.114754098
26
第70
0.114754098
27
第100
0.113821138
28
第117
0.113821138
29
第169
0.113821138
30
第28
0.106557377
31
第71
0.106557377
32
第88
0.106557377
33
第132
0.105691057
34
第147
0.105691057
35
第176
0.105691057
36
第179
0.105691057
37
第42
0.098360656
38
第47
0.098360656
39
第1
0.008264463
40
表5-2-1
对这40组序列进行氨基酸频率聚类，使用系统聚类中ward聚类方法，度量标准中区间选为平方Euclidean距离，标准为Z值。
可得聚类的案例处理摘要如下表5-2-2所示，查类表5-2-3，聚类冰挂图如图5，聚类树状图如图6：
案例处理摘要a
案例
有效
缺失
合计
N
百分比
N
百分比
N
百分比
40
100.0%
0
0.0%
40
100.0%
a.平方Euclidean距离已使用
表5-2-2
聚类表
阶
群集组合
系数
首次出现阶群集
下一阶
群集1
群集2
群集1
群集2
1
3
4
.315
0
0
17
2
7
9
3.419
0
0
26
3
17
23
7.649
0
0
20
4
6
8
11.944
0
0
31
5
14
26
17.856
0
0
30
6
11
21
23.938
0
0
27
7
18
19
30.268
0
0
8
8
15
18
37.451
0
7
25
9
12
16
45.038
0
0
19
10
10
13
52.658
0
0
20
11
25
38
60.629
0
0
24
12
30
32
68.662
0
0
13
13
30
39
77.560
12
0
22
14
22
35
86.619
0
0
25
15
24
27
95.921
0
0
21
16
28
31
105.873
0
0
19
17
2
3
116.204
0
1
26
18
33
36
126.796
0
0
23
19
12
28
138.670
9
16
31
20
10
17
150.876
10
3
28
21
20
24
165.419
0
15
29
22
29
30
180.309
0
13
30
23
33
40
195.401
18
0
32
24
25
37
210.554
11
0
29
25
15
22
226.075
8
14
28
26
2
7
242.888
17
2
34
27
5
11
262.541
0
6
34
升级会员 