第四章序列分析

第四章（1）序列分析第四章序列比较序列比较的根本任务是：发现序列之间的相似性辨别序列之间的差异目的： 相似序列相似的结构，相似的功能 判别序列之间的同源性 推测序列之间的进化关系第一节序列的相似性同源（homology）-具有共同的祖先直向同源（Orthologous

）共生同源（paralogous

）相似（similarity）

—同源序列一般是相似的

—相似序列不一定是同源的

—进化趋同（同功能）直向同源（a1inspeciesI,a1inspeciesII)共生同源（a1anda2inspeciesI）进化趋同水平转移基因复制序列的相似性描述定性的描述定量的数值相似度距离序列比较的基本操作是比对（Alignment）

两个序列的比对是指这两个序列中各个字符的一种一一对应关系，或字符的对比排列。设有两个序列：

GACGGATTAG，GATCGGAATAG

Alignment2:

GACGGATTAG GATCGGAATAGAlignment1:

GACGGATTAG GATCGGAATAG1、字母表和序列字母表4字符DNA字母表：{A,C,G,T}扩展的遗传学字母表或IUPAC编码单字母氨基酸编码符号含义说明GGGuanineAAAdenineTTThymineCCCytosineRGorAPurine

YTorCPyrimidine

MAorCAminoKGorTKeto

SGorCStronginteraction(3Hbonds)WAorTWeakinteraction(2Hbonds)HAorCorTNot-GBGorTorCnot-AVGorCorAnot-T(not-U)DGorAorTnot-CNGorAorTorCAny扩展的遗传学字母表或IUPAC编码1、字母表和序列特定的符号

—代表字母表A*—代表由字母表A中字符所形成的一系列有限长度序列或字符串或序列的集合

a、b、c—代表单独的字符

s、t、u、v—代表A*中的序列

|s|—代表序列s的长度为了说明序列s子序列和s中单个字符，在s中各字符之间用数字标明分割边界例如，设s=ACCACGTA，则s可表示为0A1C2C3A4C5G6T7A8

i:s:j指明第i位或第j位之间的子序列,

当然，0

i

j

|s|。

子序列0:s:i

称为前缀，即prefix(s,i)子序列

i:s:|s|称为后缀，即suffix(s,|s|-i+1)

i:s:i—为空序列j-1:s:j—表示s中的第j个字符，简记为sj子序列与子串子序列：选取s中的某些字符（或删除s中的某些字符）而形成s的子序列例如：TTT是ATATAT的子序列。

s的子串： 是由s中相继的字符所组成。例如：

TAC是AGTACA的子串， 但不是TTGAC的子串（是子序列）。

子串是子序列 子序列不一定是子串字符串操作字符串连接操作： 两个序列s和t的连接：

s++t

例如：

ACC++CTA=ACCCTA字符串k操作—删除字符串两端的字符其定义如下：

prefix(s,l)=sk|s|-l suffix(s,l)=k|s|-ls

i:s:j=ki-1sk|s|-j

序列比较可以分为四种基本情况：

（1）两条长度相近的序列相似

找出序列的差别

（2）判断一条序列的前缀与另一条序列的后缀相似

（3）判断一条序列是否是另一条序列的子序列

（4）判断两条序列中是否有非常相似的子序列2、编辑距离（EditDistance)GCATGACGAATCAG

TATGACAAACAGC

GCATGACGAATCAG

TATGAC-AAACAGC说明两条序列的相似程度——〉定量计算两条序列的相似程度的定量计算相似度，它是两个序列的函数，其值越大，表示两个序列越相似两个序列之间的距离。距离越大，则两个序列的相似度就越小

字符编辑操作（EditOperation）字符编辑操作可将一个序列转化为一个新序列Match（a，a）Delete（a，-）Replace（a，b）Insert（-，b）直接距离计算的不足扩展的编辑操作ACCGACAATATGCATA

ATAGGTATAACAGTCAACCGACAATATGCATA

ACTGACAATATGGATA第二条序列头尾颠倒CTAGTCGAGGCAATCTGAACAGCTTCGTTAGT？反向互补序列RNA发夹式二级结构

3、通过点矩阵进行序列比较

“矩阵作图法”或“对角线作图”

→序列1→→序列2→实例→序列1→→序列1→自我比较滑动窗口技术两条序列中有很多匹配的字符对，因而在点矩阵中会形成很多点标记。滑动窗口技术使用滑动窗口代替一次一个位点的比较是解决这个问题的有效方法。假设窗口大小为10，相似度阈值为8，则每次比较取10个连续的字符，如相同的字符超过8个，则标记基于滑动窗口的点矩阵方法可以明显地降低点阵图的噪声，并且明确无误的指示出了两条序列间具有显著相似性的区域。（a）对人类（Homosapiens）与黑猩猩（Pongo

