FFT在生物信息学中的应用(带摘要)

一、绪言 自从1990年启动人类基因组计划以来，基因测序得到极大的发展，同时蛋白质结构信息也增长迅速，外加上其它的生物数据，这一切形成了一个巨大的生物数据海洋。然而数据并不是知识，与庞大的生物数据相比，相关的生物知识人类掌握的还很少，这便催生了一门新的学科----生物信息学。美国在人类基因组计划实施五年后的总结报告中，对生物信息学作了如下定义：生物信息学是一门交叉学科，它包含了生物信息的获取、处理、储存、分发、分析和解释等在内的所有方面，它综合运用数学、计算机科学和生物学的各种工具，来阐明和理解大量数据所包含的生物学意义[1]。

快速傅立叶变换（FFT）是1965年由J.W.Cooley和J.W.Tukey提出的一种旨在加快傅立叶变换（DFT）的快速算法，它将傅立叶变换的运算复杂度由2()On

降到了(log)Onn,因而广泛应用于各个科学领域。在生物信息学中，FFT也有广泛的应用，特别是在序列比对，蛋白质结构及功能预测等领域，经常能看到FFT的身影。本文旨在通过对一些应用实例的介绍，以期使读者对FFT在生物信息学中的应用有一个入门性的认识。

二、快速傅立叶变换(FFT)简介 （一）傅立叶变换（DFT） 设()xn为有限长信号，即 ()0,1,......,1()0xnnNxn



其它

(2.1)

则称 210()(),0,1,,1NjmnNnXmxnemN





 (2.2) 为()xn的离散傅立叶变换（DFT），并称()Xm为()xn的离散频谱。相应的逆变换为 2101()()0,1,,1NjmnNmxnXmenNN

 (2.3)

若引入记号2jNNWe，则有 1010()()1()()NmnNnNmnNmXmxnWxnXmWN













 (2.4)

