DNA序列分类 作为研究DNA序列结构的尝试,提出以下对序列集合进行分类的问题:有20个已知类别的人工制造序列,其中序列标号1-10为A类,11-20为B类。请从中提取特征,构造分类方法,并用这些已知类别的序列,衡量你的方法是否足够好。然后用你认为满意的方法,对另外20个未标明类别的人工序列(标号21-40)进行分类,判断哪些属于A类,哪些属于B类,哪些既不属于A类又不属于B类。
(一)问题分析
采用DNA序列中4个字符的含量百分比对DNA序列进行分类。
(二)模型建立
(1)BP神经网络结构的确定
我们选取三层结构的BP神经网络模型。
输入层:将所提取的特征作为输入,即DNA序列中a,t,c,g的含量百分比作为BP神经网络的输入。显然,输入层有4个节点。
输出层:BP神经网络的输出为DNA序列的分类结果,将A类DNA序列的输出定义为[1,0],B类DNA序列的输出定义为[0,1],因此,输出层有2个节点。
隐含层:为确定隐含层节点个数l,我们参考下列公式:
l=sqrt(n+m)+a;
其中n为输入层节点数,m为输出层节点个数,a为1-10之间的常数。我们这里确定隐含层个数为11。
综上所述,建立了一个4-11-2结构的三层BP神经网络,并选择双型s型函数(tansig)作为隐含层的传递函数,线性函数(purelin)作为输出层的传递函数,变学习动量梯度下降算法(traingdx)作为训练函数。
(2)训练数据及测试数据的确定
训练数据用来训练BP神经网络,测试数据用来测试网络的分类能力。但由于已知类别的DNA序列只有20条(标号1-20),比较少,因此,我们将这20条数据即作为训练数据,又作为测试数据。最后,用训练好的BP神经网络对标号为21-40的DNA序列进行分类。
(3)BP神经网络的训练
(三)模型求解
利用MATLAB编程求解,在分类时某些DNA序列有时属于A类,有时属于B类。为此将程序运行100次(书上是1000次,电脑慢,就跑100次,效果也很好),统计分类结果:
DNA序号 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
A类 12 55 100 17 92 20 99 1 99 28 27 50 28 99 99 47 95 9 24 12
B类 88 45 0 83 8 80 1 99 1 72 73 50 72 1 1 53 5 91 76 88
从上表可以看出,DNA序列被分到A类或B类的次数明显的多,就可认为该DNA序列属于那一类。但是某些DNA序列被分到A、B类的次数非常相近(如标号22),因此,这些DNA序列即不属于A类又不属于B类,即无法用BP神经网络分类,需要作进一步分析。根据以上分析,我们得到最终的分类结果:
A类:23、25、27、29、34、35、37;
B类:21、24、26、28、30、31、33、38、39、40;
即不属于A类又不属于B类:22、32、36。
MATLAB程序:
clc clear all %%=========================统计字符个数以及含量==================================== fid=fopen('exp12_4_3.txt','r') %读取数据文件,返回文件标识符,文件打开成功,fid为正数,否则为-1。 i=1; %计数 while (~feof(fid)) %reof测试文件是否到最后一行,到最后一行返回1,否则为0 data=fgetl(fid); %fgetl表示读取文件的一行,不包括换行符号 a=length(find(data=='a')); %统计字符a的个数 t=length(find(data=='t')); %统计字符t的个数 c=length(find(data=='c')); %统计字符c的个数 g=length(find(data=='g')); %统计字符g的个数 e=length(find(data~='a'&data~='c'&data~='t'&data~='g')); %统计其它字符的个数 DNA_num(i,:)=[a t c g e a+c+t+g+e]; %将字符个数放到DNA_num矩阵中 DNA_HanL(i,:)=[a/(a+c+t+g) t/(a+c+t+g) c/(a+c+t+g) g/(a+c+t+g)]; %计算a,t,c,g字符的含量百分比 i=i+1; %文件行数加1 end fclose(fid); %关闭文件 %%=====================BP神经网络训练========================================== [n,m]=size(DNA_HanL); for i=1:20 %定义已知类DNA序列的输出 if i<=10 output(i,:)=[1,0]; %标号1-10为A类,输出为[1,0] else output(i,:)=[0,1]; %标号11-20为A类,输出为[0,1] end end train_output=output'; %神经网络训练的输出 train_input=DNA_HanL(1:20,:)'; %神经网络训练的输入 for LL=1:10 %程序运行1000次时,设置为1:1000 in_num=4; %输入层节点个数 mid_num=11; %隐含层节点个数 out_num=2; %输出层节点个数 TF1='tansig';TF2='purelin'; %TF1为隐含层传递函数,TF2为输出层传递函数 net=newff(minmax(train_input),[mid_num,out_num],{TF1,TF2}); %创建BP神经网络 net.trainFcn='traingdx'; %训练函数,变学习动量梯度下降算法 net.trainParam.epochs=5000; %以下为训练参数设置 net.trainParam.lr=0.1; %学习速率 net.trainParam.mc=0.75; %附加动量因子 net.trainParam.goal=0.001; % 训练目标最小误差 net=train(net,train_input,train_output); %网络训练 an=sim(net,train_input); %网络测试,此处测试数据即训练数据 for i=1:20 %测试分类结果统计 output_test_fore(i)=find(an(:,i)==max(an(:,i))); %1表示分到A类,2表示分到B类 output1(i)=find(train_output(:,i)==max(train_output(:,i))); end error=output_test_fore-output1; %BP网络分类误差 sim_input=DNA_HanL(21:40,:)'; %待分类数据 anan=sim(net,sim_input); %网络仿真,返回预测结果 for i=1:20 %预测分类结果统计 output_sim_fore(i)=find(anan(:,i)==max(anan(:,i))); %1表示分到A类,2表示分到B类 end out(LL,:)=output_sim_fore; %预测分类结果 end [nn,mm]=size(out); for ii=1:mm a=length(find(out(:,ii)==1)); b=length(find(out(:,ii)==2)); ff(ii,:)=[ii+20 a b]; end ff=ff'
Txt数据:
1.aggcacggaaaaacgggaataacggaggaggacttggcacggcattacacggaggacgaggtaaaggaggcttgtctacggccggaagtgaagggggatatgaccgcttgg 2.cggaggacaaacgggatggcggtattggaggtggcggactgttcggggaattattcggtttaaacgggacaaggaaggcggctggaacaaccggacggtggcagcaaagga 3.gggacggatacggattctggccacggacggaaaggaggacacggcggacatacacggcggcaacggacggaacggaggaaggagggcggcaatcggtacggaggcggcgga 4.atggataacggaaacaaaccagacaaacttcggtagaaatacagaagcttagatgcatatgttttttaaataaaatttgtattattatggtatcataaaaaaaggttgcga 5.cggctggcggacaacggactggcggattccaaaaacggaggaggcggacggaggctacaccaccgtttcggcggaaaggcggagggctggcaggaggctcattacggggag 6.atggaaaattttcggaaaggcggcaggcaggaggcaaaggcggaaaggaaggaaacggcggatatttcggaagtggatattaggagggcggaataaaggaacggcggcaca 7.atgggattattgaatggcggaggaagatccggaataaaatatggcggaaagaacttgttttcggaaatggaaaaaggactaggaatcggcggcaggaaggatatggaggcg 8.atggccgatcggcttaggctggaaggaacaaataggcggaattaaggaaggcgttctcgcttttcgacaaggaggcggaccataggaggcggattaggaacggttatgagg 9.atggcggaaaaaggaaatgtttggcatcggcgggctccggcaactggaggttcggccatggaggcgaaaatcgtgggcggcggcagcgctggccggagtttgaggagcgcg 10.tggccgcggaggggcccgtcgggcgcggatttctacaagggcttcctgttaaggaggtggcatccaggcgtcgcacgctcggcgcggcaggaggcacgcgggaaaaaacg 11.gttagatttaacgttttttatggaatttatggaattataaatttaaaaatttatattttttaggtaagtaatccaacgtttttattactttttaaaattaaatatttatt 12.gtttaattactttatcatttaatttaggttttaattttaaatttaatttaggtaagatgaatttggttttttttaaggtagttatttaattatcgttaaggaaagttaaa 13.gtattacaggcagaccttatttaggttattattattatttggattttttttttttttttttttaagttaaccgaattattttctttaaagacgttacttaatgtcaatgc 14.gttagtcttttttagattaaattattagattatgcagtttttttacataagaaaatttttttttcggagttcatattctaatctgtctttattaaatcttagagatatta 15.gtattatatttttttatttttattattttagaatataatttgaggtatgtgtttaaaaaaaatttttttttttttttttttttttttttttttaaaatttataaatttaa 16.gttatttttaaatttaattttaattttaaaatacaaaatttttactttctaaaattggtctctggatcgataatgtaaacttattgaatctatagaattacattattgat 17.gtatgtctatttcacggaagaatgcaccactatatgatttgaaattatctatggctaaaaaccctcagtaaaatcaatccctaaacccttaaaaaacggcggcctatccc 18.gttaattatttattccttacgggcaattaattatttattacggttttatttacaattttttttttttgtcctatagagaaattacttacaaaacgttattttacatactt 19.gttacattatttattattatccgttatcgataattttttacctcttttttcgctgagtttttattcttactttttttcttctttatataggatctcatttaatatcttaa 20.gtatttaactctctttactttttttttcactctctacattttcatcttctaaaactgtttgatttaaacttttgtttctttaaggattttttttacttatcctctgttat 21.tttagctcagtccagctagctagtttacaatttcgacaccagtttcgcaccatcttaaatttcgatccgtaccgtaatttagcttagatttggatttaaaggatttagattga 22.tttagtacagtagctcagtccaagaacgatgtttaccgtaacgtacgtaccgtacgctaccgttaccggattccggaaagccgattaaggaccgatcgaaaggg 23.cgggcggatttaggccgacggggacccgggattcgggacccgaggaaattcccggattaaggtttagcttcccgggatttagggcccggatggctgggaccc 24.tttagctagctactttagctatttttagtagctagccagcctttaaggctagctttagctagcattgttctttattgggacccaagttcgacttttacgatttagttttgaccgt 25.gaccaaaggtgggctttagggacccgatgctttagtcgcagctggaccagttccccagggtattaggcaaaagctgacgggcaattgcaatttaggcttaggcca 26.gatttactttagcatttttagctgacgttagcaagcattagctttagccaatttcgcatttgccagtttcgcagctcagttttaacgcgggatctttagcttcaagctttttac 27.ggattcggatttacccggggattggcggaacgggacctttaggtcgggacccattaggagtaaatgccaaaggacgctggtttagccagtccgttaaggcttag 28.tccttagatttcagttactatatttgacttacagtctttgagatttcccttacgattttgacttaaaatttagacgttagggcttatcagttatggattaatttagcttattttcga 29.ggccaattccggtaggaaggtgatggcccgggggttcccgggaggatttaggctgacgggccggccatttcggtttagggagggccgggacgcgttagggc 30.cgctaagcagctcaagctcagtcagtcacgtttgccaagtcagtaatttgccaaagttaaccgttagctgacgctgaacgctaaacagtattagctgatgactcgta 31.ttaaggacttaggctttagcagttactttagtttagttccaagctacgtttacgggaccagatgctagctagcaatttattatccgtattaggcttaccgtaggtttagcgt 32.gctaccgggcagtctttaacgtagctaccgtttagtttgggcccagccttgcggtgtttcggattaaattcgttgtcagtcgctcttgggtttagtcattcccaaaagg 33.cagttagctgaatcgtttagccatttgacgtaaacatgattttacgtacgtaaattttagccctgacgtttagctaggaatttatgctgacgtagcgatcgactttagcac 34.cggttagggcaaaggttggatttcgacccagggggaaagcccgggacccgaacccagggctttagcgtaggctgacgctaggcttaggttggaacccggaaa 35.gcggaagggcgtaggtttgggatgcttagccgtaggctagctttcgacacgatcgattcgcaccacaggataaaagttaagggaccggtaagtcgcggtagcc 36.ctagctacgaacgctttaggcgcccccgggagtagtcgttaccgttagtatagcagtcgcagtcgcaattcgcaaaagtccccagctttagccccagagtcgacg 37.gggatgctgacgctggttagctttaggcttagcgtagctttagggccccagtctgcaggaaatgcccaaaggaggcccaccgggtagatgccasagtgcaccgt 38.aacttttagggcatttccagttttacgggttattttcccagttaaactttgcaccattttacgtgttacgatttacgtataatttgaccttattttggacactttagtttgggttac 39.ttagggccaagtcccgaggcaaggaattctgatccaagtccaatcacgtacagtccaagtcaccgtttgcagctaccgtttaccgtacgttgcaagtcaaatccat 40.ccattagggtttatttacctgtttattttttcccgagaccttaggtttaccgtactttttaacggtttacctttgaaatttttggactagcttaccctggatttaacggccagttt