一种新的基于结构信息的双生物序列比对方法

摘 要: 用于生物序列比对的经典动态规划算法是用一个固定的替换矩阵来逐点计算生物序列间的代价, 这些方法可用来发 现具有最大计分值的比对结果, 但实际上, 则更加倾向于考虑生物序列中所隐含的结构或功能信息. 本文用可变长马尔科夫链 方法来发现生物序列中所隐含的结构或功能信息子片断并定义其权值, 最后提出一个新的基于结构信息的生物序列比对方法. 关 键 词: 序列比对; 概率后缀树; 可变长马尔科夫链 中图分类号: T P301 文献标识码: A 文 章 编 号: 100021220 ( 2006) 0120085205
( 2 ) 树中的边由字符集中的每一个字符所表示, 而树中 的节点则用一个二元组所表示 ( s, Χ s ) , 其中 s 是由此节点到
1, m j 0, m j s + 1 s + 1
定义 2. 给定一个间隙罚分函数 g : N →R 和一个在索引 二元组上的权函数 w : ∑ × ∑ →R , 则一个序列比对的分值 S core 定义如下:
S core = w ( sx 1 , ty 1 ) +
…m j = s …m j ≠s
( 4)
1, x N = s ∨y N = t ; ( 2) m in (x i+ 1 - x i , y i+ 1 - y i ) = 1, Π i= 1…, N - 1.
值来修正 N eed lem an 2 W un sch 算法. 我们的方法和 N eed le2
m an 2 W un sch 算法有相似之处: 随着比对的进行, 子序列逐个
对于任意两个记为x的生物序列来说不失一般性假定npst的构造算法可以线性时间完同时可以对其重要子片断进行预处理以得到其分值长度和权值对于每一个子片断我们可以在时间内完成因此预处理总共的时间为seglengthseg又因为seglengthseg所以这种预处理的时间复杂度也为o10的计算主要是基于动态规划方法所以其时间复杂度也应为o总之我们的算法总时间复杂度为o2算法的实验结果基于以上思想我们用一些生物序列进行了比对实验了进行分析还用clustal5b对同样的生物序列进行了比对所采用的生物序列均来自embl而参与比较的结构信息从prosite数据库中提取出来比对时采用的基本计分矩阵为pam250为其中两组比对的比较结果比对结果由genedoc输出进行的是p07146和p00765的比对比对结果表明图1p07146和p00765比对结果图的上半部分为clustal的比对结果下半部分为本算法的比对结果相同的残基数目是104进行比对表明残基数目是100属于保守的残基数目是27另外查询prosite数据库可知两序列含有相同的模式ps00134ps00135我们的算法也成功的进行了匹配
86
小 型 微 型 计 算 机 系 统
2006 年
较重要的子片断及其权值, 最后再用这些重要子片断及其权
2 序列比对和 PST 的形式化描述
设∑为一非空有限字符集, 而 s= (s i ) i= 1… s 、 t= ( t i ) i= 1… t 是在∑上的有限长的字符串 ( 本文为生物序列) , 则有, 定义 1. 序列比对一般就是符合如下条件的索引二元组 (x i , y i ) i= 1…N 的有序集: ( 1) (x i , y i ) ∈{1, …, s } ×{1, …, t }, 且 x 1 = 1 ∨y 1 =
字符进行比较. 但是一旦发现一些重要子片断匹配了, 那么其 比对的分值就要按照其权值进行修改, 以下的各小节将对其 进行详细的描述.
3. 1 PST 的构造
给定一个生物序列M = m 1m 2 …m n21m n , 其马尔科夫链的 步长最大为 L , 有一任意子序列 s = s1 s2 … s l21 s l , 有如下定 义: X j (s) = 则:
1 引 言
随着人类基因组测序的完成和其他物种基因测序工作的 相继完成, 产生了大量的基因数据, 这就要求开发一些具有高 吞吐率的技术来分析. 生物序列比对将未知序列同整个数据 库中的已知序列进行分析比较, 为生物学家预测未知生物序 列的功能提供了一条新的途径, 如果一个新的生物序列经过 比对知道和一个已知其功能的生物序列非常相似, 那么在这 个新的生物序列上很有可能也包含着一个和已知生物序列中 类似的基因并且其功能也相似. 双序列比对主要是判断两个 序列之间是否有足够的相似性, 从而判定二者之间是否具有 同源性. 但如何来计算相似性有很多方法, 例如基于固定替换 矩阵的全序列比对 [1 ] 和允许有空格或不允许有空格的局部序 列比对 [2, 3 ]. N eed lem an & W un sch 的经典算法就是用动态规 划方法来得到一个具有最大分值的序列比对, 它所产生的分 值可以作为一个依据用来判断生物序列之间是否具有同源 性. 这种类型的算法都使用了一个固定的替换矩阵, 并没有 考虑生物序列间的相对位置. 相似性的分值是采用一个固定 的替换矩阵和罚分函数逐点累加求得. 但在实际应用中, 生物 学意义上的结构或功能信息比序列的简单对齐更加重要, 所
m- 1
( 6)
如果序列比对不以 ( 1, 1 ) 开始e = w ( sx 1 , ty 1 ) 2 G ( ( 0, 0 ) , ( x 1 , y 1 ) ) +
(w ( sx i , ∑ i= 2
N
ty i )
且由 [ 8 ] 可知, 我们可以在 O ( n ) 时间内构造出其概率后 缀树, 而与其马尔科夫链的步长 L 无关, 这也就是说, 构造一 个 PST 的时间是与生物序列的长度成线性关系; 一旦 PST 被构建, 可以在 O (m ) 去预测一个长度为 m 的字符串的出现 概率.
(w ( sx i , ∑ i= 2
N
ty i ) - G ( ( x i- 1 , y i- 1 ) ,
(x i , y i ) ) )
( 1)
P rob (s) =
1
n- L
n- 1
∑ j ( s) , j= L X
( 5)
其中, G ( (x i- 1 , y i- 1 ) , (x i , y i ) ) =
以在生物学意义上序列的最优比对并不一定总是数学意义上 的具有最大分值的序列比对. 为了解决这个问题, 有两种基本 方法: 一种就是求生物序列比对的近优解或次优解, 这些并不 都是最优解的集合中就有很大可能包含了生物学意义上序列 的最优比对; 另外一种基本方法就是在进行生物序列比对记 分时不但要考虑替换矩阵, 更要考虑生物序列中所隐含的结 构或功能信息. 因此, N eed lem an & W un sch 的经典动态规划 算法在实际应用就出现了许多变种. 例如 W ilbu r and L ipm an 在 [ 4 ] 中就提出了按照生物序列的上下文来估算各种匹配计 分的思想; 在文献 [ 5 ] 中还提出了利用业已存在的包含结构或 功能信息的模式库 ( PRO S IT E ) 中的信息, 来进行生物序列比 对, 这就把序列比对算法和模式匹配结合起来, 但是各种模式 的具体权值则需要手工指定, 这就影响了其更加广泛的应用. 在本文中, 我们给出了一个自动考虑生物序列模式匹配 问题的新算法, 和文献 [ 5 ] 不同的是, 我们不要从任何模式库 中提取模式, 而是直接分析生物序列以从中获取可能包含结 构或功能信息的子片断, 采用概率后缀树 (p robab ilistic suffix
ter (H ef ei) , H ef ei 230026, C h ina )
Abstract: T he cla ssica l dynam ic p rog ramm ing a lgo rithm s fo r a lignm en ts a re ba sed on ed ition co sts com p u ted add itiona lly po si2 tion by po sition, acco rd ing to a fixed sub stitu tion m a trix. T hese m ethod s find the g loba l a lignm en t w ith m ax i m a l sco re bu t in p ractice, w e favo r the info rm a tion of the structu re o r the function of the sequences con sidered. . structu re o r the function of the sequences Key words: sequence a lignm en t; p robab ilistic suffix trees (PST ) ; va riab le m em o ry m a rkov In th is p ap er, w e p resen t a m ethod con sisting in com b in ing the techn iques of dynam ic p rog ramm ing and V a riab le M em o ry M a rkov, w h ich revea ls the
3. 2 权值的确定
- G ( (x i- 1 , y i- 1 ) , (x i , y i ) ) ) - G ( (x N , y N ) , ( s + 1, t + 1) ) ( 2)
在使用 PST 获得生物序列中的一些重要子片断后, 由于 这些子片断在进行比对时其重要性肯定是不同的, 如何确定 其权值对于我们的算法的正确性和有效性都是至关重要的一 步. 确定权值的方法是基于生物序列的统计学特性, 进一步来 说, 就是用给定特定子片断后关于每一个字符的条件概率分 布来确定, 而这种条件概率分布确实代表了生物序列中的结 构或功能信息 [11 ]. 给定一个 PST 和其上的条件概率分布 P, 我们可以预测 子片断 s = s1 s2 …s l21 s l 的概率:
g ( x i - x i- 1 ) , x i - x i- 1 > 1
而对于任意 Ρ∈∑, 有 P rob ( Ρ s) =
0, x i - x i- 1 = 1∧y i - y i- 1 = 1
g (y i - y i- 1 ) , y i - y i- 1 > 1
∑j = L X j + 1 (sΡ) m- 1 ∑j = L X j (s)
P ( s ) = P ( s1 ) ×P ( s2 s1 ) × P ( s3 s1 s2 ) ×…× P ( s l s1 s2 … s l- 1 ) =
通常, 函数 g 是一个线性函数, 对于一个长度为 n 的间 隙, 类似于 g (n ) = a + bn , 其中 a 称为间隙开放罚分, b 称为 间隙延伸罚分, 且 a > b; 函数 w 则采用一些经验的记分矩 阵, 如 PAM 250 分值矩阵、 BLO SUM 62 分值矩阵等. 定义 3. 一个 PST 树其实就是: ( 1) 在非空有限字符集∑上 ( 对于蛋白质序列, 其字符集 大小为 20, 而对于 DNA 序列则其大小为 4 ) 所构造的 ∑ 叉 树;
第 27 卷 第 1 期 2006 年 1 月
小型微型计算机系统 M I N I- M ICRO SYST EM S
V o l127 N o. 1 J an. 2006
一种新的基于结构信息的双生物序列比对方法
司秀华, 陈国良
( 中国科学技术大学 计算机科学技术系, 国家高性能中心, 安徽 合肥 230026) E2 m ail: sixh@ u stc. edu
[6, 12 ] 来定义并发现具有重要功能或结构信 tree, 简记为 PST )
息的子片断及其权值, 这样进行生物序列比对就综合考虑了 生物序列间的功能或结构信息.
收稿日期: 2004209213 基金项目: 国家 “八六三” 高技术研究发展计划基金项目 (2002 AA 104560 和 2001AA 111041) 资助 作者简介: 司秀 华, 男, 1975 年生, 博士研究生, 研究方向为生物序列比对、 并行算法; 陈国良, 男, 1938 年生, 教授, 博士生导师, 中国科学院院士, 研究方向为并 行算法、 并行数据库、 体系结构、 生物信息学等.
Novel Structure - Ba sed Pa irw ise Sequence A l ignm en t
S I X iu 2hua, CH EN Guo 2liang
(D ep a rtm en t of C om p u ter S cience and T echnology , U n iversity of S cience and T echnology of C h ina , N a tiona l H ig h P erf orm ance C om p u ting C en 2