DNA序列分类问题

DNA序列分类问题
摘要本文首先把问题中的DNA序列转换成氨基酸序列，对序列根据字符出现的频率进行特征提取，从而把DNA序列的结构转换为分析氨基酸的频率问题。

运用距离判别分析，Fisher判别，模糊聚类判别模型对DNA的分类进行了了讨论，其中方法一和方法二都各自做了检验，确定置信度均为95%，比较可靠，但方法三没有进行类似的检验，无法确定其正确率。

在最后的模型改进中提出题目中所给的人工序列都比较短，按统计特性将其分为A，B两类，和用其特征对很长的序列进行分类，可能会带来如下问题，即由于长序列的信息增加，A类或B类的特征表现应更加明显，但所使用的分类特征只是反应较短序列的规律，对长序列就不够明显了。

同时又提出可以考虑字母出现周期，序列熵值，序列片段之间具有一定的相关性等相关特征，在此方向上在进行研究，分类结果会更具有可信度。

关键字距离判别分析 Fisher判别模糊聚类
一、问题重述
DNA 序列是由4个字符a ，t ，c ，g 按一定顺序排列的序列，其中没有“断句”也没有标点符号，除了这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 类
二、问题分析
这是一个比较典型的分类问题, 为了表述的严格和方便, 我们用数学的方法来重述这个问题。