pygmaeus）的β球蛋白基因序列进行比较的完整点阵图。（b）利用滑动窗口对以上的两种球蛋白基因序列进行比较的点阵图，其中窗口大小为10个核苷酸，相似度阈值为8。(a)(b)具有连续相似区域的两条DNA序列的简单点阵图4、序列的两两比对序列的两两比对 （PairwiseSequenceAlignment）

按字符位置重组两个序列，使得两个序列达到一样的长度

s: AGCACAC

A AG

CACACAt: A

CACACTA ACACACT

A—————————————————————————— Match(A,A) Match(A,A) Delete(G,-) Replace(G,C) Match(C,C) Insert(-,A) Match(A,A) Match(C,C) Match(C,C) Match(A,A) Match(A,A) Match(C,C) Match(C,C) Replace(A,T) Insert(-,T) Delete(C,-) Match(A,A) Match(A,A)图3.6序列AGCACACA和ACACACTA的两种比对结果Alignment-1 Alignment-2不同编辑操作的代价不同为编辑操作定义函数w，它表示“代价（cost）”或“权重（weight）”。对字母表

中的任意字符a、b，定义

w(a,a)=0w(a,b)=1a

b w(a,-)=w(-,b)=1

也可以使用得分（score）函数来评价编辑操作

p(a,a)=1p(a,b)=0a

b p(a,-)=w(-,b)=-1

概念：两条序列s和t的比对的得分（或代价）等于将s转化为t所用的所有编辑操作的得分（或代价）总和；s和t的最优比对是所有可能的比对中得分最高（或代价最小）的一个比对；s和t的真实距离应该是在得分函数p值（或代价函数w值）最优时的距离。

例如：s: AGCACAC

At: A

CACACTA

cost=2

s: AGCACAC

At: A

CACACTAscore(s，t）=5序列比对的目的是寻找一个得分最大（或代价最小）的比对。5、打分矩阵（WeightMatrices）（1）核酸打分矩阵设DNA序列所用的字母表为

={A，C，G，T}a.等价矩阵b.BLAST矩阵c.转移矩阵（transition，transversion）（嘌呤：腺嘌呤A，鸟嘌呤G；嘧啶：胞嘧啶C，胸腺嘧啶T）ATCGA1000T0100C0010G0001ATCGA5-4-4-4T-45-4-4C-4-45-4G-4-4-45ATCGA1-5-5-1T-51-1-5C-5-11-5G-1-5-51表3.1等价矩阵表表3.3转移矩阵表3.2BLAST矩阵（2）蛋白质打分矩阵（i）等价矩阵（ii）氨基酸突变代价矩阵GCM（iii）疏水矩阵

（iv）PAM矩阵（PointAcceptedMutation）（v）BLOSUM矩阵（BlocksAminoAcidSubstitutionMatrices）其中Rij代表打分矩阵元素i、j分别代表字母表第i和第j个字符。PAM矩阵（PointAcceptedMutation）

基于进化的点突变模型

一个PAM就是一个进化的变异单位,即1%的氨基酸改变

这类矩阵里列出同源蛋白质在进化过程中氨基酸变化的可能性。这类矩阵式基于进化原理的证据：编码相同蛋白质的基因随着进化发生分歧，相似度降低。

科学用得多

矩阵集合-----PAM-N如，PAM120矩阵用于比较相距120个PAM单位的序列。一个PAM-N矩阵元素（i，j）的值：反应两个相距N个PAM单位的序列中第i种氨基酸替换第j种氨基酸的频率。针对不同的进化距离采用PAM矩阵序列相似度=40%50%60%

|||打分矩阵=PAM120PAM80PAM60PAM250→14%-27%

归一化打分实例BLOSUM62第二节两两比对算法1、序列两两比对基本算法直接方法—生成两个序列所有可能的比对，分别计算代价函数，然后挑选一个代价最小的比对作为最终结果。本质问题：优化动态规划寻优策略动态规划算法（DynamicProgramming）最短路经问题起点终点C1C2W1

W2路径1：C1+w1？路径2：C2+w2？

取最小值！算法求解:

从起点到终点逐层计算利用动态规划方法求解序列的两两比对起点终点ATTC………CGAAGA

AGTC………GAAGGTATTC………CGAAGAGTC………GAAGGAT+（1）ATTC………CGAAGAAGTC………GAAGG-T+（2）ATTC………CGAAGAGTC………GAAGGTA-+（3）求解过程起点终点ATTC………CGAAGA

AGTC………GAAGGT从两个序列前端开始逐步推进直到两个序列的末端。序列S:i-1ii+1序列t:j-1jj+1序列S:i-1ii+1序列t:j-1jj+1Case1:匹配（si，tj)中间过程：比对0:S:i

