装订线装订线（威海）《数学软件实践》课程设计报告题目：MATLAB实现DNA序列的分类识别学号：姓名：教师日期：论文题目论文题目 (2)摘要 (3)1 问题重述 (4)2 问题分析 (4)3 模型假设 (4)4 符号说明 (5)5 模型的建立与求解 (5)5.1 问题一的模型 (5)5.1.1 模型建立 (5)模型I (6)5.1.2 模型求解 (6)5.2 问题二的模型 (7)5.2.1 模型建立 (7)模型II (8)5.2.2 模型求解 (9)5.3 问题三的模型 (10)5.3.1 模型建立 (10)模型III (10)5.3.2 模型求解 (11)6 模型的评价与改进 (12)附录 (13)A.1 附录描述1 (13)A.2 附录描述2 (20)A.3 附录描述3 (23)摘 要1) DNA 序列矩阵a) 首先调用函数_I importdata ，将txt 文档中的字符读入到workspace 中，将其中不是ATGC 的字符用A 来替换，存到矩阵DNA 中。

b) 给定DNA 序列12N S s s s =， {}i s A C G T ∈,,,，,,,A T G C 分别表示序列S 的前i (12i N =,,...,)个元素中A C G T ,,,出现的次数.编写函数_I calculate 计算出DNA 序列每一行的,,,A T G C 。

c) 根据定义的欧几里得距离矩阵，为了求解ED 方便,写编写一个函数_1I ED ，求出四个矩阵,,,Aij Tij Gij Cij ，然后编写函数_I ED ，便可以求出AG CT AC GT AT GC ，，，，，对应的ED 。

d) 根据定义的路径距离矩阵，编写函数_I PD ，根据不同的ED ，求出不同的PD 。

e) 根据定义的商矩阵，编写函数_I EP ，由(b)和(c)求出的ED 和PD 计算出对应的EP 。

f) 根据定义的商矩阵，编写函数_I EG ，由(b)求出的ED 计算出对应的EG 。

2) 最大特征矩阵对于给定的DNA 序列， 令123456()μμμμμμμ=,,,,,,其中123456μμμμμμ,,,,,分别表示矩阵,,ED EP EG 的标准化(i N μλ=/)的最大特征根。

将求得的结果放于矩阵123EGa EGa EGa ，，中。

3) 调用给定程序作出DNA 序列的聚类树图可以根据给定的程序_I dendrogram ，对于123EGa EGa EGa ，，分别作图。

4) MATLAB 的一些基础知识关键字：DNA 序列、最大特征根、聚类树图、MATLAB 基础知识1 问题重述题目给的是关于DNA 序列的问题,首先给出一个txt 文档，是构成DNA 序列的ATGC 字符，但其中有不是ATGC 的字符，要求将txt 文档中的所有字符序列读入到矩阵DNA 中，并将其中不是ATGC 的字符用ATGC 中的任意一个换掉，得到正确的DNA 矩阵。

然后分别计算矩阵,,,A T G C (例如i T 表示前i 个字符中含有的T 的个数)。

再根据,,,A T G C 四个矩阵不同的输入顺序求出6个不同的ED 矩阵，对每个ED 矩阵，可以求出对应的PD ,再由,ED PD 求出对应的EP 。

对每个ED 矩阵，也可以求出对应的EG 矩阵。

假设输入DNA 序列的第一行，求出6个不同的,,ED EP EG 后，每个矩阵可以求得一个最大特征根，这样6个ED 有6个最大特征根，放于矩阵1EGa (1EGa 为24*6的矩阵)的第一行中，6个EP 有6个最大特征根，放于矩阵2EGa (2EGa 为24*6的矩阵)的第一行中，6个EG 有6个最大特征根，放于矩阵3EGa (3EGa 为24*6的矩阵)的第一行中。

当把DNA 序列的24行遍历完，1,2,3EGa EGa EGa 三个矩阵便赋值完成。

提示中给了画DNA 序列的聚类树图的程序，调用此程序，便可完成3个图的绘制。

2 问题分析首先调用函数_I importdata ，将txt 文档中的字符读入到矩阵DNA 中，将其中不是ATGC 的字符用A 来替换。

编写函数_I calculate 计算出DNA 序列每一行的,,,A T G C 。

根据定义的欧几里得距离矩阵,为了求解ED 方便,先由函数_1I ED ，求出四个矩阵,,,Aij Tij Gij Cij ，然后由函数_I ED ，便可以求出,,,,,AG CT AC GT AT GC 对应的ED 。

根据定义的路径距离矩阵，由函数_I PD ，根据不同的ED ，求出不同的PD 。

根据定义的商矩阵，由函数_I EP 和求出的ED 和PD 计算出对应的EP 。

根据定义的商矩阵，由函数_I EG 求出的ED 计算出对应的EG 。

编写函数_1I EGa ，假设输入DNA 序列的第一行，求出6个不同的,,ED EP EG 后，每个矩阵对应一个最大特征根，6个ED 的6个最大特征根放于矩阵1EGa 的第一行中，6个EP 的6个最大特征根放于矩阵2EGa 的第一行中，6个EG 的6个最大特征根放于矩阵3EGa 的第一行中。