已知字母序列40321,,S S S S ，n i x x x S 21=,其中
},,,{g c t a x j ∈；有字符序列集合,,B A 满足Φ=⋂B A ,并当101≤≤i 时，A S i ∈；
当2011≤≤i 时，B S i ∈。

现要求考虑当4021≤≤i 时,i S 与集合A 及集合B 的关系。

在这里,问题的关键就是要从已知的分好类的20 个字母序列中提取用于分类的特征。

知道了这些特征, 我们就可以比较容易的对那些未标明类型的序列进行分类。

下面我们将首先对用于分类的标准问题进行必要的讨论。

分类的标准及评价
首先, 我们提取的特征应该满足以下两个条件: (1) 所取特征必须可以标志A 组和B 组。

也就是说, 我们利用这些特征应该可以很好的区分已经标示分类的20个序列。

这是比较显然的一个理由。

(2) 所取特征必须是有一定的实际意义的。

这一点是决不能被忽视的。

比如,如果不考虑模型的实际意义, 我们就可以以序列的开头字母为分类标准:已知在
B类中的十个序列都是以gt开始的,而已知在A类中10 个序列没有以gt开始的, 甚至以g开始的都没有。

显然这是满足上面的第一个条件的。

如果仅因此就认为这种特征是主要的,并简单的利用这个特征将所有待分类的序列分成两类, 显然是不甚合理的。

对于这样的一个复杂的分类问题, 需要考虑的因素很多,也是就说,可供我们使用的分类特征有许多。

如何从众多的因素中提取分类的主要因素, 是我们处理这个问题的困难之处。

上面的第一个条件是我们的分类方法所必须满足的, 可以看作是个限制条件；而第二个条件是我们在设计分类方法时必须考虑到的, 可以看作是对分类方法优劣的一种衡量, 是某种意义下的目标函数。

三、模型假设
假设一任意的DNA序列属于A类和B类的概率相等；
假设二不考虑氨基酸排列顺序对DNA序列所表示的生物信息；
假设三不考虑碱基序列的非编码区与编码区的区别；
假设四题目中所给的样本信息量足够大；
假设五假设变量服从多元正态分布；
假设六假设已知总体的协方差矩阵相等。

四、符号说明
x表示碱基A在序列中所占的百分比
a
x表示碱基T在序列中所占的百分比
t
x表示碱基C在序列中所占的百分比
c
x表示碱基G在序列中所占的百分比
g
x表示样品
i
μ样品均值向量
∑协方差矩阵
d马氏距离
L欧氏距离
G表示总体
i
w线形判别函数
五、模型建立及求解
问题一的模型建立及求解
我们考虑采用序列中4个碱基A,T,C,G 的含量百分比作为DNA 序列的特征。

具体规定如下：
T g c t a i x x x x X ),,,(=表示第个DNA 序列的特征向量。

补充：1=+++g c t a x x x x
我们根据字母出现的频率这一特征建立特征矩阵，用MATLAB 计算结果如下（程序见附录一）： 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.3545 0.2636 0.2455 0.1364 0.2909 0.5000 0.1182 0.0909 0.2182 0.5636 0.1455 0.0727 0.2000 0.5636 0.1727 0.0636 0.2743 0.3628 0.1947 0.1681 0.2885 0.2212 0.2404 0.2500 0.1765 0.1863 0.2549 0.3824 0.2087 0.4087 0.1913 0.1913 0.2476 0.2190 0.2286 0.3048 0.2193 0.3860 0.2105 0.1842 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.2411 0.3571 0.1786 0.2232 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.2427 0.2039 0.2136 0.3398 0.2286 0.2095 0.3048 0.2571 0.2136 0.2039 0.2524 0.3301 0.2222 0.4359 0.1709 0.1709 0.2736 0.2358 0.2830 0.2075
0.1983 0.4310 0.1983 0.1724
方法一：距离判别分析
定义：若x,y 是来自总体G 的两个样本。

总体的均值向量为μ，协方差矩阵为∑，则x,y 之间的马氏平方距离为：)()(),(12y x x G x d T -∑-==-μ，x 与G 的马氏平方距离为：)()(),(12μμ-∑-==-x x G x d T 。

因两总体的协方差矩阵相等，即有∑=∑=∑21，所以x 到两总体的马氏平方距离的差为：
1111121212221222),()(μμμμμμ----∑-∑+∑+∑-=--T
T
T
T
x x G x d G x d
记：111)(b x a x w T
+=，其中111111121,μμμ--∑-=∑=T b a
222)(b x a x w T
+=，其中212221221,μμμ--∑-=∑=T b a
求解，当)()(21x w x w ≥时，A x ∈，当)()(21x w x w <时，B x ∈
以上面求的字母频率作为新的矩阵，对其运用判别分析方法计算,MA TLAB 编程计算（程序见附录二），结果如下：
运算结果 oldclass =
1 1 1
2 1 1 1 1 1 1 2 2 2 2
2 2 2 2 2 2
newclass =
2 2 1 2 1 2 1 2 1 2 2 2 2 1 1 2 1 2 2 2
由程序回判结果可以看出，只有4判错，所以误判率为5%，回代误判率还是很低的，故此方法分类结果比较可靠，将分类结果整理为：
A 类：23，25，27，29，34，35，37
B 类：21，22，24，、26，28，30，31，32，33，36，38，39，40 方法二：Fisher 判别
费希尔判别的基本思想：借助于方差分析的思想，利用投影将n 元的数据投影到某一个方向，使得投影后组与组之间的差异尽可能的大，然后根据一定的判别规则对新样本的类别进行判断。

(1) 设总体为n G G G ,,,21 ；
(2) 由给出的样本信息计算各总体的样本均值向量i
x 和总平均向量x ； (3) 计算矩阵∑=--=k
i T
i
i
i x x x x n B 1]))([(及]))(([11
T
k i n j i ij i ij i
x x x x V ∑∑==--= ;
(4) 求出B V V 11,--；
(5) 求出B V 1-的最大特征值λ及相应的单位特征向量U ； (6) 计算出判别函数x U x T =)(ω；
(7))(i i x ωω=，将i ω从小到大排列，则相邻两类1,+i i G G 的阀值为
2
)1,(1
++=
+i i c i i y ωω；
(8) 若)1,()(),1(+<<-i i y x i i y c c ω，则i G x ∈。

Fisher 分类的判据为：
若L （X ，A ）〈 L （X ，B ），则判定x 为A 类； 若L （X ，A ）〉 L （X ，B ），则判定x 为B 类；
若L （X ，A ）= L （X ，B ），则判定x 为不可判类 .
用交叉确定估计的方法通过MATLAB 计算出20组数据如下(程序见附录三，依次变换数据计算，计算20
通过比较各个w 前10组数据只有当w 小于y 时，才判断正确，后10组数据只有当w 大于y 时，才判断正确。

由上表得
出第4组数据判断错误，其他数据均判断正确。

误判率为5%。

所以可以对21到40组DNA 进行比较准确的分类，通过MATLAB 计算（程序见附录四），结果如下：
y =
0.5000 w =
Columns 1 through 8
0.5002 0.5000 0.4997 0.5003 0.4998 0.5003 0.4997 0.5003
Columns 9 through 16
0.4994 0.5001 0.5001 0.5001 0.5001 0.4996 0.4997 0.5000
Columns 17 through 20
0.4998 0.5003 0.5001 0.5003
当w 小于y 时，样本属于A 类，当w 大于y 时，样本属于B 类，当w=y 时，无法分类，经比较
A 类：23，25，27，29，34，35，37
B 类：21，24，26，28，30，31，32，33，38，39，40 方法三：模糊聚类判别法
在自然科学或社会科学研究中，存在着许多定义不很严格或者说具有模糊性的概念。

这里所谓的模糊性，主要是指客观事物的差异在中间过渡中的不分明性，如某一生态条件对某种害虫、某种作物的存活或适应性可以评价为“有利、比较有利、不那么有利、不利”；灾害性霜冻气候对农业产量的影响程度为“较重、严重、很严重”，等等。

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

定义目标函数为
22
1
20
1)()(),(ik m ik i k d u V U J ∑
∑
===
表示各类中样本到聚类中心的加权距离平方和，权重是样本k x 对第i 类隶度的m 次方，聚类准则取为求),(V U J 的极小值)},((min){V U J
其中，][ik u U =为模糊分类矩阵,满足10≤≤ik u 和若5.0}m ax {>=ik ik u u ，则
∈k x 第j 类。

用MATLAB 编程，计算结果如下（程序见附录五）： index1 =
2 3 5 7 9 10 14 15 16 17 19 index2 =
1 4 6 8 11 1
2 1
3 18 20
整理得分类结果为：
A类：22，23，25，27，29，30，34，35，36，37，39
B类：21，24，26，28，31，32，33，38，40
六、结果分析
在模型中，方法一和方法二的误判率仅为5%，保证了95%的置信度。

同时注意到方法二中21-40号人工序列所有可能的氨基酸序列不可分性为10%，这说明系统的可靠性会随着序列个数和长度的增加额降低。

方法三没有进行显著性检验，所以其置信度没有保证。

七、模型评价及改进
题目中所给的人工序列都比较短，按统计特性将其分为A，B两类，和用其特征对很长的序列进行分类，可能会带来如下问题，即由于长序列的信息增加，A类或B类的特征表现应更加明显，但所使用的分类特征只是反应较短序列的规律，对长序列就不够明显了。

上面的三种方法只是根据字符出现频率这个作为特征，没有更多的去进行拓展，比如还有字母出现周期，序列熵值，序列片段之间具有一定的相关性等，在此方向上在进行研究，分类结果会更具有可信度。

八、参考文献
【1】张志涌MATLAB教程北京航空航天大学出版社2010.8 【2】姜启源谢金星叶俊编数学模型高等教育出版社2009. 5
附录
附录一
clear;
A1='aggcacggaaaaacgggaataacggaggaggacttggcacggcattacacggaggacgagg taaaggaggcttgtctacggccggaagtgaagggggatatgaccgcttgg';
A2='cggaggacaaacgggatggcggtattggaggtggcggactgttcggggaattattcggttt aaacgggacaaggaaggcggctggaacaaccggacggtggcagcaaagga';
A3='gggacggatacggattctggccacggacggaaaggaggacacggcggacatacacggcggc aacggacggaacggaggaaggagggcggcaatcggtacggaggcggcgga';
A4='atggataacggaaacaaaccagacaaacttcggtagaaatacagaagcttagatgcatatg ttttttaaataaaatttgtattattatggtatcataaaaaaaggttgcga';
A5='cggctggcggacaacggactggcggattccaaaaacggaggaggcggacggaggctacacc accgtttcggcggaaaggcggagggctggcaggaggctcattacggggag';
A6='atggaaaattttcggaaaggcggcaggcaggaggcaaaggcggaaaggaaggaaacggcgg
atatttcggaagtggatattaggagggcggaataaaggaacggcggcaca';
A7='atgggattattgaatggcggaggaagatccggaataaaatatggcggaaagaacttgtttt cggaaatggaaaaaggactaggaatcggcggcaggaaggatatggaggcg';
A8='atggccgatcggcttaggctggaaggaacaaataggcggaattaaggaaggcgttctcgct tttcgacaaggaggcggaccataggaggcggattaggaacggttatgagg';
A9='atggcggaaaaaggaaatgtttggcatcggcgggctccggcaactggaggttcggccatgg aggcgaaaatcgtgggcggcggcagcgctggccggagtttgaggagcgcg';
A10='tggccgcggaggggcccgtcgggcgcggatttctacaagggcttcctgttaaggaggtgg catccaggcgtcgcacgctcggcgcggcaggaggcacgcgggaaaaaacg';
A11='gttagatttaacgttttttatggaatttatggaattataaatttaaaaatttatattttt taggtaagtaatccaacgtttttattactttttaaaattaaatatttatt';
A12='gtttaattactttatcatttaatttaggttttaattttaaatttaatttaggtaagatga atttggttttttttaaggtagttatttaattatcgttaaggaaagttaaa';
A13='gtattacaggcagaccttatttaggttattattattatttggattttttttttttttttt tttaagttaaccgaattattttctttaaagacgttacttaatgtcaatgc';
A14='gttagtcttttttagattaaattattagattatgcagtttttttacataagaaaattttt ttttcggagttcatattctaatctgtctttattaaatcttagagatatta';
A15='gtattatatttttttatttttattattttagaatataatttgaggtatgtgtttaaaaaa aatttttttttttttttttttttttttttttttaaaatttataaatttaa';
A16='gttatttttaaatttaattttaattttaaaatacaaaatttttactttctaaaattggtc tctggatcgataatgtaaacttattgaatctatagaattacattattgat';
A17='gtatgtctatttcacggaagaatgcaccactatatgatttgaaattatctatggctaaaa accctcagtaaaatcaatccctaaacccttaaaaaacggcggcctatccc';
A18='gttaattatttattccttacgggcaattaattatttattacggttttatttacaattttt tttttttgtcctatagagaaattacttacaaaacgttattttacatactt';
A19='gttacattatttattattatccgttatcgataattttttacctcttttttcgctgagttt ttattcttactttttttcttctttatataggatctcatttaatatcttaa';
A20='gtatttaactctctttactttttttttcactctctacattttcatcttctaaaactgttt gatttaaacttttgtttctttaaggattttttttacttatcctctgttat';
A21='tttagctcagtccagctagctagtttacaatttcgacaccagtttcgcaccatcttaaat ttcgatccgtaccgtaatttagcttagatttggatttaaaggatttagattga';
A22='tttagtacagtagctcagtccaagaacgatgtttaccgtaacgtqacgtaccgtacgcta ccgttaccggattccggaaagccgattaaggaccgatcgaaaggg';
A23='cgggcggatttaggccgacggggacccgggattcgggacccgaggaaattcccggattaa ggtttagcttcccgggatttagggcccggatggctgggaccc';
A24='tttagctagctactttagctatttttagtagctagccagcctttaaggctagctttagct agcattgttctttattgggacccaagttcgacttttacgatttagttttgaccgt';
A25='gaccaaaggtgggctttagggacccgatgctttagtcgcagctggaccagttccccaggg tattaggcaaaagctgacgggcaattgcaatttaggcttaggcca';
A26='gatttactttagcatttttagctgacgttagcaagcattagctttagccaatttcgcatt tgccagtttcgcagctcagttttaacgcgggatctttagcttcaagctttttac';
A27='ggattcggatttacccggggattggcggaacgggacctttaggtcgggacccattaggag taaatgccaaaggacgctggtttagccagtccgttaaggcttag';
A28='tccttagatttcagttactatatttgacttacagtctttgagatttcccttacgattttg
acttaaaatttagacgttagggcttatcagttatggattaatttagcttattttcga';
A29='ggccaattccggtaggaaggtgatggcccgggggttcccgggaggatttaggctgacggg ccggccatttcggtttagggagggccgggacgcgttagggc';
A30='cgctaagcagctcaagctcagtcagtcacgtttgccaagtcagtaatttgccaaagttaa ccgttagctgacgctgaacgctaaacagtattagctgatgactcgta';
A31='ttaaggacttaggctttagcagttactttagtttagttccaagctacgtttacgggacca gatgctagctagcaatttattatccgtattaggcttaccgtaggtttagcgt';
A32='gctaccgggcagtctttaacgtagctaccgtttagtttgggcccagccttgcggtgtttc ggattaaattcgttgtcagtcgctctrtgggtttagtcattcccaaaagg';
A33='cagttagctgaatcgtttagccatttgacgtaaacatgattttacgtacgtaaattttag ccctgacgtttagctaggaatttatgctgacgtagcgatcgactttagcac';
A34='cggttagggcaaaggttggatttcgacccagggggaaagcccgggacccgaacccagggc tttagcgtaggctgacgctaggcttaggttggaacccggaaa';
A35='gcggaagggcgtaggtttgggatgcttagccgtaggctagctttcgacacgatcgattcg caccacaggataaaagttaagggaccggtaagtcgcggtagcc';
A36='ctagctacgaacgctttaggcgcccccgggagtagtcgttaccgttagtatagcagtcgc agtcgcaattcgcaaaagtccccagctttagccccagagtcgacg';
A37='gggatgctgacgctggttagctttaggcttagcgtagctttagggccccagtctgcagga aatgcccaaaggaggcccaccgggtagatgccasagtgcaccgt';
A38='aacttttagggcatttccagttttacgggttattttcccagttaaactttgcaccatttt acgtgttacgatttacgtataatttgaccttattttggacactttagtttgggttac';
A39='ttagggccaagtcccgaggcaaggaattctgatccaagtccaatcacgtacagtccaagt caccgtttgcagctaccgtttaccgtacgttgcaagtcaaatccat';
A40='ccattagggtttatttacctgtttattttttcccgagaccttaggtttaccgtacttttt aacggtttacctttgaaatttttggactagcttaccctggatttaacggccagttt';
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
附录二
附录二：
clear;
training=[ 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.3545 0.2636 0.2455 0.1364
0.2909 0.5000 0.1182 0.0909
0.2182 0.5636 0.1455 0.0727
0.2000 0.5636 0.1727 0.0636]; group=[1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2]'; sample=[ 0.2743 0.3628 0.1947 0.1681
0.2885 0.2212 0.2404 0.2500
0.1765 0.1863 0.2549 0.3824
0.2087 0.4087 0.1913 0.1913
0.2476 0.2190 0.2286 0.3048
0.2193 0.3860 0.2105 0.1842
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.2411 0.3571 0.1786 0.2232
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.2427 0.2039 0.2136 0.3398
0.2286 0.2095 0.3048 0.2571
0.2136 0.2039 0.2524 0.3301
0.2222 0.4359 0.1709 0.1709
0.2736 0.2358 0.2830 0.2075
0.1983 0.4310 0.1983 0.1724];
oldclass=classify(training,training,group)
newclass=classify(sample,training,group)
附录三
clear;
x1=[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.3514 0.1261 0.1261 0.3964
0.3514 0.1892 0.0991 0.3604
0.2342 0.1081 0.2342 0.4234
0.2793 0.1892 0.1622 0.3694
0.2072 0.1532 0.2072 0.4324
0.1818 0.1364 0.2727 0.4091]; x2=[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.3545 0.2636 0.2455 0.1364
0.2909 0.5000 0.1182 0.0909
0.2182 0.5636 0.1455 0.0727
0.2000 0.5636 0.1727 0.0636]; b1=sum(x1);%计算每列的总和
b2=sum(x2);
m1=size(x1,1);%计算矩阵的行数
m2=size(x2,1);
k1=size(x1,2);%计算矩阵的列数
k2=size(x2,2);
x11=b1/m1;%计算均值向量
x21=b2/m2;
X=(x11+x21)/2;%计算总平均向量
B=m1*(x11-X)'*(x11-X)+m2*(x21-X)'*(x21-X); c1=0;
c2=0;
for i=1:m1
c1=c1+(x1(i,:)-x11)'*(x1(i,:)-x11); end
for i=1:m2
c2=c2+(x2(i,:)-x21)'*(x2(i,:)-x21);
end
v=c1+c2;
E=eye(k1);
V=E/v;
V1=V*B;
[C,u]=eig(V1);
eigenvalue=diag(u);
lamda=eigenvalue(1);
u_lamda=C(:,1);
W1=u_lamda'*x11';
W2=u_lamda'*x21';
y=(m1*W1+m2*W2)/(m1+m2)
z1=[0.2973 0.1351 0.1712 0.3964]; w1=u_lamda'*z1'
附录四
clear;
x1=[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]; x2=[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.3545 0.2636 0.2455 0.1364
0.2909 0.5000 0.1182 0.0909
0.2182 0.5636 0.1455 0.0727
0.2000 0.5636 0.1727 0.0636]; b1=sum(x1);%计算每列的总和
b2=sum(x2);
m1=size(x1,1);%计算矩阵的行数
m2=size(x2,1);
k1=size(x1,2);%计算矩阵的列数
k2=size(x2,2);
x11=b1/m1;%计算均值向量
x21=b2/m2;
X=(x11+x21)/2;%计算总平均向量
B=m1*(x11-X)'*(x11-X)+m2*(x21-X)'*(x21-X); c1=0;c2=0;
for i=1:m1
c1=c1+(x1(i,:)-x11)'*(x1(i,:)-x11); end
for i=1:m2
c2=c2+(x2(i,:)-x21)'*(x2(i,:)-x21); end
v=c1+c2;
E=eye(k1);
V=E/v;
V1=V*B;
[C,u]=eig(V1);
eigenvalue=diag(u);
lamda=eigenvalue(1);
u_lamda=C(:,1);
W1=u_lamda'*x11';
W2=u_lamda'*x21';
y=(m1*W1+m2*W2)/(m1+m2)
z=[0.2743 0.3628 0.1947 0.1681
0.2885 0.2212 0.2404 0.2500
0.1765 0.1863 0.2549 0.3824
0.2087 0.4087 0.1913 0.1913
0.2476 0.2190 0.2286 0.3048
0.2193 0.3860 0.2105 0.1842
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.2411 0.3571 0.1786 0.2232
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.2427 0.2039 0.2136 0.3398
0.2286 0.2095 0.3048 0.2571
0.2136 0.2039 0.2524 0.3301
0.2222 0.4359 0.1709 0.1709
0.2736 0.2358 0.2830 0.2075
0.1983 0.4310 0.1983 0.1724]; for i=1:20
w(i)=u_lamda'*z(i,:)';
end
w
附录五
clear;
A1='aggcacggaaaaacgggaataacggaggaggacttggcacggcattacacggaggacgagg taaaggaggcttgtctacggccggaagtgaagggggatatgaccgcttgg';
A2='cggaggacaaacgggatggcggtattggaggtggcggactgttcggggaattattcggttt aaacgggacaaggaaggcggctggaacaaccggacggtggcagcaaagga';
A3='gggacggatacggattctggccacggacggaaaggaggacacggcggacatacacggcggc aacggacggaacggaggaaggagggcggcaatcggtacggaggcggcgga';
A4='atggataacggaaacaaaccagacaaacttcggtagaaatacagaagcttagatgcatatg ttttttaaataaaatttgtattattatggtatcataaaaaaaggttgcga';
A5='cggctggcggacaacggactggcggattccaaaaacggaggaggcggacggaggctacacc accgtttcggcggaaaggcggagggctggcaggaggctcattacggggag';
A6='atggaaaattttcggaaaggcggcaggcaggaggcaaaggcggaaaggaaggaaacggcgg atatttcggaagtggatattaggagggcggaataaaggaacggcggcaca';
A7='atgggattattgaatggcggaggaagatccggaataaaatatggcggaaagaacttgtttt cggaaatggaaaaaggactaggaatcggcggcaggaaggatatggaggcg';
A8='atggccgatcggcttaggctggaaggaacaaataggcggaattaaggaaggcgttctcgct tttcgacaaggaggcggaccataggaggcggattaggaacggttatgagg';
A9='atggcggaaaaaggaaatgtttggcatcggcgggctccggcaactggaggttcggccatgg aggcgaaaatcgtgggcggcggcagcgctggccggagtttgaggagcgcg';
A10='tggccgcggaggggcccgtcgggcgcggatttctacaagggcttcctgttaaggaggtgg catccaggcgtcgcacgctcggcgcggcaggaggcacgcgggaaaaaacg';
A11='gttagatttaacgttttttatggaatttatggaattataaatttaaaaatttatattttt taggtaagtaatccaacgtttttattactttttaaaattaaatatttatt';
A12='gtttaattactttatcatttaatttaggttttaattttaaatttaatttaggtaagatga atttggttttttttaaggtagttatttaattatcgttaaggaaagttaaa';
A13='gtattacaggcagaccttatttaggttattattattatttggattttttttttttttttt tttaagttaaccgaattattttctttaaagacgttacttaatgtcaatgc';
A14='gttagtcttttttagattaaattattagattatgcagtttttttacataagaaaattttt ttttcggagttcatattctaatctgtctttattaaatcttagagatatta';
A15='gtattatatttttttatttttattattttagaatataatttgaggtatgtgtttaaaaaa aatttttttttttttttttttttttttttttttaaaatttataaatttaa';
A16='gttatttttaaatttaattttaattttaaaatacaaaatttttactttctaaaattggtc tctggatcgataatgtaaacttattgaatctatagaattacattattgat';
A17='gtatgtctatttcacggaagaatgcaccactatatgatttgaaattatctatggctaaaa accctcagtaaaatcaatccctaaacccttaaaaaacggcggcctatccc';
A18='gttaattatttattccttacgggcaattaattatttattacggttttatttacaattttt tttttttgtcctatagagaaattacttacaaaacgttattttacatactt';
A19='gttacattatttattattatccgttatcgataattttttacctcttttttcgctgagttt ttattcttactttttttcttctttatataggatctcatttaatatcttaa';
A20='gtatttaactctctttactttttttttcactctctacattttcatcttctaaaactgttt gatttaaacttttgtttctttaaggattttttttacttatcctctgttat';
A21='tttagctcagtccagctagctagtttacaatttcgacaccagtttcgcaccatcttaaat ttcgatccgtaccgtaatttagcttagatttggatttaaaggatttagattga';
A22='tttagtacagtagctcagtccaagaacgatgtttaccgtaacgtqacgtaccgtacgcta ccgttaccggattccggaaagccgattaaggaccgatcgaaaggg';
A23='cgggcggatttaggccgacggggacccgggattcgggacccgaggaaattcccggattaa ggtttagcttcccgggatttagggcccggatggctgggaccc';
A24='tttagctagctactttagctatttttagtagctagccagcctttaaggctagctttagct agcattgttctttattgggacccaagttcgacttttacgatttagttttgaccgt';
A25='gaccaaaggtgggctttagggacccgatgctttagtcgcagctggaccagttccccaggg tattaggcaaaagctgacgggcaattgcaatttaggcttaggcca';
A26='gatttactttagcatttttagctgacgttagcaagcattagctttagccaatttcgcatt tgccagtttcgcagctcagttttaacgcgggatctttagcttcaagctttttac';
A27='ggattcggatttacccggggattggcggaacgggacctttaggtcgggacccattaggag taaatgccaaaggacgctggtttagccagtccgttaaggcttag';
A28='tccttagatttcagttactatatttgacttacagtctttgagatttcccttacgattttg acttaaaatttagacgttagggcttatcagttatggattaatttagcttattttcga';
A29='ggccaattccggtaggaaggtgatggcccgggggttcccgggaggatttaggctgacggg ccggccatttcggtttagggagggccgggacgcgttagggc';
A30='cgctaagcagctcaagctcagtcagtcacgtttgccaagtcagtaatttgccaaagttaa ccgttagctgacgctgaacgctaaacagtattagctgatgactcgta';
A31='ttaaggacttaggctttagcagttactttagtttagttccaagctacgtttacgggacca gatgctagctagcaatttattatccgtattaggcttaccgtaggtttagcgt';
A32='gctaccgggcagtctttaacgtagctaccgtttagtttgggcccagccttgcggtgtttc ggattaaattcgttgtcagtcgctctrtgggtttagtcattcccaaaagg';
A33='cagttagctgaatcgtttagccatttgacgtaaacatgattttacgtacgtaaattttag ccctgacgtttagctaggaatttatgctgacgtagcgatcgactttagcac';
A34='cggttagggcaaaggttggatttcgacccagggggaaagcccgggacccgaacccagggc tttagcgtaggctgacgctaggcttaggttggaacccggaaa';
A35='gcggaagggcgtaggtttgggatgcttagccgtaggctagctttcgacacgatcgattcg caccacaggataaaagttaagggaccggtaagtcgcggtagcc';
A36='ctagctacgaacgctttaggcgcccccgggagtagtcgttaccgttagtatagcagtcgc agtcgcaattcgcaaaagtccccagctttagccccagagtcgacg';
A37='gggatgctgacgctggttagctttaggcttagcgtagctttagggccccagtctgcagga aatgcccaaaggaggcccaccgggtagatgccasagtgcaccgt';
A38='aacttttagggcatttccagttttacgggttattttcccagttaaactttgcaccatttt acgtgttacgatttacgtataatttgaccttattttggacactttagtttgggttac';
A39='ttagggccaagtcccgaggcaaggaattctgatccaagtccaatcacgtacagtccaagt caccgtttgcagctaccgtttaccgtacgttgcaagtcaaatccat';
A40='ccattagggtttatttacctgtttattttttcccgagaccttaggtttaccgtacttttt aacggtttacctttgaaatttttggactagcttaccctggatttaacggccagttt';
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=y(21:40,:);
[center,U1,obj_fcn]=fcm(Y,2);%Y(data)要聚类的数据函数，每一行为一个样本maxU1=max(U1);%2：聚类数 center：最终的聚类中心矩阵，其每一行为聚类中心的坐标值
index1=find(U1(1,:)==maxU1) %U：最终的模糊分区矩阵
index2=find(U1(2,:)==maxU1) %obj_fcn：在迭代过程中的目标函数值。