基于模糊聚类及BP神经网络的DNA序列分类

1基于模糊聚类及BP神经网络的DNA序列分类8队建模：张初旭编程：朱马写作:朱柳2摘要本文通过对已知分类的DNA序列的分析，对未标明类别的第21-40组以及182组DNA序列进行了分类。首先运用4种常用的模糊聚类方法：相关系数法、夹角余弦法、海明距离法、欧氏距离法对第1-20组序列分类。根据分类结果的正确率，发现相关系数法的分类效果很不理想，所以将除相关系数法外的3种方法所得到的分类结果进行综合，分别得到了模糊聚类方法在以碱基含量和氨基酸含量为序列特征下的分类结果及正确率。再构造三层BP神经网络，应用Matlab神经网络工具箱，将已知分类的1-20组DNA序列中的1—6,11—16组作为学习样本确定较优的网络参数，并以7—10,17—20组序列作为检验样本，检验此神经网络算法在以碱基含量和氨基酸含量为序列特征下的分类结果及正确率。然后将这四种不同聚类方案的分类效果进行分析比较，选出正确率较高的方案对第21-40组以及182组DNA序列进行了分类。最后对分类结果进行了合理性分析，客观地指出了模型的优缺点，还针对建模过程中简化、忽略的因素给出了模型改进的方向。关键词DNA序列分类碱基含量氨基酸含量模糊聚类BP神经网络3一、问题重述2000年6月，人类基因组计划中DNA全序列草图完成，预计2001年可以完成精确的全序列图，此后人类将拥有一本记录着自身生老病死及遗传进化的全部信息的“天书”。这本大自然写成的“天书”是由4个字符A，T，C，G按一定顺序排成的长约30亿的序列，其中没有“断句”也没有标点符号，这4个字符表示4种碱基。破译这部世界上最巨量信息的“天书”要研究DNA全序列具有什么结构，由这4个字符排成的看似随机的序列中隐藏着什么规律，这是生物信息学（Bioinformatics）最重要的课题之一。人类发现了DNA序列中的一些规律性和结构。例如，在全序列中有一些是用于编码蛋白质的序列片段，即由这4个字符组成的64种不同的3字符串，其中大多数用于编码构成蛋白质的20种氨基酸。又例如，在不用于编码蛋白质的序列片段中，A和T的含量特别多些，于是以某些碱基特别丰富作为特征去研究DNA序列的结构也取得了一些结果。此外，利用统计的方法还发现序列的某些片段之间具有相关性，等等。这些发现让人们相信，DNA序列中存在着局部的和全局性的结构，充分发掘序列的结构对理解DNA全序列是十分有意义的。目前在这项研究中最普通的思想是省略序列的某些细节，突出特征，然后将其表示成适当的数学对象。这种被称为粗粒化和模型化的方法往往有助于研究规律性和结构。需完成的工作是：1.从题中所给的已知类别的人工制造序列（序列标号1—10为A类，11-20为B类）中提取特征，构造分类方法，并用这些已知类别的序列衡量本文所采用的方法是否足够好。2.选取最适宜的方法对另外20个未标明类别的人工序列（标号21-40）进行分类。要详细描述所采用的方法，给出计算程序，并准确注明所使用的现成的分类方法的名称。3.用本文所采用的方法对数据文件Nat-model-data中给出的182个较长的自然DNA序列进行分类，给出分类结果。二、问题分析2.1背景知识（1）转录：以DNA双链中的一条为模板合成mRNA的过程。（2）翻译：将mRNA中的碱基序列翻译为蛋白质的氨基酸序列的过程。（3）密码子：mRNA链上决定一个氨基酸的相邻的三个碱基叫做一个“密码子”。共有64种密码子，其中有61种能合成氨基酸的密码子(包括起始密码子)及3个终止密码子，由它们决定多肽链的氨基酸种类和排列顺序的特异性以及翻译的起始和终止。（4）起始密码子：mRNA翻译起始时的第一个密码子。4（5）由61种能合成氨基酸的密码子能够合成20中氨基酸。（6）DNA，mRNA，密码子及氨基酸的关系：蛋白质合成氨基酸代表密码子翻译转录mRNADNA（7）密码子与氨基酸对照表1.表1.密码子与氨基酸对照表2.2问题的分析要构造对题中所给的1-20条人工制造序列的分类方法，就要首先确定从序列中提取何种特征来进行分类。而序列特征需要满足以下两个条件：1.可以标志A组和B组；2.有一定的生物学意义。其中，第二个条件可以看做是在分类正确率达到要求后，对分类方法是否有实用性的一种衡量标准。运用excel对题中已分类序列（第1-20组）进行分析，得出碱基含量分布图（如图1）：500.10.20.30.40.50.60.712345678910111213141516171819A0.2973C0.1712G0.3964T0.1351图1碱基含量分布表发现在不同段的DNA中，每个碱基出现的概率不同，从大体分布来看，序列所以考虑将4种碱基各自的含量作为序列的特征。又因为有科学研究结果表明：在不用于编码蛋白质的序列片断中，A和T的含量特别多些,因此，以碱基含量作为特征去研究DNA序列的分类是具有一定的生物学意义的。而这种只考虑了碱基含量的方法并没有考虑到碱基排列顺序上的不同而造成翻译产生的氨基酸不同。如表1中，序列(UCA)与序列(CUA)有着相同的碱基含量，却因排列顺序的不同而可以分别转录成丝氨酸及亮氨酸这两种不同的氨基酸。因此，以氨基酸含量作为序列的特征更具有生物学上的意义。由于题中所给DNA序列并未说明是全部的还是只截取了DNA链的某部分，因此无法确定将mRNA翻译成氨基酸的起始点，所以将所给序列的碱基依次作为密码子的起始点更为合理，如对给定序列aggcacgg，密码子依次为agg，ggc，gca…。确定了序列的特征，就可以建立分别以碱基含量和氨基酸含量为数据源的分类模型。下面要考虑分类方法的选择。由于对DNA序列分类后，期望该分类在生物学上具有一定的意义，也就是分类后的两部分DNA序列在被转录翻译成多肽或蛋白质后，能体现出不同的生物学性状和功能，而性状和功能的体现程度可能会所不同，因此考虑采用模糊聚类分析方法。又因为对DNA序列的分类与诸多生物学因素有关，是复杂且不确定的非线性系统，因此考虑用神经网络来解决。而目前,在神经网络中应用最多的是BP网络，所以决定采用BP神经网络方法对DNA序列进行分类。这样，对两种不同的数据源分别采用两种不同的分类方法就构成了四种分类方案，然后就可以选取正确率最高的一种方案对未标明类别的第21-40组DNA序列以及第二问中的182组自然DNA序列进行分类。2.3工作流程图本文将根据下面的工作流程图展开对DNA序列的分类工作。6图1.工作流程图三、基本假设1.转录及翻译过程中不发生基因突变，碱基缺失错位等情况。2.uaa,uag,uga这三种密码子为终止密码子，即不能被翻译成氨基酸。为简便起见，将这三种命名为终止氨基酸，所以假设由64种密码子可以生成21种氨基酸。3.各序列中DNA碱基构成三联体（密码子）的起始位置和排列顺序不影响分类的结果。4.较长的182个自然序列与已知类别的20个样本序列具有共同的特征。四、变量定义及说明iS：第i组DNA序列ijx：第i个DNA序列的第j个指标7ijr：样本is与js之间的相似程度衡量指标：置信水平H：使所有]1,0[ijr),,2,1,(nji的确定常数E：使得所有),,2,1,(]1,0[njirij的确定常数K：实数值样本个数五、模型建立及求解5.1基于碱基含量特征分类的模型首先，计算出题中所给的序列1-20中的A,C，G,T的含量百分比如表2：ACGT10.29730.17120.39640.135120.27030.16220.41440.153230.27030.21620.45050.063140.42340.10810.18020.288350.23420.23420.42340.108160.35140.12610.39640.126170.35140.09910.36040.189280.27930.16220.36940.189290.20720.20720.43240.1532100.18180.27270.40910.1364110.35450.04550.10000.5000120.32730.02730.14550.5000130.25450.10000.12730.5182140.30000.08180.11820.5000150.29090.00000.06360.6455160.36360.08180.09090.4636170.35450.24550.13640.2636180.29090.11820.09090.5000190.21820.14550.07270.5636200.20000.17270.06360.5636表2.序列1-20中的碱基含量表5.1.1采用模糊聚类方法对序列1-20分类为了表述的严格和方便，用数学方法来重述此问题。已知字母序列1S,2S,...,40S，injixxxxxS......321，其中gctaxj,,,；有字符序列A，B，满足BA，并当101i时，ASi；当2011i时，BSi。现要考虑当4021i时，iS与集合A及集合B的关系。(1)原始数据标准化8首先对样本数据进行预处理，并将数据压缩到[0,1]闭区间内。niijjxnx11（1）nijijjxxnS12)(1（2）其中：ijx表示第i个DNA序列的第j个指标jjijijSxxx'（3）minmaxmin''''jjjijijxxxxx（4）其中，min'jx，max'jx分别表示jx1'，jx2'，jx3'…njx'中的最小值和最大值。当min''jijxx时，则0ijx；当max''jijxx时，则1ijx。将N个样本的第j个指标的平均值公式（1）及标准差公式（2）带入原始数据标准化公式（3），即可得到标准化数据。然后再运用如下极值标准化公式（4），将公式（3）得到标准化数据压缩到[0.1]内，得到原始数据标准化并压缩到[0,1]范围后的输出数见表3.表3.标准化后的第1-20组DNA序列样本指标ACGT10.48150.63330.86050.125020.37040.60000.90700.156230.37040.80001.00000.000041.00000.40000.30230.390650.22220.86670.93020.078160.70370.46670.86050.109470.70370.36670.76740.218780.40740.60000.79070.218790.11110.76670.95350.1562100.00001.00000.88370.1250110.70370.16670.09300.7500120.59260.10000.20930.7500130.29630.36670.16280.7813140.48150.30000.13950.7500150.44440.00000.00001.0000160.74070.30000.06980.6875170.70370.90000.18600.3437180.44440.43330.06980.75009190.14810.53330.02330.8594200.07410.633300.8594(2)模糊相似矩阵的确定引入模糊相似矩阵如下：mmmmmmrrrrrrrrrR212221211211其中，ijr表示样本is与js之间的相似程度，当ijr接近于1时，表明这两个样本相似程度较高。对应于本文中，20m，is，js分别表示第i个和第j个DNA样本序列。由于模糊相似矩阵R的确定方法有很多，经过对数量积法、相关系数法等11种常用方法的演算，本文从中选取了分类正确率较高的4种方法进行详细阐述并给出了分类结果。a.相关系数法由公式),,2,1,(,)()()()(12121njixxxxxxxxrmkjjkmkiikmkjjkiikij得到模糊相似矩阵