DNA序列分类

DNA序列分类

实验目的

学习利用MATLAB提取DNA序列特征建立向量的方法，掌握利用FCM命令进行DNA分类的方法，学会做出分类图形直接给出分类结果的MATLAB编程。

知识扩展

DNA序列分类

DNA(Deoxyribonucleic acid)，中文译名为脱氧核苷酸，是染色体的主要化学成分，同时也是基因组成的，有时被称为“遗传微粒”。DNA是一种分子，可组成遗传指令，以引导生物发育与生命机能运作。主要功能是长期性的资讯储存，可比喻为“蓝图”或“食谱”。DNA分子是由两条核苷酸链以互补配对原则所构成的双螺旋结构的分子化合物。其中两条DNA链中对应的碱基A-T以双键形式连接，C-G以三键形式连接，糖-磷酸-糖 形成的主链在螺旋外侧，配对碱基在螺旋内侧。

FCM算法中样本点隶属于某一类的程度是用隶属度来反映的，不同的样本点以不同的隶属度属于每一类；但是算法中的概率约束∑uij=1使得样本的典型性反映

不出来，不适用于有噪音，样本分布不均衡，存在两个或者两个以上样本分别距两个类的距离相等的样本等等。

欧氏距离（ Euclidean distance）也称欧几里得距离，它是一个通常采用的距离定义，它是在m维空间中两个点之间的真实距离。 公式

在二维和三维空间中的欧式距离的就是两点之间的距离，二维的公式是 d = sqrt((x1-x2)^+(y1-y2)^) 三维的公式是

d=sqrt(x1-x2)^+(y1-y2)^+(z1-z2)^) 推广到n维空间，欧式距离的公式是 d=sqrt( ∑(xi1-xi2)^ ) 这里i=1,2..n

xi1表示第一个点的第i维坐标,xi2表示第二个点的第i维坐标

n维欧氏空间是一个点集,它的每个点可以表示为(x(1),x(2),...x(n)),其中x(i)(i=1,2...n)是实数,称为x的第i个坐标,两个点x和y=(y(1),y(2)...y(n))之间的距离d(x,y)定义为上面的公式. 欧氏距离判别准则如下：

若dAdB,则将Xi点判为B类

若dA=dB,则将Xi点判为不可判别点。 欧氏距离看作信号的相似程度。 距离越近

1.问题的提出

2000年6月，人类基因组计划中DNA全序列草图完成，预计2001以完精确的全序列图，此后人类将拥有一本记录着自身生老病死及遗的全部信息的“天书”，这本大自然写成的“天书”，是由4个字符A,T,C,G按一定顺序排成的长约30亿的序列，其中没有“断句”也没有标点符号，除了这4种碱基以外，人们对它包含的“内容”知之甚少，难以读懂，破译这部世界上最巨量信息的“天书”是21实际最重要的任务之一。在这个目标中，研究DNA全序列具有什么结构，由这4个字符排成的看似随机的序列中隐藏着什么规律，又是解读这部天书的基础，是生物信息学最重要的课题之一。

2．问题的分析

这是一个比较典型的分类问题，为了表述的严格和方便，我们用数学的方法来重述这个问题。在这里问题的关键就是要从已知的20个字母序列中提取用于分类的特征。知道了这些特征，我们就可以比较容易的，对那些未标明类型的序列进行分类，下面我们将首先对用于分类的标准问题进行必要的讨论。

3．分类的方法

为了在众多可能的分类中寻求合理的分类结果，为此，就要确定合理的聚类准则。定义目标函数为

J(U,V)(uik)(dik)k1i1202m2

显然，J(U,V)表示了各类中样本到聚类中心的加权距离平方和，权重是样本XK对第i类隶

度的m次方，聚类准则取为求J(U,V)的极小值(min){J(U,V)}。

其中，U=[uik]为模糊分类矩阵，i=1，2；k=1，2，···，20；且满足0≤ uik ≤1和

若uik=max{uik}>0.5，则xk∈第j类。

在MATLAB中，我们只要直接调用如下程序即可：

[Center,U,obj_fcm]=fcm(data,cluster_n) data:要聚类的数据函数，每一行为一个样本 cluster_n:聚类数（大于1）

Center:最终的聚类中心矩阵，其每一行为聚类中心的坐标值 U:最终的模糊分区矩阵

obj_fcm:在迭代过程中的目标函数值

4. 对DNA序列组合分类的分析

（1）提取DNA序列特征建立两类序列的特征向量 （2）确定两类序列的中心 （3）分类方法 （4）回代误判率

（5）未知的20个序列判别结果

5.提取DNA序列特征建立两类序列的特征向量