遍历完DNA 序列的24行，得到1,2,3EGa EGa EGa 三个最大特征根矩阵。

由提示中给定的画DNA 序列聚类树图的程序画出对应的聚类树图。

3 模型假设1) 假设实验数据足够准确。

2) 假设实验过程中不存在仪器误差。

3) 假设提示中给的程序完整且正确。

4 符号说明表 1矩阵，存放从文档中取出的数据，并将非正常字符转换。

A 向量，其中()A i 表示序列前i 个中含有的A 的个数 T 向量，其中()T i 表示序列前i 个中含有的T 的个数 G 向量，其中()G i 表示序列前i 个中含有的G 的个数 C向量，其中()C i 表示序列前i 个中含有的C 的个数A ij 矩阵，ij j i A A A =- ij T 矩阵，ij j i T T T =- ij G 矩阵，ij j i G G G =- ij C 矩阵，ij j i C C C =-ED 矩阵，根据ACGT 输入顺序的不同，产生不同的ED PD 矩阵，1121[][][][][]ji ij i i i i j j PD PD ED ED ED ,++,+-,==+++ EP 矩阵，[][][][][]0ij ij ij ij ii EP E P ED PD i j EP =/=/,≠;=. EG 矩阵，[][][][]0ij ij ij ii EG E G ED i j i j EG =/=/|-|,≠;=. 1EGa 矩阵，存放ED 的最大特征根 2EGa矩阵，存放EP 的最大特征根 3EGa 矩阵，存放EG 的最大特征根5 模型的建立与求解5.1 问题一的模型读取数据到矩阵DNA 中，并将其中不是ATGC 的字符用A 替换掉，得到处理过的DNA 矩阵。

并根据输入的DNA 序列的每一行，计算出4个矩阵,,,A T G C 。

5.1.1 模型建立首先，txt 文档中的数据未经过处理，必须编写一个函数_I importdata ，将txt 中的所有字符读入到工作空间中，其中共有24个序列，每一个序列都以>‘’开始，此时函数中要有判断语句，如果遇到>‘’就将它的下一行序列读入到矩阵DNA 中，读入完成以后，要运用一段程序，判断其中的字符若不是ATGC ，就用A 替换。

这样就很清楚的得到了DNA 序列矩阵。

提示中给定的i i i A C G ,,,i T 分别表示序列S 的前i (12i N =,,...,)个元素中A C G T ,,,出现的次数。

编写函数_I calculate ，根据DNA 中某一行的输入，得出需要的ATGC 矩阵。

模型I1. 读取部分：P=importdata('TF24.txt'); %读取TF 中的数据到a 中，P 是cell 类型的列向量2. 判断部分： for i=1:mif double(Q{i}(1))==62A(j)=i; %A 数组存放的的是首字符为'>'的行数 j=j+1; endend3. 取数据部分：for i=1:nDNA{i}=cell2mat((Q(A(i)+1:A(i+1)-1))');%将介于两个'>'之间的段取出，然后转置，然后分开成单个字符，放于cell 的第一个中end4. 替换非正常数据部分：if DNA{i}(j)~='A'&&DNA{i}(j)~='T'&&DNA{i}(j)~='G'&&DNA{i}(j)~='C'disp('此处有一个字符不是ATGC');DNA{i}(j)='A';disp('转换完成！');end5.1.2 模型求解 部分结果展示图表 1DNA 矩阵图表 2 输入DNA 的第一行后得出的矩阵A 的部分结果5.2 问题二的模型要求解,,ED EP EG 的最大特征根矩阵，首先要把ED 求解出来，而在求解ED 时，涉及到AG CT AC GT AT GC ，，，，，这6个，而在求解它们是还有,,,Aij Tij Gij Cij 这4个矩阵，AG CT AC GT AT GC ，，，，，的不同只是因为,,,Aij Tij Gij Cij 输入顺序的不同而产生的。

所以先设计函数把,,,Aij Tij Gij Cij 求解出来。

再由函数求解不同的ED ，根据求解出来的ED ，解相应的,,PD EP EG 就变得简单多了。

取出其中的ED 及其对应的,EP EG ，分别求它们的最大特征根矩阵，放到三个数组1,2,3EGa EGa EGa 中，再根据给定的作图程序，做出DNA 序列的聚类树图。

5.2.1 模型建立1. 为了求解ED 方便,写编写一个函数_1I ED ，两层for 循环遍历求出,,,Aij Tij Gij Cij 。

2. 由函数_I ED ，其中输入4个形参，便可以求出AG CT AC GT AT GC ，，，，，对应的ED 。

由函数_I PD ，根据不同的ED ，求出不同的PD ，我们注意到，PD 是对称矩阵，为了减少运算量，只求取PD 的上三角矩阵便可。

3. 由函数_I EP ，由于已经求出ED 和PD ，便可计算出对应的EP 。

此时求解中注意到有除法运算，而且最后对角线上全是0。

为了避免除数为0的情况，将对角阵的值全部变为1，运算完成以后再将对角线赋值为0即可。

4. 由函数_I EG ，根据求出的ED 计算出对应的EG ，只需要一个if else -判断即可。