《DNA序列分类.doc》由会员分享，可在线阅读，更多相关《DNA序列分类.doc（13页珍藏版）》请在三一文库上搜索。
1、DNA序列分类 推荐精选实验目的学习利用MATLAB提取DNA序列特征建立向量的方法，掌握利用FCM命令进行DNA分类的方法，学会做出分类图形直接给出分类结果的MATLAB编程。知识扩展DNA序列分类DNA(Deoxyribonucleic acid)，中文译名为脱氧核苷酸，是染色体的主要化学成分，同时也是基因组成的，有时被称为“遗传微粒”。DNA是一种分子，可组成遗传指令，以引导生物发育与生命机能运作。主要功能是长期性的资讯储存，可比喻为“蓝图”或“食谱”。DNA分子是由两条核苷酸链以互补配对原则所构成的双螺旋结构的分子化合物。其中两条DNA链中对应的碱基A-T以双键形式连接，C-G以三键形
2、式连接，糖-磷酸-糖 形成的主链在螺旋外侧，配对碱基在螺旋内侧。FCM算法中样本点隶属于某一类的程度是用隶属度来反映的，不同的样本点以不同的隶属度属于每一类；但是算法中的概率约束uij=1使得样本的典型性反映不出来，不适用于有噪音，样本分布不均衡，存在两个或者两个以上样本分别距两个类的距离相等的样本等等。推荐精选欧氏距离（ Euclidean distance）也称欧几里得距离，它是一个通常采用的距离定义，它是在m维空间中两个点之间的真实距离。公式在二维和三维空间中的欧式距离的就是两点之间的距离，二维的公式是 d = sqrt(x1-x2)+(y1-y2) 三维的公式是 d=sqrt(x1-x
3、2)+(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)定义为上面的公式. 欧氏距离判别准则如下： 若dA<dB,则将Xi点判为A类 若dA>dB,则将Xi点判为B类 若dA=dB,则将Xi点判为不可判别点。欧氏距离看作信号的相似程度。 距离越近
4、1.问题的提出2000年6月，人类基因组计划中DNA全序列草图完成，预计2001以完精确的全序列图，此后人类将拥有一本记录着自身生老病死及遗的全部信息的“天书”，这本大自然写成的“天书”，是由4个字符A,T,C,G按一定顺序排成的长约30亿的序列，其中没有“断句”也没有标点符号，除了这4种碱基以外，人们对它包含的“内容”知之甚少，难以读懂，破译这部世界上最巨量信息的“天书”是21实际最重要的任务之一。在这个目标中，研究DNA全序列具有什么结构，由这4个字符排成的看似随机的序列中隐藏着什么规律，又是解读这部天书的基础，是生物信息学最重要的课题之一。2问题的分析这是一个比较典型的分类问题，为了表述
5、的严格和方便，我们用数学的方法来重述这个问题。在这里问题的关键就是要从已知的20个字母序列中提取用于分类的特征。知道了这些特征，我们就可以比较容易的，对那些未标明类型的序列进行分类，下面我们将首先对用于分类的标准问题进行必要的讨论。3分类的方法为了在众多可能的分类中寻求合理的分类结果，为此，就要确定合理的聚类准则。定义目标函数为推荐精选显然，J(U,V)表示了各类中样本到聚类中心的加权距离平方和，权重是样本XK对第i类隶 度的m次方，聚类准则取为求J(U,V)的极小值(min)J(U,V)。 其中，U=uik为模糊分类矩阵，i=1，2；k=1，2，···，20；且
6、满足0 uik 1和 若uik=maxuik>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序列特征建立两类序列的特征向
7、量为了对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×46. 确定两类序列的中心推荐精选为第k个序列到第i类中
8、心的欧式距离，实际计算时要对取定的初始值进行迭代计算直至max|uikt-uikt-1|< ，为事先指定的精度。回代误判率 用欧氏距离作为判据虽然简洁直观，但是存在着明显的缺陷：从概率统计的角度来看，用欧氏距离描述随机点之间的距离并不是很好。因此对待分类样本是随机样本，具有一定的统计性质时，这个模型并不是很能很好的描述两个随机点之间的接近程度。 我们对于已知的A，B两类序列利用上述方宣传法进行判别，结果如3.15表所示 未知的20个序列判别结果 利用模糊均值聚类方法，对于未知的20个序列进行判别，结果如表3.16推荐精选MATLAB编程(1) 提取特征建立特征矩阵A1='aggc
9、acggaaaaacgggaataacggaggaggacttggcacggcattacacggaggacgaggtaaaggaggcttgtctacggccggaagtgaagggggatatgaccgcttgg'A2='cggaggacaaacgggatggcggtattggaggtggcggactgttcggggaattattcggtttaaacgggacaaggaaggcggctggaacaaccggacggtggcagcaaagga'A3='gggacggatacggattctggccacggacggaaaggaggacacggcggacatacacg
10、gcggcaacggacggaacggaggaaggagggcggcaatcggtacggaggcggcgga'A4='atggataacggaaacaaaccagacaaacttcggtagaaatacagaagcttagatgcatatgttttttaaataaaatttgtattattatggtatcataaaaaaaggttgcga'A5='cggctggcggacaacggactggcggattccaaaaacggaggaggcggacggaggctacaccaccgtttcggcggaaaggcggagggctggcaggaggctcattacggg
11、gag'A6='atggaaaattttcggaaaggcggcaggcaggaggcaaaggcggaaaggaaggaaacggcggatatttcggaagtggatattaggagggcggaataaaggaacggcggcaca'A7='atgggattattgaatggcggaggaagatccggaataaaatatggcggaaagaacttgttttcggaaatggaaaaaggactaggaatcggcggcaggaaggatatggaggcg'A8='atggccgatcggcttaggctggaaggaacaaatagg
12、cggaattaaggaaggcgttctcgcttttcgacaaggaggcggaccataggaggcggattaggaacggttatgagg'A9='atggcggaaaaaggaaatgtttggcatcggcgggctccggcaactggaggttcggccatggaggcgaaaatcgtgggcggcggcagcgctggccggagtttgaggagcgcg'A10='tggccgcggaggggcccgtcgggcgcggatttctacaagggcttcctgttaaggaggtggcatccaggcgtcgcacgctcggcgcgg
13、caggaggcacgcgggaaaaaacg'推荐精选A11='gttagatttaacgttttttatggaatttatggaattataaatttaaaaatttatattttttaggtaagtaatccaacgtttttattactttttaaaattaaatatttatt'A12='gtttaattactttatcatttaatttaggttttaattttaaatttaatttaggtaagatgaatttggttttttttaaggtagttatttaattatcgttaaggaaagttaaa'A13='gtattacaggc
14、agaccttatttaggttattattattatttggattttttttttttttttttttaagttaaccgaattattttctttaaagacgttacttaatgtcaatgc'A14='gttagtcttttttagattaaattattagattatgcagtttttttacataagaaaatttttttttcggagttcatattctaatctgtctttattaaatcttagagatatta'A15='gtattatatttttttatttttattattttagaatataatttgaggtatgtgtttaaaaaaaat
15、ttttttttttttttttttttttttttttttaaaatttataaatttaa'A16='gttatttttaaatttaattttaattttaaaatacaaaatttttactttctaaaattggtctctggatcgataatgtaaacttattgaatctatagaattacattattgat'A17='gtatgtctatttcacggaagaatgcaccactatatgatttgaaattatctatggctaaaaaccctcagtaaaatcaatccctaaaccttaaaaaacggggcctatcc'A18
16、='gttaattatttattccttacgggcaattaattatttattacggttttatttacaattttttttttttgtcctatagagaaattacttacaaaagttattttacatactt'A19='gttacattatttattattatccgttatcgataattttttacctcttttttcgctgagtttttattcttactttttttcttctttatataggatctcatttaatatcttaa'A20='gtatttaactctctttactttttttttcactctctacattttcatct
17、tctaaaactgtttgatttaaacttttgtttctttaaggattttttttacttatcttgttat'A21='tttagctcagtccagctagctagtttacaatttcgacaccagtttcgcaccatcttaaatttcgatccgtaccgtaatttagcttagatttggatttaaaggatttagattga'A22='tttagtacagtagctcagtccaagaacgatgtttaccgtaacgtacgtaccgtacgctaccgttaccggattccggaaagccgattaaggaccgatcg
18、aaaggg'A23='cgggcggatttaggccgacggggacccgggattcgggacccgaggaaattcccggattaaggtttagcttcccgggatttagggcccggatggctgggacc'A24='tttagctagctactttagctatttttagtagctagccagcctttaaggctagctttagctagcattgttctttattgggacccaagttcgacttttacgatttattttgaccgt'A25='gaccaaaggtgggctttagggacccgatgctttagtc
19、gcagctggaccagttccccagggtattaggcaaaagctgacgggcaattgcaatttaggcttaggcca'A26='gatttactttagcatttttagctgacgttagcaagcattagctttagcaatttcgcatttgccagtttcgcagctcagttttaacgcgggatctttagcttcaagctttttac'A27='ggattcggatttacccggggattggcggaacgggacctttaggtcgggacccattaggagtaaatgccaaaggacgctggtttagccagtc
20、cgttaaggcttag'A28='tccttagatttcagttactatatttgacttacagtctttgagatttcccttacgattttgacttaaaatttagacgttagggcttatcagttatggattaatttagcttattttcga'A29='ggccaattccggtaggaaggtgatggcccgggggttcccgggaggatttaggctgacgggccggccatttcggtttagggagggccgggacgcgttagggc'A30='cgctaagcagctcaagctcagtcagtc
21、acgtttgccaagtcagtaatttgccaaagttaaccgttagctgacgctgaacgctaaacagtattagctgatgactcgta'A31='ttaaggacttaggctttagcagttactttagtttagttccaagctacgtttacgggaccagatgctagctagcaatttattatccgtattaggcttaccgtaggtttacgt'A32='gctaccgggcagtctttaacgtagctaccgtttagtttgggcccagccttgcggtgtttcggattaaattcgttgtcagtc
22、gctcttgggtttagtcattcccaaaagg'推荐精选A33='cagttagctgaatcgtttagccatttgacgtaaacatgattttacgtacgtaaattttagccctgacgtttagctaggaatttatgctgacgtagcgatcgactttagcac'A34='cggttagggcaaaggttggatttcgacccagggggaaagcccgggacccgaacccagggctttagcgtaggctgacgctaggcttaggttggaacccggaaa'A35='gcggaagggcgta
23、ggtttgggatgcttagccgtaggctagctttcgacacgatcgattcgcacacaggataaaagttaagggaccggtaagtcgcggtagcc'A36='ctagctacgaacgctttaggcgcccccgggagtagtcgttaccgttagtatagcagtcgcagtcgcaattcgcaaaagtccccagctttagccccagagtcgacg'A37='gggatgctgacgctggttagctttaggcttagcgtagctt tagggccccagtctgcaggaaatgcccaaaggaggcc
24、accgggtagatgccasagtgcaccgt'A38='aacttttagggcatttccagttttacgggttattttcccagttaaactttgcaccattttacgtgttacgatttacgtataatttgaccttattttggacactttagtttgggttac'A39='ttagggccaagtccccaggcaaggaattctgatccaagtccaatcacctacagtccaagtcaccgtttgcagctaccgtttaccgtacgttgcaagtcaaatccat'A40='ccattclgg
25、gtttatttacctctttattttttcccgagaccttaggtttaccgtactttttaacggtttacctttgaaatttttggactagcttaccctggatttaacggccagttt'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 &#
26、39;c' x1(i,3)=x1(i,3)+1; case 'g' x1(i,4)=x1(i,4)+1; end endend 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.39
27、64 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
28、.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
29、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
30、 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)
31、,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
32、9;,'none','marker','o','color','g');line(Y(index2,1),Y(index2,2),'linestyle','none','marker','*','color','r');hold onplot(center(1,1),center(1,2),'kpentagram','markersize', 15,'LineWidth',
33、2);plot(center(2,1),center(2,2),'ksquare','markersize', 15,'LineWidth',2);box on,title('模糊C均值聚类');推荐精选r11=0;r12=0t=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 endendsprintf('属于A类的DNA序列序号')t2sprintf(
34、9;属于B类的DNA序列序号')t1运行结果：属于A类的DNA序列序号t2 = Columns 1 through 15 0 0 3 0 5 0 7 0 9 0 0 0 0 14 15Columns 16 through 19 0 17 0 19ans =属于B类的DNA序列序号t1 =Columns 1 through 15 1 2 0 4 0 6 0 8 0 10 11 12 13 0 0Columns 16 through 20 16 0 18 0 20 推荐精选 未知的20个序列聚类图（ × 表示A类， 表示B类， ， 表示两类中心） （注：可编辑下载，若有不当之处，请指正，谢谢!） 推荐精选