为了对DNA序列进行分类，我们首先对已知的两类DNA序列进行研究，从中找到两类序列的特征。由于在不用于编码蛋白质的序列片段中，A和T的含量特别多些，于是我们利用MATLAB软件，通过编程计算出A，B两类序列中4中碱基对含量的百分比，对每个序列构造四维向量 xk =( xk1，xk2，xk3，xk4 ) (k=1,2···,20)

其中， xk1，xk2，xk3，xk4分别表示第k个序列所含有的碱基对A,T,C,G含量的百分比，利用MATLAB软件，我们得到A,B两类序列的特征矩阵A=(xkj)20×4

6. 确定两类序列的中心

v(uik120mik)x/(ukk1220mik)(i1,2;1m)dik1/(uik)j12m1

为第k个序列到第i类中心的欧式距离，实际计算时要对取定的初始值进行迭代

tt-1

计算直至max｛|uik-uik|｝<ɛ ，ɛ为事先指定的精度。 回代误判率

用欧氏距离作为判据虽然简洁直观，但是存在着明显的缺陷：从概率统计的角度来看，用欧氏距离描述随机点之间的距离并不是很好。因此对待分类样本是随机样本，具有一定的统计性质时，这个模型并不是很能很好的描述两个随机点之间的接近程度。

我们对于已知的A，B两类序列利用上述方宣传法进行判别，结果如3.15表所示

属于A类的隶属值序号12345678910111213141516171819200.99690.99570.96770.43560.97730.96230.91680.98420.97380.94200.02310.03180.01150.00220.06700.03190.49850.00450.04330.06550.00310.00430.03230.56410.02270.03770.08320.01580.02620.05800.97690.96820.98850.99780.93300.96810.50150.99550.95670.9345AAABAAAAAABBBBBBBBBB属于B类的隶属值聚类结果djkdikxkvi

未知的20个序列判别结果 利用模糊均值聚类方法，对于未知的20个序列进

行判别，结果如表3.16

序号2122232425262728293031323334353637383940属于A类的隶属值0.04790.79720.92300.02010.97860.01240.95630.07410.83360.58080.05440.39870.11370.95760.97860.84920.99190.04280.64440.0488属于B类的隶属值0.95210.20280.07700.97990.02140.98760.04370.92590.16640.41920.94560.60130.88630.04240.02140.15080.00810.95720.35560.9512聚类结果BAABABABAABBBAAAABAB

MATLAB编程

(1) 提取特征建立特征矩阵

A1='aggcacggaaaaacgggaataacggaggaggacttggcacggcattacacggaggacgaggtaaaggaggcttgtctacggccggaagtgaagggggatatgaccgcttgg';

A2='cggaggacaaacgggatggcggtattggaggtggcggactgttcggggaattattcggtttaaacgggacaaggaaggcggctggaacaaccggacggtggcagcaaagga';

A3='gggacggatacggattctggccacggacggaaaggaggacacggcggacatacacggcggcaacggacggaacggaggaaggagggcggcaatcggtacggaggcggcgga';

A4='atggataacggaaacaaaccagacaaacttcggtagaaatacagaagcttagatgcatatgttttttaaataaaatttgtattattatggtatcataaaaaaaggttgcga';

A5='cggctggcggacaacggactggcggattccaaaaacggaggaggcggacggaggctacaccaccgtttcggcggaaaggcggagggctggcaggaggctcattacggggag';

A6='atggaaaattttcggaaaggcggcaggcaggaggcaaaggcggaaaggaaggaaacggcggatatttcggaagtggatattaggagggcggaataaaggaacggcggcaca';

A7='atgggattattgaatggcggaggaagatccggaataaaatatggcggaaagaacttgttttcggaaatggaaaaaggactaggaatcggcggcaggaaggatatggaggcg';

A8='atggccgatcggcttaggctggaaggaacaaataggcggaattaaggaaggcgttctcgcttttcgacaaggaggcggaccataggaggcggattaggaacggttatgagg';

A9='atggcggaaaaaggaaatgtttggcatcggcgggctccggcaactggaggttcggccatggaggcgaaaatcgtgggcggcggcagcgctggccggagtttgaggagcgcg';

A10='tggccgcggaggggcccgtcgggcgcggatttctacaagggcttcctgttaaggaggtggcatccaggcgtcgcacgctcggcgcggcaggaggcacgcgggaaaaaacg';

A11='gttagatttaacgttttttatggaatttatggaattataaatttaaaaatttatattttttaggtaagtaatccaacgtttttattactttttaaaattaaatatttatt';

A12='gtttaattactttatcatttaatttaggttttaattttaaatttaatttaggtaagatgaatttggttttttttaaggtagttatttaattatcgttaaggaaagttaaa';