与0:T:j序列S:i-1ii+1序列t:j-1jj+1序列S:i-1ii+1序列t:j-1jj+1Case2:删除（si，

-）

序列S:i-1ii+1序列t:j-1jj+1序列S:i-1ii+1序列t:j-1jj+1Case3:插入（

-，tj

）

设序列s、t的长度分别为m和n。考虑两个前缀

0:s:i

0:t:j

假如已知序列0:s:i和0:t:j所有较短子列的最优比对，即已知： （1）0:s:(i-1)

和0:t:(j-1)

的最优比对 （2）0:s:(i-1)

和0:t:j

的最优比对 （3）0:s:i

和0:t:(j-1)

的最优比对则0:s:i和0:t:j的最优比对一定是上述三种情况之一的扩展（1）替换（si，tj）或匹配（si，tj)，这取决于si

是否等于tj

；（2）删除（si，

-）；（3）插入（-，tj

）。

令：为序列0:s:i和与序列0:t:j比对的得分按下述方法求解

其初值为：fori=1,2,......,mforj=1,2,......,n距离矩阵按照上述方法，对于给定的得分函数p，两个序列所有前缀的得分定义了一个(m+1)

(n+1)的距离矩阵

D=(di,j)

其中di,j=S(0:s:i,0:t:j)di,j的计算公式如下：di,j

最小值的三种选择决定了各矩阵元素之间的关系，用下图表示：di,jdi,j-1di-1,jdi-1,j-1距离矩阵元素di,j

的计算

S(0:s:i,0:t:j)S(0:s:i-1,0:t:j)S(0:s:i-1,0:t:j-1)S(0:s:i,0:t:j-1)动态规划算法计算过程：计算过程从d0,0开始可以是按行计算，每行从左到右，也可以是按列计算，每列从上到下。当然，任何计算过程，只要满足在计算di,j时

di-1,j、di-1,j-1、和di,j-1都已经被计算这个条件即可。在计算di,j后，需要保存di,j是从di-1,j、di-1,j-1、或di,j-1中的哪一个推进的，或保存计算的路径，以便于后续处理。上述计算过程到dm,n结束。

最优路径求解： 与计算过程相反从dm,n开始，反向前推。假设在反推时到达di,j，根据保存的计算路径判断di,j究竟是根据di-1,j、di-1,j-1、和di,j-1中的那一个计算而得到的。找到这个点以后，再从此点出发，一直到d0,0为止。走过的这条路径就是最优路径（即代价最小路径），其对应于两个序列的最优比对。

计算过程：（1）初始化

计算过程：（2）反复计算按列计算

计算过程：（2）反复计算按行计算其他方式

计算过程：（3）求最佳路径

t

s

ACACACTAAGCACACA例：s=AGCACACAt=ACACACTA得分矩阵D（9×9）

t

s

ACACACTA0-1-2-3-4-5-6-7-8A-1G-2C-3A-4C-5A-6C-7A-8初始化

计算d（2，2）t

s

ACACACTA0-1-2-3-4-5-6-7-8A-110-1-2-3-4-5-6G-201C-3A-4C-5A-6C-7A-8

最终的得分矩阵及序列比对t

s

ACACACTA0-1-2-3-4-5-6-7-8A-110-1-2-3-4-5-6G-2010-1-2-3-4-5C-3-11110-1-2-3A-4-2021210-1C-5-3-1132321A-6-4-2024333C-7-5-3-113543A-8-6-4-202455AGCACAC

A|||||||A

CACACTA序列长度的影响：令cw(s,t)表示两个长度分别为m和n的序列的相似性得分设cw(s,t)=99如果

m=n=100->则可以说这两个序列非常相似但如果m=n=1000，则仅有10%相同相对长度的得分

sim（s，t）=2*cw(s,t)/(m+n)算法分析： 数据结构di,j

空间复杂度：O(mn)

时间复杂度：O(mn)2、子序列与完整序列的比对----AGCT----ATGCAGCTGCTT目标： 使S(s,i:t:j)最大序列S:序列t:ij不计前缀0:t:i的得分,也不计删除后缀的j+1:t:|t|得分不计前缀0:t:i的得分——处理第一行t

s

ACACACTA000000000AGCACACA不计删除后缀的j+1:t:|t|得分

——处理最后一行dm,,jdm,,j-1dm-1,,jdm-1,,j-1

S(0:s:i,0:t:j)S(0:s:i-1,0:t:j)S(0:s:i-1,0:t:j-1)S(0:s:i,0:t:j-1)不计代价距离矩阵初始化时，对第一行进行如下处理：

