（信号与信息处理专业论文）基于信号处理方法的基因识别算法研究.pdf

中文摘要 进入2 1 世纪，随着计算机技术的不断进步以及其在生物信息领域的应用， d n a 测序工作正接近尾声。面对庞大的基因组数据，对数据进行解读，完成基 因序列不同部分的解读成为基因组分析中至关重要的部分。本论文通过将信号分 析领域的算法与基因识别相结合，使用信号领域的方法来发现生物基因组中新的 特性，并利用该特性实现基因识别。 论文介绍了生物信息学的基本理论、发展和研究的主要内容。通过分析原核 生物和真核生物基因组的结构特征，得出终止密码子的对编码蛋白质的限制作用 以及编码区的周期三特性。 在此特性基础上，介绍了将d n a 序列转化为数值序列的常用映射方法，并 提出一种降维的映射方法。然后利用现有的数字信号处理方法如短时傅立叶变 换、g a b o r 变换及小波变换等分析不同物种的基因组，通过观察编码区的周期三 特性进行基因识别，达到了预测编码区的目的。通过对大量d n a 序列进行分析， 比较不同方法的预测效果。并做出相应的改进。 论文的最后一部分通过将d n a 序列的基因识别过程转化成对d n a 序列的 分类，利用统计学原理中的支持向量机的学习方法结合d n a 序列中启动子的结 构特征，使用已有的数据库生成训练集，实现对未知序列的分类。 关键词；基因识别；编码区；谱分析；支持向量机；序列分析 a b s t r a c t c u r r e n t l y , w i t ht h ed e v e l o p m e n to fc o m p u t e rt e c h n o l o g ya n di t sa p p l i c a t i o n si n b i o i n f o r m a t i c s ，d n as e q u e n c i n gw o r ki sn e a r i n gc o m p l e t i o n f o rt h el a r g eg e n o m i c d a t a ，i d e n t i f i c a t i o no fd i f f e r e n tp a r t so ft h es e q u e n c e si so n eo ft h em o s ti m p o r t a n t t a s k si nb i o i n f o r m a t i c s t h i sd i s s e r t a t i o ni sc o m m i t t e dt oc o m b i n et h es i g n a la n a l y s i s a n dg e n ei d e n t i f i c a t i o nm e t h o dt od i s c o v e rt h ed i f f e r e n tc h a r a c t e r i s t i c so fd i f f e r e n t s t r u c t u r e si ng e n es e q u e n c e s ，a n dt h e nw ec a nu s es o p h i s t i c a t e ds i g n a lp r o c e s s i n g m e t h o d sf o rt h i sf l j e a $ f i r s t , t h i sd i s s e r t a t i o ni n t r o d u c e st h eb a s i ct h e o r ya n dd e v e l o p m e n to f b i o i n f o r m a t i c s a n dt h e nb a s e do nt h eg e n o m es t r u c t u r eo ft h ep r o k a r y o t e sa n d e u k a r y o t e s ，i ti sc o n c l u d e dt h a tt h ec o d i n gs e q u e n c e sh a v et h ef e a t u r eo f3 - p e r i o d , w h i c hc a nb er a r e l yo b s e r v e di nt h en o n - c o d i n gs e q u e n c e s s e c o n d ，b yt h ed i f f e r e n c eo fc o d i n ga n dn o n - c o d i n g ，w ei n t r o d u c es e v e r a l m e t h o d st or e p r e s e n tt h ed n as e q u e n c e si n t od i g i t a ls i g n a la n dp r o p o s ean e w m a p p i n gw a y t or e d u c et h ed i m e n s i o no ft h en u m e r i c a ls e q u e n c e t h e nw ec a nu s et h e m e t h o d so fs i g n a lp r o c e s s i n gt oa n a l y z ei t , s u c ha ss t a t i s t i c a lc o r r e l m i o n sa n a l y s i s ， f o u r i e rt r a n s f o r m ，w a v e l e tt r a n s f o r ma n dd i g i t a lf i l t e r b yo b s e r v i n gt h e3 - p e r i o d f e a t u r eo ft h es e q u e n c e ，w ec a ni d e n t i f yt h ec o d i n gr e g i o n s f i n a l l y ，w et u r nt h eg e n ei d e n t i f i c a t i o np r o c e s si n t ot h ec l a s s i f i c a t i o no ft h ed n a s e q u e n c e w ec a nu s et h es t a t i s t i c a lt h e o r yo fs u p p o r tv e c t o rm a c h i n e t os u m m a r i z e t h ec h a r a c t e r i s t i c so f d i f f e r e n tp a r t so fd n a s e q u e n c ei nt h ee x i s t i n gd a t a b a s e t h e n t h es v mc a na c h i e v et h ec l a s s i f i c a t i o no fu n k n o w ns e q u e n c e k e yw o r d s ：g e n ei d e n t i f i c a t i o n ；c o d i n gr e g i o n ；s p e c t r u ma n a l y s i s ；s u p p o r t v e c t o rm a c h i n e ；s e q u e n c ea n a l y s i s 第一章绪论 第一章绪论 继达尔文进化论的发表、遗传定律以及d n a 双螺旋结构的发现，使生物学 进入了分子生物的时代。如今，随着人体基因组测序工作完成，专家已经着手开 始破译基因组中的全部遗传信息。 本章通过介绍生物信息学的研究内容、发展情况，对原核生物和真核生物的 基因组特性以及基因识别一些主要算法进行分析，进而阐述了本文研究的目的、 意义以及研究的内容和方法。 1 1 生物信息学及其发展现状 由于目前针对生物体的计算机测序技术的不断改进，不同生物的基因组测序 工作逐步完成，测序后的数据都被存放在基因数据库中，这也使基因数据库中的 数据变得越来越庞大。而我们对d n a 序列的了解却显得十分匮乏。一边是大量 的基因数据以及其背后隐藏着人类目前尚不知道的生物学知识：另一边是我们渴 望利用这些新的发现及新的生物特性帮助人们来了解自身和周围的生物环境。这 是如今科学家们所面临的一个严峻的挑战。“生物信息学”正是为了迎接这种挑 战而发展起来的一门新兴学科。 生物信息学是2 l 世纪的核心学科之一，它是由数学、计算机科学、生物学 等不同学科交叉而成的，如今蛋白质序列以及基因序列的研究成为它的两个主要 研究方面，即通过分析蛋白质和核苷酸序列，使人们了解序列中所包含的遗传信 息。 1 2 生物信息学主要研究内容 基因组信息学，主要是指通过不同的技术手段将基因中的信息提取出来并做 相应的解读、储存、处理等，是生物信息学的核心部分。其首要前提是“解读 基因组的所包含的信息，即d n a 序列中不同部分在生物体中所起的作用。目前， 生物信息学的可以分为以下几个研究方向： ( 1 ) d n a 序列比对。其处理思想是将测试序列与数据库中的序列进行比较， 通过遍历数据库，找出与之有相似结构的序列。由于要读取d n a 序列的各个部 分，而且需要搜索庞大的基因组数据库，所以当序列长度较小且结构简单时，使 第一章绪论 用该方法可以很好的找到同源基因，而对于长的结构复杂的序列，往往需要耗费 很长的时间而且准确率也很低。 ( 2 ) 蛋白质空间结构及生物功能预测及序列比对。与序列比对的原理相同， 只是研究对象从d n a 序列的角度上升到蛋白质链的角度。通过搜索数据库，找 出在空间结构以及生物功能上与之相似的蛋白质链。由于蛋白质是负责生物体机 能的，所以在生物体中，空间结构相似的蛋白质链，其表现的生物功能往往相同。 ( 3 ) d n a 序列信息的获取，解读及分析。通过对整个d n a 序列分析，实 现对序列中不同区域的识别和鉴定，包括对编码区、非编码区、启动子、终止子 等结构特性的分析。并在此基础上得出与之相符的理论模型或对d n a 序列进行 标记和功能注释。 ( 4 ) 生物信息学技术方法研究。因为生物信息学研究的对象从最小的分子 结构到大的生物体，其中包含的信息多种多样，所以对它研究的过程不仅仅要利 用其他领域的研究方法如物理学、数学、化学等。还需要像神经网络分析、聚类 分析、机器学习等在计算机领域得到广泛应用的分析方法。 1 3 研究的意义及目的 本论文的研究内容主要放在生物信息学中的基因识别部分，由于基因库中的 数据繁多，且不同物种间的基因序列差异很大，基因识别的要解决的问题就是面 对一个未知生物的核苷酸序列，通过统计、变换等方法识别出该序列中是否存在 能编码蛋白质的基因以及该基因区在序列中的位置。由于真核生物的d n a 序列 中存在大量非编码区，在转录表达编码蛋白质的过程中该编码区会被丢弃。例如， 人体中能够编码蛋白质的基因区仅是基因组总数的3 - 5 1 2 1 。所以，面对如此庞 大的数据，选择合适的识别算法尤为重要。 序列比对是目前对未知d n a 序列进行处理分析的基本方法，通过搜索数据 库，找到与待测序列结构相似的序列，并认为该序列与待测序列有相似的遗传信 息，c l u s t a l 、f a s t a 、b l a s t 等1 3 】是序列比对常用的软件。但由于d n a 序列一般都 是很长的序列，而且基因库的数据繁多，所以使用比对的方法来分析速度比较慢， 其次由于数据库之间的差异，目前很多生物的基因组信息没有被录入，所以不同 数据库进行序列比对的结果也各不相同。 由于d n a 序列可以被划分为不同的区域，所以针对d n a 序列的基因识别 可以分为启动子、终止子的识别和编码区的识别等。由于d n a 序列是字符序列， 这很大程度上使我们无法使用其它领域的研究方法来对它进行分析。所以如果通 过一定的映射规则将它映射为数值序列，这样一条d n a 序列就可以看作一个离 第一章绪论 散的时间信号，而d n a 序列的信息也全部包含其中。在此基础上我们就可以使 用信号领域的一些分析方法，如自相关分析、傅立叶变换、小波变换等来对离散 时间信号进行分析，观测d n a 序列潜在的周期特性。尽管d n a 序列的谱分析 取得了很好的成果，但还没有形成完整的理论。如果能将数字信号处理中的方法 很好地运用到d n a 序列分析中，提高其在d n a 分析的准确度，那么凭借已经 成熟的信号处理一些快速变换方法，一定能够大大提高物种基因分析的速度和准 确性。 1 4 本文的主要工作 本论文可以分为三个部分。 第一部分主要对原核生物和真核生物的基因组的基本结构，以及编码区的核 苷酸分布特性和周期性进行分析，这也是进行基因识别的理论基础。 第二部分介绍了将d n a 序列映射为数值序列的一些方法，分为实数域映射 方法和复数域映射方法。并介绍了信号处理中的一些谱分析方法，主要是短时傅 立叶变换、g a b o r 变换和小波变换。利用这些方法进行序列的基因识别，比较不 同算法的优点和存在的问题。 第三部分介绍了统计学理论中的支持向量机的学习方法，并将它应用到基因 识别中，通过分析编码区启动子的结构特性，利用这个特性实现对d n a 序列进 行分类。 第_ 二章基冈识别的基本原理 第二章基因识别的基本原理 本章首先对原核生物和真核生物的基因组结构以及已有的基因识别算法进 行了介绍，并进一步分析d n a 序列中编码区的核苷酸的分布特征，从中提取出 基因区的周期特征。 2 1 原核生物及真核生物基因组 不管是原核生物还是真核生物，它们全部的遗传信息都包含在基因组中，而 这些基因组基本都是由脱氧核糖核酸也就是d n a 构成( 只有特殊的一小部分病 毒的基因组是由r n a 组成) ，所以要想分析了解生物的遗传信息，首先应对这些 生物的d n a 序列的结构进行解读。 图2 1 给出了原核生物的基因结构，图中的“编码区域 就是基因识别问题 中需要寻找的能够编码蛋白质的区域，它是由起始密码子( a t g , g t g , t t g ) 和 终止密码子( t g a ，t a a ，t a g ) 围起的一段d n a 序列。由于一个蛋白质由三联 密码子完成编码，所以一段完整的编码区长度都能被三整除，且不会被内含子隔 断。而且原核生物的编码区编码的蛋白质长度一般都大于6 0 个氨基酸。这样可 以总结出原核生物基因组的结构特征：( 1 ) 基因组结构简单，只存在少数的间隔， 与非编码区相比编码区域较大，一般能达到8 5 左右【2 】。( 2 ) 几乎所有原核生物 的编码区域都是完整的一段编码区，不会出现真核生物中外显子和内含子交错出 现的情况。 5 ， 启动子 终止区 3 ， 图2 1 原核生物的基因结构 图2 2 和图2 3 给出了真核生物的基因组序列以及编码区的结构，其序列由 不同的部分组成。由于真核生物基本都属于多细胞生物，所以虽然同一生物体中 不同功能的生物细胞都有该生物的全部遗传信息，但由于不同的细胞发挥的功能 不同，所以不同细胞内的基冈表达方式也不相同。真核生物体内的基因组在长度 4 第二章基冈识别的基本原理 方面受到的限制很少，存在大量的“垃圾d n a9 9 0 正是由于这些特点，真核生物 基因组的转录、翻译和表达等过程都要比原核生物复杂得多。 在真核生物的基因组中。基因序列被划分为不同的区域，其中不能编码蛋白 质的序列为内含子( i n t r o n ) ，能编码蛋白质的序列为外显子( e x o n ) 。在真核生 物体基因组序列中，内含子要比外显子多得多，上一章就讲到人体基因组中仅有 3 - 5 属于编码区，而剩余部分均为内含子即非编码区【2 】。且基因区也不是连续 存在的，而是被内含子打散分布在不同的区域，这一点也是对真核生物基因识别 的难点所在。 从图2 3 中还可以看出真核生物基因组基因区的核苷酸分布特征，大多数真 核生物的内含子都遵照g t - a g 规则，之所以这样命名原因是每个编码区和非编 码区的连接部分都有一段特殊结构的序列，即非编码区由g t - a g 围成。 上游增强子t a t a 盒 外显子 内含子 n 终止区 起 图2 - 2 真核生物基因结构示意图 图2 3 外显子和内含子交错存在 2 2 原核及真核生物基因识别算法 由于基因数据的庞大以及基因结构的复杂性，计算机依靠其速度的优势就成 为完成基因识别的最有力的工具。基因识别工作最初只是将重点放在编码区，而 不考虑非编码区的影响，随着生物医学的发展，以及对生物体的认识逐渐加深， 基因识别的工作开始对完整的基因进行分析。如今已经发展出数十种计算机识别 算法来预测蛋白质编码基因。 5 第二章基因识别的基本原理 由于真核生物基因结构的复杂性，实现其基因识别的过程要比原核生物复杂 很多。因为真核生物的基因识别包括了对不同结构的功能序列的识别，如对外显 予内含予的识别、启动子和终止子的识别、5 和3 剪接位点的识别，而有些生物 往往内含子非常多，很难实现对各个部分的正确的识别。目前比较流行的针对真 核生物基因识别的软件主要是基于隐马尔可夫模型、神经网络、决策树、曲线判 别等理论。 虽然这些识别算法在一定程度上能完成基因识别的工作，但由于理论模型的 不完善或侧重点不同，这些算法都或多或少存在一些缺陷。如基于同源性的识别 算法由于受到数据库的局限，很多新的基因无法在数据库中找到与之相似的同源 基因，在这种情况下就会导致基因的识别的准确率很低。而对于真核生物的基因 识别算法，由于真核生物基因组复杂性的原因，其识别效果比原核生物的识别效 果还要差很多。往往在整条基因序列上对外显子的识别成功率还不足3 0 。基于 这种情况，就要求我们能够挖掘出更多编码区和非编码区的特征区别，希望能利 用在别的领域如数学、统计学、信号处理领域的已经成熟的方法来完成基因识别 的工作，来提高基因识别的准确率和计算效率。 2 3 蛋白质编码区的结构特性 本节将介绍在物种的基因组中终止密码子中的碱基分布特性以及其对限制 蛋白质编码区长度等方面的重要作用，还介绍了三联密码子中的碱基分布的周期 特性。蛋白质编码区的这些结构特性，也为后期对蛋白质编码区的预测奠定了理 论基础。 2 3 1 终止密码子碱基的分布特性 这一节我们将重点放在作为o r f 结束点的终止密码子上。因为作为蛋白质 编码区的结束位置，为何终止密码子会对蛋白质的形成起到抑制作用，我们可以 通过分析核苷酸的分布特性来得到答案。因为虽然物种之间的基因组结构特性各 不相同，但经过对大量物种基因组进行统计发现，对于三联密码子，不同的核苷 酸在密码子中不同位置上的分布具有一定的倾向性。如核苷酸t 、g 、a 分别在 三联子的第一、二和第三位置上出现的概率要比在其他位置上出现的概率小很多 0 7 。而且这个特性在大多数物种的基因区都存在。而我们又发现，蛋白质编码区 的终止子t g a 正好是由这三种核苷酸组成，将这两点联系起来，就可以说明为 何终止子会限制d n a 序列进行蛋白质的合成。在达尔文的物种进化论中提到， 不同物种经过多年的演变，它们的基因组结构往往会受到外部环境、自身突变等 6 第二章基因识别的基本原理 多方面不利因素，而蛋白质在生物体中又扮演着至关重要的位置，生物体需要它 来完成不同的功能，那么这在客观上就要求蛋白质的空间结构及编码蛋白质的序 列具有一定的稳定性，以此来抵抗内部选择以及突变对其功能造成的改变。 通过上述讨论，我们可以得出终止密码子最大的作用就是限制蛋白质的形 成。当然，t a g ，t 从这两个密码子也属于终止密码子，但使用频率远低于t g a 的使用频率。因为经过对能编码蛋白质的基因区统计分析，发现在密码子的二、 三碱基位出现a g 或从的含量比出现g a 的含量要高，这也说明t a g 、t 从 密码子对蛋白质的形成所起的限制作用要小于密码子t g a 。 下图给出了碱基t 、g 、a 在密码子不同位置上的分布倾向性。 眦鼍r l i 图2 - 4 不同碱基在密码子不同位置的分布 正是由于终止密码子在单核苷酸和双核苷酸水平上有很多分布特性，那么我 们在基因识别的过程中，可以利用它的这些特点进行识别，也可以对d n a 序列 进行分类。 2 3 2 编码区核苷酸的周期性 本节的目的是发现不同的核苷酸在三联密码子中的位置的非对称性。通过不 同的统计参数，来观察编码区核苷酸的周期特性。 位置的非对称性总和 通过对已知的编码区和非编码区进行测量，该估计值是一个线性加权的组 合，并且计算出不同的核苷酸在不同位置的频率特性。用，( 6 ，) 来表示d n a 序 第二章基冈识别的基本原理 列s 中，核苷酸b 在相应的三联之密码中的位置r 的频率。则 六( 6 ) = 二。l ( b ，r ) 3 表示核苷酸b 在三联密码子中不同位置的平均频率。然后 计算此平均值的方差，来表示核苷酸b 位置的非对称性，即 a s y m ( b ) = 二。( 工( 6 ，) 以( 6 ) ) 2 ，该值是独立于相应的测量序列的然后进行做 和，计算出不同核苷酸位置的非对称性的总和p a ( s ) ： p a ( s ) = a s y m ( a ) + a s y m ( t ) + a s y m ( g ) + a s y m ( c ) ( 2 - 1 ) 通过对人体1 7 6 1 个外显子和1 7 5 3 个内含子进行统计分析，计算出在外显子 和内含子区核苷酸的位置的非对称性指数p a ( s ) 。图2 5 给出了针对测量序列的 图2 5 外显子和内含子中核苷酸的p a 指数 对不同核苷酸在编码区的位置的周期性进行测量，如果给定了一个d n a 序 列，我们就可以计算出在核苷酸f 后相距k 个位置出现核苷酸的次数，( 七) ，例 如，给定的序列为 s a g g a c g g g 峨c 入 则按照上面的定义可以计算出：( 1 ) = 2 ，n r ( 0 ) = 1 ，g g ( 0 ) = 3 ， 。( 7 ) = 2 等。图2 - 6 给出了针对人体的一段编码区序列进行测量，在相距k 位 置上出现a a 的次数： 8 第二章基因识别的基本原理 夕 ： 对 出 现 、犬 数 a 0 湖a o 】o n s m td n s 攻 * 、 一、，、一n 一、- ，。， l ， 、 ， 、， 。- ，- 、一 、“ 。、- o g b 甘 ，j 1l u uvu 曰l l 斗 一 - 、 f 鬯 j u 1 uu l 2561 11 41 7却为a3 23 53 64444 7 图2 - 6 人体外显子和内含子a a 对出现次数 通过上图可以看出，在外显子序列区域存在明显的周期特征。在核苷酸a 后的位置k = 2 ，5 ，8 出现核苷酸a 的概率远大于其他位置。并且注意到在该位置 上核苷酸a 都是在相同的三联密码子中的相同位置。这样的周期规律，也反映 出了在编码区的核苷酸的位置关系。 周期不对称指数( n u ) 给定一个d n a 序列s ，考虑三种不同的概率，概率p m 表示在距离k = 2 , 5 ，8 找到相同的核苷酸对的概率，概率兄，表示在位置七= o ，3 ，6 找到相同核苷酸对 的概率，概率艺，表示在位置k = 1 , 4 ，7 找到相同核苷酸对的概率。基于编码区 的周期三特性，可以推出在编码区只。要大于心和p 三，而在非编码区这三个概 率应该相近，则核苷酸这种周期三的非对称性可以用周期不对称指数( p a l ) 来 表示 p a l ( s ) ：竺堕矧i2(2-2) m i n ( p 加，) 我们就可以利用此参数来对一个d n a 序列s 的编码区进行预测，结果如图2 - 7 。 9 第二章基i 天1 识别的基本原理 一 一二二。 周朔不对称指数 图2 - 7 外显子和内含子中核苷酸的周期不对称指数( p a l ) 给定一个d n a 序列s ，用概率既( 七) 来表示在核苷酸i 后相距七个位置出现 核苷酸，的概率，该值可以根据上面提到的n 口( 七) 来计算得出。则去掉序列s 的 影响，来计算核苷酸i 后相距七个位置出现核苷酸的概率： p l = p t | u o p | pj(2-3- 式中a 和p ，表示核苷酸i 与核苷酸_ ，在序列s 中出现的概率。因此，针对不同的 距离七，1 6 种不同的核苷酸对的户信息概率关系可以通过上式得到。进而得到针 对不同七的户信息量。 m ) = “刚觚州圳0 9 ( 名 口4 ) 图2 - 8 给出了针对人体的1 7 6 1 外显子和1 7 5 3 内含子序列的互信息量计算结 果。 1 0 第二章基冈识别的基本原理 2 4 小结 平均互信息量 图2 8 外显子和内含子的平均互倍息量 本节首先介绍了原核生物以及真核生物的基因组的基本结构，并介绍了现有 的针对原核生物和真核生物基因组的基因识别算法。为了能和数字信号处理相结 合，并在基因识别的过程中使用信号处理中的一些时频分析方法，分析了生物基 因组中编码区的外显子、内含子以及终止密码子的碱基结构特性、分布特性和周 期特性。 第三章d n a 序列的映射及谱分析原理 第三章d n a 序列的映射及谱分析原理 上一章分析在编码区外显子中的核苷酸存在位置的非对称性以及周期三特 性，而在内含子中不具有这种特性。本章将介绍将d n a 字符序列映射为数值序 列的一些方法，包括实数域映射和复数域映射的方法，提出一个在不丢失频谱信 息的情况下通过系数矩阵将4 维子序列转化为3 维序列的方法。最后介绍了数字 信号领域常用的谱分析方法的基础理论和性质。 3 1d n a 序列的映射为数值序列 在用数学方法、物理方法和数字信号处理方法对d n a 序列进行分析的研究 中，首先需要把d n a 序列按照一定的规则映射成相应的数值序列。通过将该字 母序列数值化，数值化后的序列可以看作一条离散的时间信号，通过上一章我们 了解到d n a 序列的蛋白质编码区隐藏着很多周期性质，通过合适的映射关系， 这些周期性也保留到了离散的时间信号中。我们就可以使用信号处理中的方法对 该离散信号进行分析，挖掘信号的时频域特征，从而可以快速的对基因序列进行 周期性分析、基因识别和同源性等方面的分析。而且d n a 序列数值映射方法的 优劣会直接影响到最终分析结果以及生物学意义的解释。 3 1 1 常用的基因序列数值化方法 卟忙絮 t t 法 卟忙搿 g g 法 卟 s 翳 c 百法 卟 糍： 1 2 ( 3 - i ) ( 3 - 2 ) ( 3 - 3 ) ( 3 - 4 ) 第三章d n a 序列的映射及谱分析原理 h y d r o g e nb o n de n e r yr u l e ( s w 规则) ： 卟锰篇笏 嘌呤一嘧啶法( 1 w 法) ： 蜘船：搿 复数域映射法 互补表示法 ( 3 - 5 ) ( 3 - 6 ) ( 3 7 ) ( 3 - 8 ) 其中刀= 1 , 2 ，n 表示表示d n a 序列在位置刀的相应符号的映射值。这样 d n a 字母序列就映射成为了采样频率为1 h z 的时间离散序列。与采用实数域的 映射规则相比，采用复域的映射规则可以减少序列的直流分量，使交流分量性质 更加突出。 3 1 2d n a 的4 d 表示法 上一章在讨论终止密码子的结构特征的过程中，我们曾用双核苷酸的分布特 性进行研究。因此在映射过程中，如果我们只考虑单个碱基的映射，可能会把隐 含在双核苷酸甚至三核苷酸中的生物信息丢失，所以这里讨论一种基于双核苷酸 的序列映射方法【2 。具体的变换过程如下：首先我们依据化学特性的区别，将组 成d n a 序列的核苷酸a ，t ，g ，c 分进行分组：嘌呤r = a ，g 嘧啶y = c ，t ) ；氨基 m = a ，c 酮基k = g ，t ；弱氢键w = a ，t 强氢键s = c ，g ，。然后列出双核苷酸 所有的不同组合结果： c t ，g a ，a g ，t c ，a c ，c a ，g t ，t g ，a t ，t a ，c g ，g c ，a a ， t t , c c ，g g ) 。依据上面三种碱基的不同分组，进一步把所有双核苷酸组合按化 学特性分为三组：分别为 a g ，g a ，c t ，t c ； a c ，c a ，g t ，t o ；( a t ，t a ，c g ，g c ； 还有剩余的四种组合为一组 从，t t ，c c ，g g 。这样我们就可以将一条字符序列 通过四维空间中的一系列点描述出来，通过曲线来对d n a 序列进行分析。一个 序列通过以下四维坐标( 五，y ，z ，s ，) 来表示： 彳r g c = = = = m 帕功功“双“ ， ，l l 1 f 一 一 ，j、_l = 一 u 彳r g c = = = = 力力帕m“双娟“ ， ，l l - “叫卜 一 一 ，-fl-li_l【 = 一 v 第三章d n a 序列的映射及谱分析原理 彳g + g + c z + 丁c ， x 。= - 二_ _ _ _ _ - 二- 一l _ 二 以一l 彳c + c 彳。+ g t , + 刀g ， v = = 一 a t , + 7 ，+ c g ，+ g c z 。= 二l 以一l 4 彳+ c c ，+ g g 。+ 冠 胃 = 一 。 以一1 ( 3 - 9 ) 式中对于每个碱基在空间中坐标( x ，乃，z ，s ，) 可以用来表示不同双核苷酸在d n a 序列中的频率分量。 3 1 3 一种新的降维的方法 通过将d n a 字符序列映射为一个4 维的子学列，那么序列x - 】就可以表示 为： 研”】= ”月【刀】+ “7 【刀】+ 材c 【玎】+ “g 【刀】，刀= o ，l ，2 ，刀 ( 3 - 1 0 ) 在这种情况下，如果对子序列使用0 ，l 映射的方法，则存在： 甜【刀】+ “7 【疗】+ 球c 【刀】+ “g 【疗】= l ，扫，a l l 疗 ( 3 1 1 ) 说明四维的子序列是线性相关的，如果对子序列进行傅立叶变换观察频谱： 研七】_ 二研， e i n m , k = o ，l ，2 ，n - i ( 3 1 2 ) 通过以上两式可以推出： u c 七，+ u r c 七，+ u cc 七，+ u gc 七，2 主三： 。3 。3 ， 也是线性相关的，那么就可以通过减少一个参数来降低维数，节省计算量。通过 一组系数矩阵，与已知的四维子序列相乘，生成一个三维的子序列向量，具体处 理过程： t x ，x 窖 石。，= 三乏茎；差 r l 】 刀】 刀】 刀】 ( 3 - 1 4 ) 通过上述变换后，取一组变换系数：( a t ，d e , a b ) = ( 。，0 ，1 ) ，( t ，, t e , t b ) = ( 竽，o 一三) ， ( g , , g g , g b ) = ( 一等，一鱼3 ，一争，( ：( 一孚，孚，一争则可以将四维序列 1 4 第三章d n a 序列的映射及谱分析原理 化成。一个三维序列： 层 x ，= 半( 2 u r 【门卜“。，【刀】一z f g 胛】) - ， 广 x 。【刀】= 华( “f 【川一“g 玎】) 1 x 6 门】= ( 3 ”爿【聆卜d t 【胛】一u c 玎】一”g 【玎】) ， ( 3 15 ) 这样就可以把对四个序列的处理转换成对三个序列的处理，在序列很长的情况 下，既能保证不丢失频谱信息，又能很好地节省计算的开支。 为了实现对d n a 序列的频谱分析，对上述经过降维处理的三个序列进行傅 立叶变换得到x ， 七】，x g 【七】，x 。【足】，然后我们可以通过彩色图谱对该序列进行分 析，将三个频谱序列结合起来放在一张彩色图谱上进行观察。将x ，【玎 序列用红 色表示，k 行 用绿色表示，【”】用蓝色表示。图( 3 1 ) 是对一段线虫染色体 ( g e n b a n k a c c e s s i o nn u m b e r n c0 0 0 9 6 7 ) 的分析，使用d f t 窗口长度n = 6 0 。 该染色体有三段可编码蛋白质的基因区( 9 5 3 1 0 6 6 ，1 6 6 8 1 7 2 7 ，1 8 0 7 2 0 2 8 ) 。可 以看到在图中对应于基因区有明显的高亮部分。且对应的频率为1 0 ( 因为d n a 编码区有6 聚体的特征，k ：_ n ：1 0 ) 。 图3 1 线虫染色体基因中截取4 0 0 0 个核苷酸 第三章d n a 序列的映射及谱分析原理 3 2 常用的谱分析方法 通过将字符序列映射为数值序列，把它看作一条离散的时间信号，就可以使 用数字信号领域的序列分析方法来对它进行分析。如自相关分析、傅立叶变换、 短时傅立叶变换、小波变换等都可以用来观察信号的时频特性。 上一章中我们已经讨论过在基因中蛋白质编码区的外显子具有周期3 特性， 而内含子不具有该特性，这是由于在密码子转变成氨基酸过程中存在着编码倾向 性。所谓的周期3 特性【l 羽，即在序列的频谱在1 3 b p 处存在一个峰值。周期3 特 性在大部分蛋白质编码序列中存在。我们可以利用这个特性来对未知的d n a 序 列进行蛋白质编码区的预测。具体做法就是通过时频变换，来观察d n a 序列的 频谱或功率谱。 3 2 1 直接傅立叶变换 该方法是在基因序列数值化的基础上，直接对离散时间序列进行离散傅立叶 变换，求其功率谱。对于长度为n 的序列x ( 刀) ，其傅立叶变换为： 2 j r x 嘲= = x ( 玎) e 叫- 肼，k = o 1 ，2 ，n - i ( 3 一1 6 ) 其功率谱为 s ( 七) = l x 0 2 、n u ( 3 - 2 6 ) 1 8 第三章d n a 序列的映射及谱分析原理 | = 审竺砧l 仔2 7 ， 护2 赤脚列2d 引 。孑劢 一l 旌l 旌x ( a ) 窗口函数g f o ) 的图形( b ) a 取值不同时高斯函数的图形 图3 - 3 窗口函数与高斯函数的图形 利用上一节定义，计算g a b o r 变换的高斯窗的“时域持续宽度 和“频域持 续宽度”： 1 9 第三章d n a 序列的映射及谱分析原理 卜2 壶灿。2 i g ) | 2 j ( 3 - 2 9 ) 卜赤工面- - ) 2 嗣陬遗 通过式( 3 - 2 9 ) ，可以得出一个窗口在时频域对待测信号的分析面积为 i t + a - a g ，t + 口+ 4 】【万+ w o - a 雪，万+ 砜+ 4 】。时域宽度和频域宽度为2 a 窖 与2 以。那么通过平移该窗口( 图3 - 4 ) ，使窗口覆盖整个时频域，以此来分析 待测信号，就相当于将待测信号在时频域分解为一个一个的窗口来分析。这就类 似于使用放大镜来放大观察信号的每个部分的特性。 圈匪匿 ! ! ：l 圈匪 匪圈 圈圈圈圈1 匪匪 匪圈 图3 - 4 短时傅里叶变换分辨元胞的相空间表示 这种加窗的傅立叶变换就是将信号划分成小的部分来分析，但是由于每次变 换使用的是同一个窗口，而窗口的分析面积保持不变，这就使得当分析信号包含 很多高频、低频部分时，分析窗无法同时对高低频部分进行精确分析。什么形状 的窗函数才能最大限度得满足对非平稳信号的分析昵? 我们希望当窗口面对信 号的高频部分时，窗口的时域持续宽度很窄，而频域持续宽度很宽。当窗口面对 信号的低频部分时，窗口的时持续宽度很宽，而频域持续宽度很窄。图3 5 给出 了我们想要的理想分析窗口： 第三章d n a 序列的映射及谱分析原理 ，l。l “，ll 巳- il 巳。i l l气1 i l吒吒1 l j f。j f 77 3 2 4 小波变换 图3 - 5 理想的时频分辨率 上一节讲到理想的窗函数在时域和频域分析应该满足的形式，即当我们面对 的非平稳信号中包含很多高频和低频部分时，我们希望所使用的分析窗口可以随 着待测信号的频率变化调整大小，小波变换就是在这种要求下提出的。 小波分析理论在上世纪初就已经提出，然而直到8 0 年代之后，随着m a l l a t 将正交小波基的概念与通信领域、图像处理领域系统地结合起来，小波分析才进 入高速发展的阶段【5 】。虽然也是在傅立叶变换的基础上进行分析，但在变换的过 程中却能根据待测信号的不同做出相应的调整。小波变换与g a b o r 变换一样采用 加窗分析的思想，但它使用的窗口却可以根据所测信号的不同进行改变。 小波变换的具体实现是通过窗口虬。( f ) = 1 1 y ( 兰) 对信号进行变换 ia a w f ( a , b ) 2 击亡巾旷( 等) d ， ( 3 _ 3 0 ) 式中，6 为时间平移参数，通过它来确定小波函数的中心位置，a r ，a 0 ，a 为尺度参数，由它确定了小波函数的中心位置。正是由于引入这两个参数，才使 得小波函数的窗口实现可调。 当小波满足窗口条件和容许性条件时，可以得出小波的均值为零，且是具有 带通性质的窗口。通过对小波函数分析，在平移之前分析窗所占的矩形面积是 lf 竽，竽i ，而经过平移之后分析窗所占的矩形区域为： l 二 j la t 坐，竺士箜l 【 2口2 aj ( 3 3 1 ) 2 1 第三章d n a 序列的映射及谱分析原理 3 2 5 连续小波变换的计算 给定一个待测信号，我们如何使用小波变换来分析它呢? 从式( 3 3 0 ) 可以 得出，可以通过以下几个过程来实现连续小波的计算。如图3 - 6 ，首先选择一个 小波基，使其与待测信号进行积分，得出一个该时段的相关系数c 。该系数c 可 以理解为所该区间的待测信号与所选小波基在波形的相似评价。然后调整平移因 子b ( 图3 7 ) ，使小波基与不同时段的待测信号进行积分，直到该小波滑过整个 待测信号的时段。重新调整小波的伸缩因子a ( 图3 8 ) ，也就是改变小波的时频 分辨率，重复上面的过程。通过不断变换伸缩因子a ，就可以完成不同小波基下 待测信号的小波系数。 信号 小波 信号 小波 c = 0 0 1 0 2 图3 计算小波变换系数示意图 信号 小泼 信号与小没 c = o 2 2 7 图3 7 向右平移小波图3 - 8 调整参数丑重新计算小波变换系数 由小波变换的定义式( 3 3 3 ) ，有 町( 口，6 ) = ( ( f ) ，。( ，) ) = 二厂( ，) ，虻。( t ) d t = 二厂( r ) 下1 沙( t - b ) d t ( 口 o ，厂r ( 尺) ) (332)q 。 a a 、 其中，y ( r 净万1y ( 字) ，并设，( ，) = f ( k a t ) ，f ( k , k + 1 ) ，则 第三章d n a 序列的映射及谱分析原理 帅，6 ) = 军n ( f ) i 口一p 户 = 军拈l 口l - 1 2y , ( 等) d r ( 3 - 3 3 ) = i 口l 叫2 莓厂c 七，( = 1y ( 等) d r 一y ( 等) d ，) 那么小波变换可以通过以上5 步来实现，在m a t l a b 中，连续小波变换可 通过以下代码实现： m a t l a b 实现连续小波变换的代码 p r e c i s = 1 0 ；小波函数积分精度控制 s i g n a l = s i g n a l ( ：) ； i e n = l e n g t h ( s i g n a l ) ； c o e f s = z e r o s ( 1 e n g t h ( s c a l e s ) ，l e n ) ； n b s c a l e s = l e n g t h ( s c a l e s ) ； p s i _ i n t e g ，x v a l = i n t w a v e ( w n a m e ，p r e c i s ) ；计算从一到k 的小波积分序列 w t y p e = w a v e m n g r ( 。t y p e ，w n a m e ) ； i f w t y p e 一5 ，p s i i n t e g = c o n j ( p s i i n t e g ) ；e n d 判断是否为复小波，对复小波 取共轭 x v a l = x v a l - x v a l ( 1 ) ； d x = x v a l ( 2 ) ； i n d = l ； f o rk = l ：n b s c a l e s 循环计算各尺度的小波系数 a = s c a l e s ( k ) ； j = 【l + f l o o r ( o ：a x m a x ( a d x ) ) 】； i f l e n g t h ( j ) 一l ，j = 【ll 】；e n d f = f l i p l “p s i i n t e g ( j ) ) ； e o e f s ( i n d ，：) = s q r t ( a ) * w k e e p ( d i f f ( c o n y ( s i g n a l ，d ) ，i e n ) ；计算公式( 3 - 3 3 ) i n d = i n d + l ； 小波函数的功率谱为i 野( d ，b ) 1 2 ，则其平均功率为 僻盟b i - 坐b o 侮3 4 ，