A13='gtattacaggcagaccttatttaggttattattattatttggattttttttttttttttttttaagttaaccgaattattttctttaaagacgttacttaatgtcaatgc';

A14='gttagtcttttttagattaaattattagattatgcagtttttttacataagaaaatttttttttcggagttcatattctaatctgtctttattaaatcttagagatatta';

A15='gtattatatttttttatttttattattttagaatataatttgaggtatgtgtttaaaaaaaatttttttttttttttttttttttttttttttaaaatttataaatttaa';

A16='gttatttttaaatttaattttaattttaaaatacaaaatttttactttctaaaattggtctctggatcgataatgtaaacttattgaatctatagaattacattattgat';

A17='gtatgtctatttcacggaagaatgcaccactatatgatttgaaattatctatggctaaaaaccctcagtaaaatcaatccctaaacc~ttaaaaaacgg~ggcctat~cc';

A18='gttaattatttattccttacgggcaattaattatttattacggttttatttacaattttttttttttgtcctatagagaaattacttacaaaa~gttattttacatactt';

A19='gttacattatttattattatccgttatcgataattttttacctcttttttcgctgagtttttattcttactttttttcttctttatataggatctcatttaatatcttaa';

A20='gtatttaactctctttactttttttttcactctctacattttcatcttctaaaactgtttgatttaaacttttgtttctttaaggattttttttacttatc~t~tgttat';

A21='tttagctcagtccagctagctagtttacaatttcgacaccagtttcgcaccatcttaaatttcgatccgtaccgtaatttagcttagatttggatttaaaggatttagattga';

A22='tttagtacagtagctcagtccaagaacgatgtttaccgtaacgtacgtaccgtacgctaccgttaccggattccggaaagccgattaaggaccgatcgaaaggg';

A23='cgggcggatttaggccgacggggacccgggattcgggacccgaggaaattcccggattaaggtttagcttcccgggatttagggcccggatggctggga~cc';

A24='tttagctagctactttagctatttttagtagctagccagcctttaaggctagctttagctagcattgttctttattgggacccaagttcgacttttacgattta~ttttgaccgt';

A25='gaccaaaggtgggctttagggacccgatgctttagtcgcagctggaccagttccccagggtattaggcaaaagctgacgggcaattgcaatttaggcttaggcca';

A26='gatttactttagcatttttagctgacgttagcaagcattagctttag~caatttcgcatttgccagtttcgcagctcagttttaacgcgggatctttagcttcaagctttttac';

A27='ggattcggatttacccggggattggcggaacgggacctttaggtcgggacccattaggagtaaatgccaaaggacgctggtttagccagtccgttaaggcttag';

A28='tccttagatttcagttactatatttgacttacagtctttgagatttcccttacgattttgacttaaaatttagacgttagggcttatcagttatggattaatttagcttattttcga';

A29='ggccaattccggtaggaaggtgatggcccgggggttcccgggaggatttaggctgacgggccggccatttcggtttagggagggccgggacgcgttagggc';

A30='cgctaagcagctcaagctcagtcagtcacgtttgccaagtcagtaatttgccaaagttaaccgttagctgacgctgaacgctaaacagtattagctgatgactcgta';

A31='ttaaggacttaggctttagcagttactttagtttagttccaagctacgtttacgggaccagatgctagctagcaatttattatccgtattaggcttaccgtaggttta~cgt';

A32='gctaccgggcagtctttaacgtagctaccgtttagtttgggcccagccttgcggtgtttcggattaaattcgttgtcagtcgctcttgggtttagtcattcccaaaagg';

A33='cagttagctgaatcgtttagccatttgacgtaaacatgattttacgtacgtaaattttagccctgacgtttagctaggaatttatgctgacgtagcgatcgactttagcac';

A34='cggttagggcaaaggttggatttcgacccagggggaaagcccgggacccgaacccagggctttagcgtaggctgacgctaggcttaggttggaacccggaaa';

A35='gcggaagggcgtaggtttgggatgcttagccgtaggctagctttcgacacgatcgattcgcacacaggataaaagttaagggaccggtaagtcgcggtagcc';

A36='ctagctacgaacgctttaggcgcccccgggagtagtcgttaccgttagtatagcagtcgcagtcgcaattcgcaaaagtccccagctttagccccagagtcgacg'; A37='gggatgctgacgctggttagctttaggcttagcgtagctt

tagggccccagtctgcaggaaatgcccaaaggaggccaccgggtagatgccasagtgcaccgt'; A38='aacttttagggcatttccagttttacgggttattttcccagttaaactttgcaccattttacgtgttacgatttacgtataatttgaccttattttggacactttagtttgggttac';

A39='ttagggccaagtccccaggcaaggaattctgatccaagtccaatcacctacagtccaagtcaccgtttgcagctaccgtttaccgtacgttgcaagtcaaatccat';

