DNA序列的分类模型 汤诗杰周亮王晓玲 (中国科技大学,合肥230026 指导老师孙广中 编者按本文提出了DNA序列分类的三种模型,其一,基于A、G、TC四种碱基出现的频 率;其二利用了同一碱基在序列中的间隔,这一信息是单纯考虑频率所不能包含的;在第三种模 型中,作者把DNA序列视为一个信息流,考虑每增加一个字符所带来的信息增量.尽管文中倍 息量的定义方式仍可讨论,但本文思想新颖活跃,有其独特之处.本文最后的分类方法,是以上 三种的综合使用 摘要本文针对DNA序列分类这个实际问题,提出了相应的数学模型,为了很好的体 现DNA序列的局部性和全局性的特征,我们给出了衡量分类方法优劣的标准,即在满足一定限 制条件的情况下,是否能充分反映序列的各方面特性 依据我们提出的判别标准,单一标准的分类是无法满足要求的.我们的方法是侧重点不 的三种方法的综合集成、这三种方法分别体现了序列中元素出现的概率,序列中元素出现的周 期性,序列所带有的信息含量.利用这个方法,完成了对未知类型的人工序列及自然序列的分 类工作,最后,对分类模型的优缺点进行了分析,并就模型的推广作了讨论 1问题的提出(路) 2问题的分析 这是一个比较典型的分类问题,为了表述的严格和方便,我们用数学的方法来重述这个 问题.已知字母序列S1,S2,S3……S40,S1=x1x2x3…xn,其中x∈|a,t,c,gl;有字符 序列集合A,B,满足A∩B=g,并当1≤≤10时,S∈A;当1120时,S;∈B.现要 求考虑当21≤i≤40时,S1与集合A及集合B的关系 在这里,问题的关键就是要从已知的分好类的20个字母序列中提取用于分类的特征 知道了这些特征,我们就可以比较容易的对那些未标明类型的序列进行分类.下面我们将 首先对用于分类的标准问题进行必要的讨论 3分类的标准及评价 首先,我们提取的特征应该满足以下两个条件: (1)所取特征必须可以标志A组和B组.也就是说,我们利用这些特征应该可以很好 的区分已经标示分类的20个序列.这是比较显然的一个理由 (2)所取特征必须是有一定的实际意义的.这一点是决不能被忽视的.比如,如果不考 虑模型的实际意义,我们就可以以序列的开头字母为分类标准:已知在B类中的十个序列 都是以gt开始的,而已知在A类中10个序列没有以gt开始的,甚至以g开始的都没有 显然这是满足上面的第一个条件的.如果仅因此就认为这种特征是主要的,并简单的利用 这个特征将所有待分类的序列分成两类,显然是不甚合理的
DNA序列的分类模型 499 A 对于这样的一个复杂的分类问题,需要考虑的因素很多,也就是说,可供我们使用的分 类特征有许多如何从众多的因素中提取分类的主要因素,是我们处理这个问题的困难之 处.上面的第一个条件是我们的分类方法所必须满足的,可以看作是个限制条件;而第二个 条件是我们在设计分类方法时必须考虑到的,可以看作是对分类方法优劣的一种衡量,是 某种意义下的目标函数 4模型的建立及分析 由上面的分析可知,由于DNA序列本身的复杂性,我们很难在不知道确切的分类标准 的情况下,使用单一的方法来处理这个分类问题,由于,DNA序列同时具有局部性和全局 性的特征,我们尝试综合使用几种设计思想不同的方法来处理这个问题,以使该分类方法具 有好的分类性能和相当的健壮性 下面我们先从不同的角度出发,提出三种侧重点不同的分类方法,第一种从频率角度出 发,第二种从字母出现的周期性的角度出发,第三种从序列所带的某方面的信息量出发,并 给出它们单独使用时的分类结果.我们认为,这三方面综合考虑,可以较好的体现出序列各 个方面的特征,最后,从这三种方法出发,得到一个综合系统的分类方法,并利用它得到了最 终的182个序列的分类结果 方法1基于字母出现频率 不同段的DNA中,每个碱基出现的概率并不相同,从生物理论中,我们知道,编码蛋白 质的DNA中G、C含量偏高,而非编码蛋白质的DNA中A、T含量偏高.因此,A、G、TC的 频率中会含有很多的信息,下面给出A、B组的频率统计,见表1,表2(略) 由统计的数字可以看出,A组的碱基构成与B组的碱基构成有较大的不同A组的G 含量较高,B组的T含量较高.为做定量化的分析,引入数学中的内积概念,即将A、T、G、C 的频率分别作为四维向量的四个分量(PA,PG,Px,PC),现在我们得到两组向量A、B(i 1,2,3…10),然后将未知的序列21-40作为一个新的向量C,要将它归入A组或B组 我们可以尝试在 Hilbert空间中将向量归一化后求C与A组和B组的平均距离.记C、A、 B为归一化后的向量.为此我们计算内积和∑C·A与∑C,B,其中内积定义为欧氏 度量引导出的内积(c1,c2,3,C4).(a1,a2,a3,a4)=c1a1+c242+c3a3+c4a4,即 内积=(PA,P,Pr,Pn)A(PA,Pa,Px,P)短角 内积小的两个序列,我们可以认为它们的相关性小,而内积大的序列,我们就认为其相 关性大因此如果∑CA>∑CB,则认为C应归入A类,否则认为它应归入B类 计算结果如表3所示 由此,我们找到了区分C组的一种方法,这种比较∑C·A和∑C·B1的方法,我们 可以归纳为一个目标函数F1(1),即 F1(4)= ∑C·B
500 全国大学生数学建模竞赛优秀论文汇编 表3 二未知的与A组的与B组的属于的未知的与A组的与B组的属于的 内积 内积 类别 内积 类别一 0.938814 0.920957 2 0.9269220.803952 12 0.8669760,853967 A 0.9397270.656827 130.8609550.917122 0.7885240.937135 0.9616890.67678 456789 0.9481940.772073 150.9603220.73909A 0.8012010.930121 16 0.90428 0.747578 0.9530190.76695 A 170,9447240.723664 0.7460710.968035 180 0.9546521B 0.9310070.613193 19 0.8856310.811837 A 0.8977740.844082 0.75584 0.941 方法一讨论这种方法是从概率统计的角度分析问题,通过对每个字母出现频率的计 算,找出A,B两类DNA链中的频率特性,建立四维向量空间,然后对待求分类的序列统计 频率,与已知分类的向量进行内积运算,找出量化的关联性,从而将其分类,但这种方法也 有其局限性,在统计字母出现的频率时,忽略了字母所在位置以及各个字母之间的相互关 系,造成用这种方法对已知分类的序列进行检验时,个别频率特性不明显的序列不太容易分 类.所以,这种方法虽然有其科学性,但还不够完善,不能完全体现序列的所有特征 方法二甚于字母出现周期性 在以上进行了基于字母出现频率的分类之后,我们认为,一个序列所含的信息远不止每 个字母出现的频率,还有字母出现和它前后若干个字母的相关联性,字母在序列中出现的规 律性等等.前一个问题我们留到下面讨论,现在我们想办法处理后一个问题 对于某单个字母,以a为例,假设它在序列中第t1,t2,……,t+1,个位置出现,我们试 图找出这些数字之间的关联.首先,可以认识到考查t1的分布及绝对值是意义不大的,因为 序列是一大段DNA中的一个片断,片断的起始段不同会导致t1的不同.于是为了抵消t 的线性位移,考虑下面一组值 即字母a出现的间距 可以看出,序列51,s2,"…sn的大小包含的信息是a的“稠密度”,也可看成一个与频 率有关的量,前面已经处理过,所以我们可以考虑序列s1,2…,的波动幅度,幅度越 小,说明s(i=1,2,……,k)的值越趋于统一,即a的出现周期性越大而表征波动幅度的 量在统计中是中心矩.现求s的二阶中心矩,即方差 Var (s n-12(:-3)2,s= 同理,可以求出var、Var,、Va 由所得数据知,对varx与Var,上述方法对A、B组的区分率很高,就有良好的可分辨 性,为了强调这种待征的显著性,我们用F2=Var/Var作为这种方法的目标函数 由图1可以看出点与原点连线的斜率在A组中和B组中有显著差别,根据这个特征,A
DNA序列的分类模型 50 组和B组可以很好地区分开来,并且较好地弥补了方法一中的不全面之处 车个1, 出的 g间隔方差 方法二讨论这种方法是从序列中相邻相同字母之间的距离即字母出现的周期性着手 分析的,它统计了每个字母在序列中两次出现的间隔,并且用方差度量这种间隔的波动大 小,由此找到了一个能较好区分A,B组的目标函数,综合地考虑了序列全局和局部的性质 方法三基于序列熵值 我们可以把一串DNA序列看成一个信息流,这与生物学的基础知识是相应的.关于 A、B的分类,可以考虑其单位序列所含信息量(即熵)的多少,从直观上来看,我们可以认 为,重复得越多,信息量越少.,这是我们通过观察A、B组的特点而归纳出的方法 设序列为L=(a1,a2,a3,…,an);前m个字符所带的信息量为/m(1),记 :gm(1)=fm(1)-fm-1(1), 即gm(1)为加上第m个字母之后所增加的信息量.然后,由gm(1)=f=(1)-fm1(1),得 fn、1)(1),则f(D)为整个序列所带的信息量、F3(1)2即为单位长度所带 的信息量.现在的问题就归结为如何找出一个合适的gm(4) 我们有理由认为:g具有以下性质 性质1:gm(1)>0.即任意加上一个字符,它或多或少带有一定信息量; 性质2:第m个字符(或者是以它结尾的较短序列)与前面的序列(信息流)重复得越 多,gn()的值必然越小 性质3:第m个字符(或者是以它结尾的较短序列)如果和与它靠得越近的重复 (1)的值越小;和与它离得越远的重复,g=(1)的值越大;( 性质4:f0(1)=0 对此,我们可以构造如下函数 gn (L) b+t11+r202
全国 大学生 数学建模竞赛优秀论文汇编 其中b为防止分母为零而设的一个小正数;与 o,=2aaa; n 1以第m-t个字符结尾的字串且与以第t个字符结尾的i字串完全相同 否则 a为一个小于1的数,其存在体现了g的性质3.即如果越近的位置出现重复,认为字 串信息量越少,反之较多 G;的表达式中,t表示两个相同字串之间的距离,i表示字串长度,这个表达式定量的 给出距离和信息量之间的关系 又由于长度不同的字串重复对信息量的影响是不同的,所以必须在a前乘上一个权值 ,由概率统计的知识可知,这种影响是呈指数上升的,则可选择一适当的常数c>1,使得t =c-1,这个表达式定量的给出长度和信息量之间的关系 可以认为,字串长度太大的重复非常少见,则可将p取为某一固定的正数.那么,给出 a、b、c、四个参数,就可以把f严格确定下来.通过反复上机搜索,我们认为,取p=6,即 只检查长度为1到6的字串即可3可 另外,取a=0.392,b=0.1,c=3可以将AB组F3(1)值分得较开,并可以用来处理 未知数据 方法三讨论这种方法从序列的信息量(熵)入手,认为当序列中有大量的重复元素时 信息量就会比重复少的序列所含有的信息少.所以,其侧重点是是序列前后的重复性,也就 是序列元素的相关性.从所给的A,B两类中可以很清楚地看到B中序列重复量大,所含的 信息明显少于A组,而这个特征就被我们定义的熵函数凸显出来.将DNA序列看成一个信 流的方法由于其在实际问题中的广泛背景,将会是一个很有价值的想法,统计学和信息论 的一套非常成熟的强大工具也会在DNA研究中发挥巨大的作用.其 综合模型的建立 以上我们分别用三种方法得出了分类方案,这三种方案分别基于三种不同的方面对问 题进行分析.第一种方法主要考虑的是单个字母出现的频率;第二种方法主要考虑每个字 母的出现是否具有周期性;而第三种方法则考虑的是每条DNA所蕴含的信息量.我们将这 令人满意的结 考虑序列某一方面的特征,所以,总有一些不尽如人意的地方,于是,我们认为应该把三种 方法综合起来考虑,使序列各方面的特征都能得到体现,以使分类更加科学 下苗就是我们将几种方法综合考虑得到最后结果 以上我们用三种方法得到了三个目标函数:F1(1),F2(1),F3(1),这三个目标函数可以作 为分类的判别标准.将它们看成定义在序列空间L=|是由a,g,t,c四个字母组成的序 列!上,作用于实轴上的函数,现在,我们必须找到一个函数F,使得F可以体现序列的各个特 征 由于F1(1),F2(D),F3(1)的值域范围差别很大,为了有效的比较这三个函数,我们必 须将它们归一化,将2=f(1)(i=1,2,3,以下同)看成一定义在L空间上的随机变量,A 为L的子集,则将f归一化得