d0,j=0for0

j

n

最后一行的计算应该是：同样，dm,n依然是最优局部比对的得分，而匹配的子列i:t:j

按如下方式寻找： （1）j=min{k

dm,k=dm,n}

（2）反推比对路径，最终通过斜线（非空位）到达（0，i）。（3-10）（3-11）3、寻找最大的相似子序列目标： 使dw(i:s:j,i’:t:j’)最大序列S:序列t:i’j’ij数据结构：

（m+1）

（n+1）的矩阵D

但是，对数组元素含义解释与基本算法有所不同

每个元素的值代表序列0:s:i某个后缀和序列0:t:j

某个后缀的最佳比对。

这种局部比对不计前缀的得分，所以新的边界条件是：d0,j=0 for0

j

n （3-12）di,0=0 for1

i

m另外，由于0:s:i

和0:t:j

总有一个得分为“0”的空后缀比对，因此矩阵D中的所有元素大于或等于“0”，于是，新的递归计算公式为：(3-13)寻找最佳比对的子序列在矩阵中找最大值该值就是最优的局部比对得分该值对应的点为序列局部比对的末点然后反向推演前面的最优路径，直到局部比对的起点。

TATA||||TATA4、准全局比较所谓准全局比较就是在评价序列比对时不计终端“空缺”（endspace，或空位）的得分或代价序列1长度为8序列2长度为18（a）6个匹配，1个失配，1个空位（b）8个匹配情况1：不记s后面的空位与t后缀比对的得分 在矩阵di,j中取最后一行的最大值，即：(3-14)序列S:序列t:i’j’ij空位后缀情况2：不记s前面的空位与t前缀比对的得分 将矩阵di,j中的第一行各元素值置为“0”序列S:序列t:i’j’ij空位前缀 情况3：……

情况4：……

半全局比较算法与基本算法在计算di,j时的区别归纳为下列四个方面：（1）第一行初始值为“0”，表示不计第一个序列的前端空位；（2）寻找最后一行的最大值，表示不计第一个序列的末端空位；（3）第一列初始值为“0”，表示不计第二个序列的前端空位；（4）寻找最后一列的最大值，表示不计第二个序列的末端空位。

对于最后一行和最后一列的另一种处理办法是：——最后一行的横向移动不被空位罚分——最后一列的纵向移动也不被罚分这样，就可以允许在两条序列终端自由存在空位。当矩阵D所有元素计算完以后，其右下角得值即为两条序列最终准全局比对的得分。ACACTGATCG||||||ACACTG5、关于连续空位的问题K阶空位

—K个连续的空位字符“-” ATG-A-T-C-A-G ATG-----ATCAG ATGCAGTGCAATG ATGTTTTTATCAG生物学意义“插入”或“删除”突变突变次数连续空位可能对应于一次突变非连续空位对应于多次突变对于连续空位的代价是一个线性的函数。设p（k）代表空位得分函数，其中k是连续空位的个数，则：

p（k）=-bk

这里b（>0）是单个“空位”得分的绝对值。处理方法：任何一个比对可以被唯一地分为若干个相继的块。有三类块： （1）两个字符的比对 （2）与序列s空位进行比对的t的最大连续字符序列 （3）与序列t空位进行比对的s的最大连续字符序列为比较序列s（长度为m）和序列t（长度为n），我们使用三个（m+1）

（n+1）的矩阵各矩阵第一行和第一列初始值的设定如下：（3-16）（3-17）（3-18）递归计算过程如下：(3-20)

for1

k

jfor1k

jfor1k

ifor1k

i(3-19)(3-21)上述算法的时间复杂度为O(n3)。比起标准算法，其多花的时间主要用于处理连续的空位。那么，是否可以改进连续空位的得分函数，而使得算法的时间复杂度降低为O(n2)呢？如果认为k个连续空位比k个孤立空位出现的可能性更大，则

p(k)

kp(1) （3-22）或更一般地，

p(k1+k2+…+kn)

p(k1)+p(k2)+…+p(kn) (3-23)可以用下式重新计算连续“空位”的得分：

p(0)=0 （3-24）

p(k)=–h–g(k-1),k

1 （3-25）

h

0，g

0，h

g。依然用A、B、C三个矩阵，各自的意义如下：ai,j

——0:s:i

与0:t:j

最优比对的得分，该比对以si

和tj

匹配结束bi,j——0:s:i

与0:t:j

最优比对的得分，该比对以空位和tj

匹配结束ci,j——0:s:i与0:t:j

最优比对的得分，该比对以si

和空位匹配结束对个矩阵元素的初始化工作按以下公式进行：

a0,0=0ai,0=-

for1i