（二）快速傅立叶变换（FFT） 考虑离散傅立叶变换的运算量，由其定义式(2.2)知需要2)1(N次复乘法运算，复杂度为2()ON，随着N的增大，运算量将飞速增加，这是令人难以接受的。1965年，美国的J.W.Cooley和J.W.Tukey提出了一种新的计算频谱的算法，将运算复杂度降到(log)ONN，使运算速度大大增加，此即所谓的快速傅立叶变换（FFT）。 FFT巧妙地运用了奇偶对分的思想。当N为偶数时，Cooley和Tukey发现，若记

12021202()(2)()(21)NlkNlNlkNlYkxlWZkxlW













 (2.5)

则有 ()()()0,1,,12()()(),0,1,,122kNkNNXkYkWZkkNNXkYkZkWk







 (2.6)

将式(2.6)用图表示如下，称之为蝶形图。 YkZkXkXN2kW

N

k

WN

k

图2. 1 FFT的蝶形 通过奇偶对分，将求一个N长序列的频谱化为了求两个/2N长的序列的频

谱，减少了运算量。如法炮制，若假定2KN，这样一而再、再而三地奇遇对分下去，直到长度为1时停止，如此以来运算量将会锐减。这便是所谓的快速傅立叶变换(FFT)，可以证明，这种算法的复杂度为(log)ONN。

（三）FFT加快相关系数的运算 设011(,,,)nxxxx,011(,,,)nyyyy为两个n长的数值序列，x和y的圆周相关系数定义为 10110(,,,)ncycnjkjkkzzzzxy

z (2.7)

其中下标kj表示的是模n后的值。易知计算两个等长序列的圆周相关系数的复杂度是2()On,而用FFT则可以降低其运算复杂度，加快运算。具体地，令

120(,,,)Rnnyyyy,X、RY、cycZ分别表示x、Ry、cycz经傅立叶变换后的频

谱，则有cycRZXY。即经过两次正的傅立叶变换和一次逆的傅立叶变换便可将相关系数计算出来。若在傅立叶变换中采用FFT快速算法，则可将复杂度降为(log)Onn。

不过，在生物信息学中，线性相关比圆周相关更常用。设有两个序列011(,,,)lxxxx和011(,,,)myyyy，令011(0,0,,0,,,,)lxxxx，

011(,,,,0,0,,0)myyyy，x、y均为长为1nlm序列，则x和y的线性

相关系数011(,,,)linnzzzz等于x和y的圆周相关系数cycz,于是线性相关系数的运算也可以通过用FFT将其复杂度降为(log)Onn。

三、FFT在序列比对方面的应用 （一）生物序列的表示 在生物信息学中，为简化生物序列的表示，用取自特定字母表的字符串来表示生物序列。以DNA序列和蛋白质序列为例，DNA序列由四种核苷酸组成，可用 “A”，“T”，“C”，“G” 四个字母组成的字符串表示，如“CCATGCTAGTA”。蛋白质序列比较复杂，一般由二十种常见的氨基酸组成，故可用二十个字母{A,R,N,D, C,Q,E,G,H,I,L,K,M,F,P,S,T,W,Y,V}组成的字符串来表示。

（二）序列比对的定义 序列比对（sequence alignment）简单地说就是对不同的生物序列进行分析比较。诸多生物学现象表明很多生物序列可能是由同一序列经过不同的变异演化而来，通过序列比对，找到它们的稳定区域及发生突变的位点和类型，进而可以研究这些突变的生物与医学效应。

在分子生物学中，序列A中的若干基本分子单元发生变异，使序列A变成另外一个序列B，那么称B是A的突变序列。序列的突变有以下四种基本类型：

1.差错型突变：例DNA序列中核苷酸的变异，如ACCGT中的A变为G 2.换位型突变：例DNA序列中核苷酸片段的换位，如ACCGT变为GTACC 3.插入型突变：如在DNA序列GGTTGG的TT间插入AA，变为GGTAATGG 4.删除型突变：如在DNA序列ACAGTTA中删除AG，变为ACTTA 为了便于表示后两种突变类型，我们在序列比对中引入空位“-”,比如 序列1：ACAGTTA 序列2：AC--TTA 表示从序列1中删除AG变异为序列2，也可以说是从序列2中插入AG变异为序列1。序列比对，通俗地讲，就是在两个要比对的序列（可能不一样长）中的适当位置插入空位“-”，以得到两个序列的最优匹配。

为了能定量描述匹配的程度，这就涉及到打分矩阵的概念。以DNA序列为例，考虑到空位“-”，序列的字母表扩大为{,,,,}Vatcg，打分矩阵可表示为

,(,)abVDdab，是一个5阶方阵。

打分矩阵分为得分矩阵与罚分矩阵两种。所谓得分矩阵，就是两序列匹配时取值大，不匹配时取值小；罚分矩阵则相反。在生物信息学中，最简单的打分矩阵为Hamming矩阵或WT矩阵。例

Hamming得分矩阵：1(,)0HabVdababV 相应的Hamming罚分矩阵：0(,)1HabVdababV

WT罚分矩阵：00.770.450.7710.7700.770.4510.450.7700.7710.770.450.770111110Wd 无论是罚分矩阵还是得分矩阵，都在序列比对中起优化作用，因此统称为优化矩阵，记为,(,)abVWwab。如果用A、B分别表示序列A、B在某些位置插入空位“-”后的序列，并且保证插入“-”后的序列长度一样长（记长度为n），

则称1(,)(,)niiiwABwab为序列比对的优化函数。而所谓序列比对就是寻找使优化函数取得最优解的插入“-”后的序列00,AB（取得分矩阵时，最优解乃是取得最大值；取罚分矩阵时，最优解乃是取得最小值）。 另外，我们还可以定义多条序列的序列比对（MSA），有兴趣的读者可参阅相关文献，这里不再赘述。

（三）动态规划算法与MAFFT 二重序列比对的经典算法是由Needleman S B和Wunsch C D于1970年提出的动态规划算法，但Needleman-Wunsch算法只能用于全局序列比对，1981年Smith T F和Waterman M S提出了改进的动态规划算法，使之不仅能用于全局序列比对，也能用于局部序列比对。

Needleman-Wunsch动态规划算法的主要步骤[2]如下： 1.将比对的两条序列制成比对二维表 如果序列为123(,,,)nAaaaa和123(,,,)mBbbbb，那么所构成的二维表如下表。 表3. 1 二维表

2.计算二维表中的元素(,)sij 首先，给出边界上即(0,)sj和(,0)si的分值。采用得分矩阵,(,)abVWwab，假设空位得分为d，而连续空位得分等于dgap，gap为连续空位的长度，则(0,),(,0),(0,0)0sjjdsiids；此时其余位置(,)sij的计算公式如下