A40='ccattclgggtttatttacctctttattttttcccgagaccttaggtttaccgtactttttaacggtttacctttgaaatttttggactagcttaccctggatttaacggccagttt'; s=['A'];

x1=zeros(40,4); for i=1:40

n=eval([s(1,1),num2str(i)]); [u,y]=size(n); for j=1:y; m=n(:,j); switch m case 'a';

x1(i,1)=x1(i,1)+1; case 't';

x1(i,2)=x1(i,2)+1; case 'c';

x1(i,3)=x1(i,3)+1; case 'g';

x1(i,4)=x1(i,4)+1; end end end

b=sum(x1'); b1=b';

c=[b1,b1,b1,b1]; y=x1./c

运行结果： y =

0.2973 0.1351 0.1712 0.3964 0.2703 0.1532 0.1622 0.4144

0.2703 0.0631 0.2162 0.4505 0.4234 0.2883 0.1081 0.1802 0.2342 0.1081 0.2342 0.4234 0.3514 0.1261 0.1261 0.3964 0.3514 0.1892 0.0991 0.3604 0.2793 0.1892 0.1622 0.3694 0.2072 0.1532 0.2072 0.4324 0.1818 0.1364 0.2727 0.4091 0.3545 0.5000 0.0455 0.1000 0.3273 0.5000 0.0273 0.1455 0.2545 0.5182 0.1000 0.1273 0.3000 0.5000 0.0818 0.1182 0.2909 0.6455 0 0.0636 0.3636 0.4636 0.0818 0.0909 0.3645 0.2710 0.2243 0.1402 0.2936 0.5046 0.1101 0.0917 0.2182 0.5636 0.1455 0.0727 0.2037 0.5741 0.1574 0.0648 0.2743 0.3628 0.1947 0.1681 0.2885 0.2212 0.2404 0.2500 0.1782 0.1881 0.2475 0.3861 0.2105 0.4123 0.1930 0.1842 0.2476 0.2190 0.2286 0.3048 0.2212 0.3894 0.2035 0.1858 0.2308 0.2308 0.2019 0.3365 0.2564 0.4444 0.1453 0.1538 0.1485 0.1881 0.2178 0.4455 0.2897 0.2523 0.2430 0.2150 0.2432 0.3604 0.1802 0.2162 0.1743 0.3303 0.2294 0.2661 0.2703 0.3333 0.1892 0.2072 0.2353 0.1667 0.2353 0.3627 0.2451 0.2059 0.2059 0.3431 0.2286 0.2095 0.3048 0.2571 0.2157 0.2059 0.2451 0.3333 0.2222 0.4359 0.1709 0.1709 0.2736 0.2358 0.3019 0.1887 0.1897 0.4310 0.2155 0.1638

(2) 回代判别 Y=y(1:20,:);

[center,U,obj_fcn]=fcm(Y,2); maxU=max(U);

index1=find(U(1,:)==maxU); index2=find(U(2,:)==maxU); line(Y(index1,1),

Y(index1,2),'linestyle','none','marker','o','color','g');

line(Y(index2,1),Y(index2,2),'linestyle','none','marker','*','color','r')

运行结果：

小于0.5的为A类,大于0.5的为B类.

(3) 未知的20个序列判别 Y=y(21:40,:);

[center,U1,obj_fcn]=fcm(Y,2); maxU1=max(U1);

index1=find(U1(1,:)==maxU1); index2=find(U1(2,:)==maxU1);

line(Y(index1,1),Y(index1,2),'linestyle','none','marker','o','color','g');

line(Y(index2,1),Y(index2,2),'linestyle','none','marker','*','color','r'); hold on

plot(center(1,1),center(1,2),'kpentagram','markersize', 15,'LineWidth',2);

plot(center(2,1),center(2,2),'ksquare','markersize', 15,'LineWidth',2);

box on,title('模糊C均值聚类');

r11=0;r12=0

t=U1(1,:)-U1(2,:); for i=1:20, if t(i)<0

r11=r11+1; t1(i)=i; else if t(i)>0 r12=r12+1; t2(i)=i; end end end

sprintf('属于A类的DNA序列序号') t2

sprintf('属于B类的DNA序列序号') t1

运行结果：

属于A类的DNA序列序号 t2 =

Columns 1 through 15

0 0 3 0 5 0 7 0 9 0 0 0 0 14 15

Columns 16 through 19 0 17 0 19 ans =

属于B类的DNA序列序号 t1 =

Columns 1 through 15

1 2 0 4 0 6 0 8 0 10 11 12 13 0 0

Columns 16 through 20

16 0 18 0 20

未知的20个序列聚类图（ × 表示A类， 表示B类， 示两类中心）

表，