东南大学chapter5-1DNA序列分析

主讲人：孙啸

制作人：刘志华东南大学吴健雄实验室第五章DNA序列分析第五章DNA序列分析

DNA序列分析

——基因序列

——基因表达调控信息

寻找基因牵涉到两个方面的工作：识别与基因相关的特殊序列信号预测基因的编码区域结合两个方面的结果确定基因的位置和结构

基因表达调控信息隐藏在基因的上游区域，在组成上具有一定的特征，可以通过序列分析识别这些特征。

第一节DNA序列分析步骤和分析结果评价在DNA序列中，除了基因之外，还包含许多其它信息，这些信息大部分与核酸的结构特征相关联，通常决定了DNA与蛋白质或者DNA与RNA的相互作用。存放这些信息的DNA片段称为功能位点如启动子（Promoter）、基因终止序列（Terminatorsequence）、剪切位点（Splicesite）等。发现重复元素数据库搜索分析功能位点序列组成统计分析综合分析一个基本的DNA序列分析方案功能序列分析的准确性来自于对“功能序列”和“非功能序列”的辨别能力。两个集合：训练集（trainingset）用于建立完成识别任务的数学模型。 测试集或控制集（controlset）用于检验所建模型的正确性。用训练集中实例对预测模型进行训练，使之通过学习后具有正确处理和辨别能力。然后，用模型对测试集中的实例进行“功能”与“非功能”的判断，根据判断结果计算模识别的准确性。收集已知的功能序列和非功能序列实例（这些序列之间是非相关的）训练集（trainingset）测试集或控制集（controlset）建立完成识别任务的模型检验所建模型的正确性对预测模型进行训练，使之通过学习后具有正确处理和辨别能力。进行“功能”与“非功能”的判断，根据判断结果计算模识别的准确性。识别“功能序列”和“非功能序列”的过程

Sn

——敏感性Sp——特异性Tp是正确识别的功能序列数，Tn为正确识别的非功能序列数，Fn是被错误识别为非功能序列的功能序列数，Fp是被错误识别为功能序列的非功能序列数。敏感性和特异性的权衡对于一个实用程序，既要求有较高的敏感性，也要求有较高的特异性。如果敏感性很高，但特异性比较低，则在实际应用中会产生高比率的假阳性；相反，如果特异性很高，而敏感性比较低，则会产生高比率的假阴性。对于敏感性和特异性需要进行权衡，给出综合评价指标。对于一个识别程序准确性可按下式进行综合评价：另一个综合评介指标为相关系数，其计算计算公式为：选择训练集和测试集在检测算法的可行性时，需要从已知的数据中按照不同的方式选择训练集和测试集测试集的构成非常关键在不同的测试集上进行测试可能会得到不同的准确性结果，甚至准确性相差很大。建立标准的功能序列测试集合。如基因转录剪切位点的测试集合、编码区域的测试集合等。第二节核苷酸关联分析对于一个给定的基因组，最简单的计算就是统计DNA序列中各类核苷酸出现的频率。对于随机分布的DNA序列，每种核苷酸的出现是均匀分布的出现频率各为0.25。而真实基因组的核苷酸分布则是非均匀的核苷酸

频率

A0.3248693727808C0.1751306272192G0.1751306272192T0.3248693727808酵母基因组核苷酸出现频率在统计过程中，如果同时计算DNA的正反两条链，则根据碱基配对原则，A和T、C和G的出现频率相同。如果仅统计一条链，则虽然A和T、C和G的出现频率不同，但是非常接近。核苷酸

频率

A0.344C0.155G0.157T0.343单链核苷酸出现频率

基因和其它功能区域在正反两条链上出现的可能性通常一样核苷酸出现频率也不应该有偏差正反两条链在信息的组织结构方面不应该有差别单链上A和T、C和G的出现频率相近。正反两条链碱基互补的原则

单链上A和T、C和G的出现频率相近的解释两联核苷酸频率不同基因组中两个连续核苷酸出现的频率也是不相同的4种核苷酸可以组合成16种两联核苷酸酵母基因组两联核苷酸频率表对酵母基因组两联核苷酸的统计结果其中核苷酸对出现频率最高的达到0.119而出现频率最低的只有0.028令:

Pij

——代表两联核苷酸（i，j）的出现频率

Pi——代表核苷酸i的出现频率则：

Pij’=Pij/(PiPj)

的值反应核苷酸i和j的关联关系如果Pij’=1，则在两个连续的位置上，核苷酸i和j的出现是相对独立的。关联性分析

对于酵母基因组

PA=0.3248PAA=0.1193PAA’=0.1193/（0.3248*0.3248） =1.131>1

表明在两个连续位置上“A”的出现不是独立的，而是相关的。关联性分析

同样，对于相隔一定距离k（k代表核苷酸个数）的两个核苷酸，也可能具有一定的相关性。假设Pij(k)代表核苷酸j出现在核苷酸i之后第k个位置的频率，则可定义一个反应统计相关性的互信息I(k)I(k)值得大小实际上反应了距离为k的两个核苷酸之间的相关性的程度三联核苷酸——基因密码子在进行编码区域识别时，常常需要对三联核苷酸进行统计分析，这实际上是分析密码子的使用偏性。由于密码子的简并性（degeneracy），每个氨基酸至少对应1种密码子，最多有6种对应的密码子。在基因中，同义密码子的使用并不是完全一致的。不同物种、不同生物体的基因密码子使用存在着很大的差异基因密码子的使用与基因编码的蛋白的结构和功能有关，与基因表达的生理功能有着密切的联系蛋白的三级结构与密码子使用概率有密切的关系通过对密码子的聚类分析，可以很清晰地将具有不同三级结构蛋白质的编码基因分成不同的类，而具有相似三级结构蛋白的编码基因则大致聚在同一类中，从而证明基因密码子的使用偏性与蛋白质三级结构具有密切的相关性。在不同物种中，类型相同的基因具有相近的同义密码子使用偏性对于同一类型的基因由物种引起的同义密码子使用偏性的差异较小针对酵母第一染色体的分析结果第三节功能位点分析功能位点（functionalsite）与特定功能相关的位点，是生物分子序列上的一个功能单元，或者是生物分子序列上一个较短的片段。功能位点又称为功能序列（functionalsequence）、序列模式（motif）、信号（signal）等。核酸序列中的功能位点包括转录因子结合位点、转录剪切位点、翻译起始位点等。在蛋白质序列分析中，常使用序列模式这个名词，蛋白质的序列模式往往与蛋白质结构域或者作用部位有关。功能位点示意基因组序列中若干个相邻的功能位点组合形成功能区域（functionalregion）。功能位点分析的任务发现功能位点特征识别功能位点1、利用共有序列搜索功能位点共有序列（consensus）又称一致性片段共有序列是关于功能位点特征的描述，它描述了功能位点每个位置上核苷酸进化的保守性

例如:NTATN利用共有序列进行功能位点分析牵涉到两个方面的问题，如何构造共有序列如何利用共有序列在给定的核酸序列上搜索寻找功能位点，并计算所找到的功能位点的可靠性共有序列具有以下几个方面的特征：（1）共有序列中既有保守的位置，也有可变的位置；（2）任何位置上的核苷酸可以用15种类型之一来表示：核苷酸表示符号符

号含

义说

明GG腺嘌呤AA鸟嘌呤TT胸腺嘧啶CC胞嘧啶RGorA嘌呤YTorC嘧啶MAorC氨基KGorT羧基SGorC强氢键(3个氢键)WAorT弱氢键(2个氢键)HAorCorT非GBGorTorC非AVGorCorA非T(非U)DGorAorT非CNGorAorTorC任意碱基共有序列构造过程:(1)初始化共有序列为一系列可变位置，以“N”代表；(2)在可变位置寻找出现次数最多的核苷酸，并将该位置转化为保守位置；(3)对当前所得到的共有序列进行特异性检查，若通过检查，转（5），否则转（4）；(4)形成与当前共有序列一致的位点子集，转（2）；(5)从原位点集合中删除与当前共有序列一致的位点，若还有剩余位点，则转（1），构造另外的共有序列。TTATGATATATACGCTTGTC

TCCAC

TTATGATATATACGCTTGTC

TCCAC

TNNNNtTATG

tACGC

tTGTC

tCCAC

tTATG

tACGC

tTGTC

tCCAC

TNNNC

[1]

[2]

[3]

[4]

[2]

[3]

NNNNNTNNNN非特异

TNNNC非特异

tACGc

tTGTc

tCCAc

[4]

[2]

tACGc

tTGTc

tCCAc

[3]

TNSNC特异

[5]

Consensus1：

TNSNC剩余位点：

TTATG

ATATA

[5]

Consensus2：

NTATNTNNSC在给定的序列中搜索与共有序列一致的序列片段数据库搜索共有序列表示方法的缺点：是关于序列特征的一种定性描述，对于DNA序列，它能够说明序列每个位置可能出现的碱基类型，但是不能准确地说明各位置上不同类型碱基出现的可能性大小。2、用感知矩阵分析功能位点用权系数描述功能位点各位置上每种核苷酸的相对重要性感知矩阵（或加权矩阵）根据一系列功能位点的多重对比排列结果而建立的其大小为4n4代表碱基的种类数目，n代表功能位点的长度矩阵的每一个元素M(a,j)的值代表第a种核苷酸在功能位点第j个位置上出现的得分,a{A,T,G,C}。123456A18227-319T26142-10G3110-50-19C5-916880感知矩阵示例对于一个序列s=a1a2…an，根据对应位置上核苷酸的类型，取感知矩阵中对应的权值，加和以后得到该序列的得分设S=ATTGCA，则

Ws=1+6+14-5+8+19=43T——功能位点阈值T‘——非功能位点阈值如果WsT，则S是功能位点；如果WsT'，则S是非功能位点。感知矩阵M的构造算法令A+代表功能位点集合

A-代表非功能位点集合过程如下：（1）初始化M为零矩阵；（2）执行过程（3）-（6）的循环；（3）逐步取训练集合中的每个实例Si，如果SiA+，转 过程（4）；如果SiA-，转过程（5）；（4）如果W（Si）T，M不变，否则根据Si的核苷酸分布 将M中所有对应元素的值加1；转（6）；（5）如果W（Si）T‘，M不变，否则根据Si的核苷酸分 布将M中所有对应元素的值减1；转（6）；（6）若训练集合中的所有实例都处理过，则循环结束， 转（7），否则继续执行循环体，直到处理完所有实 例；（7）如果M稳定，则结束；否则转（2）。上述算法反复调整感知矩阵M的元素值，直到M矩阵能够正确识别训练集中的所有功能位点和非功能位点。对于最终得到的感知矩阵，要求其具有敏感性和特异性，每一列上的元素值应该尽可能地有明显的差别，以便反应功能位点各个位置上的特点。与感知矩阵类似，如果令矩阵每一个元素M(a,j)的值代表第a种核苷酸在功能位点第j个位置上出现的概率，则M是一个概率矩阵。假设各个位置上出现的碱基是相互独立的，即任何两个位置上的碱基是不相关的，那么对于给定一个序列s=a1a2…an，可以计算出功能位点序列为s的概率：如果分别统计功能位点和非功能位点，通过计算可以形成两个矩阵M和M’，进一步计算可以判断一个给定的序列究竟属于功能位点，还是属于非功能位点。给定一个序列s=a1a2…an，定义似然比LR(M,M’,s)：在进行功能位点检测时，计算LR(M,M’,s)，并与给定的阈值L比较，如果LR(M,M’,s)>L，则序列s可能是一个功能位点。概率矩阵M和M’的每个元素是一个0和1之间的正数。如果令一个4n新矩阵U的元素(a,j)的值为

log2(M(a,j)/M’(a,j))则矩阵U的每个元素值可能是正值，也可能是负值。实际上，矩阵U就是感知矩阵。第四节隐马尔柯夫模型1、马尔柯夫链（Markovchain）考虑一个具有多个状态的系统S，S={s1,s2,…,s|s|},令S0、S1、…、St为一系列在各个时刻系统状态的变量，即状态链。对于每个1到|S|的整数，它们分别与状态链中的一个状态相联系，并且在任何时刻，这条链都处于一个特殊的状态。当且仅当对于任何t有

则St形成一条马尔柯夫链。简单地说，就是系统未来的状态仅依赖于当前状态。St称为在时刻t系统链的状态。一条马尔柯夫链完全决定于初始分布P(S0)和转换概率Pt=P(St+1|St)。令状态转换矩阵为

F=(fij)

fij代表从状态si移动到状态sj的概率。生物序列可以被描述为一个随机过程的输出，其中对于一个给定的核酸在位置p出现的概率依赖于已占据前面k个位置的核酸，这样一种表示称为k阶马尔柯夫模型。ATCGTAGCAT…….一个序列具有不同的统计性质 （如二目频率或三目周期性）不同的功能区域（如编码区域、非编码区域）对应于不同的马尔柯夫模型。马尔柯夫链在识别CpG岛中的应用CpG岛是一类长度在几百bp的特殊DNA序列，其中CG核苷酸对出现的频率非常高。

ACGCGCGTACGCGAATCpG岛在基因组中有重要的生物学意义，而识别CpG岛有助于在基因组序列中确定我们感兴趣的区域。CpG岛的识别问题表述为：给定一段DNA序列

X=(x1，x2，…,xL)，确定X是否是一个CpG岛。设字母表A={a,t,c,g}，对于字母表中的任何两个字符s、t，定义转换概率为fst=p(xi=t|xi-1=s)，即字符s后面出现字符t的概率。假设{xi}是一个随机过程，随机变量xi的取值仅依赖于xi-1，即对于所有x1,x2,…,xiA，整个序列X的发生概率为为了处理方便，添加两个特殊的字符‘B’（begin）和‘E’（end），使得x0=‘B’,xL+1=‘E’，则上述公式简化为：令fst+为CpG

岛内的字符转换概率

fst-为CpG

岛外的字符转换概率 则X的对数似然得分为上述计算值越大，则X越可能是CpG岛。

CpG岛内部和外部的转换概率

另外一个待解决的问题是：给定DNA序列，确定CpG岛的位置。 直接的方法：对窗口内的子序列计算得分Score(Xk)，具有正值的Xk

就是可能的CpG岛子序列起始位置为k+1，长度为l问题： 事先不知道CpG岛的长度 但是假设CpG岛的长度为l如果l比较大，而真实的CpG岛又比较小，则上述概率计算值不足以证实CpG岛；如果l取值比较小，则难以找出整个CpG岛。这是该算法的最大不足之处，需要考虑其他的算法。HMM2、隐马尔柯夫模型（HMM）功能位点的正则表达式来表示相当于一致性序列这里的正则表达式描述了一个功能位点的构成规律，或者说描述了功能位点各个位置上核苷酸的组成。TGCC—AGG ???ACAC—ATC问题： 对于每个位置，仅仅说明可能的取值，而没有说明各种取值出现的可能性大小例如，用这样的方法无法区分下面两条序列究竟哪一个更可能属于功能位点：

TGCC--AGG ACAC—ATC第一个序列中，假设所有位置上都是取已知出现次数最少的字符而对于第二个序列，所有位置上都是取已知出现次数最多的字符。显然，第一个序列几乎肯定不是功能位点，而第二个序列几乎可以肯定是功能位点，但是用正则表达式表达却无法将两种极端的情况分开。隐马尔柯夫模型可以用于生物序列分析，该模型在生物信息分析方面有重要的应用。一阶隐马尔柯夫模型包括有限数目的系统状态、离散的字母表、状态转换矩阵和字符释放概率。一个HMM模型是一个三元组M=（A，S，Θ）A是字母表S是有限状态集合，每个状态可以释放字母表中的字符。Θ为概率集合，包括两个部分：状态转换概率fkl

k,lS，表示从状态k转换到状态l的概率；字符释放概率，记为ek(b)

kS,bA

表示在状态k下释放出字符b的概率。令路径

=（1,2,…,L）是一个相继状态序列

X=(x1,x2,…,xL)是一个字符序列按下述方式定义状态转换概率和字符释放概率：对于给定的路径，可以按下面的公式计算出产生序列X的概率：

这里，令0为起始状态，L+1为终止状态。例如，对于前面给出的两个序列ACACATC和TGCTAGG，它们的得分分别为：P（ACACATC）=0.81.00.81.00.80.60.40.61.01.00.81.00.8=4.710-2P(TGCTAGG)=0.21.00.21.00.20.60.20.61.01.00.21.00.2=0.002310-2从上述计算结果可以看出，两个序列差别非常大

一个功能位点的HMM模型是通过对一系列的功能位点实例进行机器学习而形成的用这样的模型可以定量的计算一个序列片段是功能位点的可能性计算方法是从模型的第一个状态出发，根据序列的核苷酸组成，将相应的状态值与状态转换值连乘，结束于最后一个状态一个检测CpG岛的HMM模型有8个状态，状态名称和释放的字符为：状态：A+C+G+T+A-C-G-T-

释放字符：ACGTACGT

其中，带有“+”号的状态表示在CpG岛内部，用“-”号标记的状态代表CpG岛外部。假设字符处于CpG岛内的概率是p

处于CpG岛外的概率是q

可以得到状态转换概率CpG岛HMM模型中的状态转换概率解码问题： 给定一个隐马尔柯夫模型M=（A，S，Θ）和一个字符序列X，在M中为X寻找一条最优路径*，在路径中的每一个状态都选择释放一个字符，要求使得P(X|*)最大，记为：

在处理CpG岛问题中，最优路径可以帮助我们寻找CpG岛所在的位置。如果找到最优路径*，则这条路径穿过的“+”状态将对应于CpG岛。3、Viterbi算法求解HMM模型的最优路径基本思想： 动态规划算法给定一个字符序列X=(x1,x2,…,xL)

，以vk(i)代表序列前缀（x1,x2,…,xi）终止于状态k（kS,1≤i≤L）的最可能路径的概率。求解过程如下：（1）初始化

（2）对于每个i=0,…,L-1及每个lS，按下式进行递归计算：

（3）最后，计算序列X终止于状态“end”最可能的路径概率，即P(X|*)的值在正向的递归计算过程中，保持向前推进的反向指针，这样，在正向计算完成后，根据反向指针重构最优路径*。算法的时间复杂度为O(L|S|2),空间复杂度为O(L|S|)。在概率的计算过程中，需要使用大量的乘法运算，在有限计算精度的情况下，会产生误差。如果使用对数值，可以解决这个问题。因此，以vk(i)代表序列前缀（x1,x2,…,xi）终止于状态k（kS,1≤i≤L）的最可能路径的对数得分值，则初值按如下方式设置递归计算及最终得分计算改为（5-26）4、前向概率和反向概率给定一个隐马尔可夫模型M=（A，S，Θ） 一个字符序列X=(x1,x2,…,xL)

要求计算模型M产生X的概率P(X|M)与最优路径问题不一样 前面的问题是在可以产生序列X的各种路径中，选择一条最优路径*，使得P(X|*)最大。 而现在的问题是： 既然有多条路径可以产生序列X，那么模型M产生序列X总的可能性有多大？如果有一条从状态“begin”出发，终止于状态“end”的路径=（0,1,2,…,L,L+1），其中0=“begin”，L=1=“end”,该路径中各状态所释放的字符组成的序列与X相同，则模型M产生X的概率为这里代表所有那些从状态“begin”出发、终止于状态“end”的路径。（5-27）

（5-28）

由于一个HMM模型中可能的路径非常多，穷举每条路径显然是不合适的。下面介绍解决该问题的前向算法（forwardalgorithm）与反向算法(backwardalgorithm)。算法的根本任务是对于每个1≤i≤L及kS,计算概率P(i=k|X，M)。定一个序列X=(x1,x2,…,xL)，令k(i)为释放前缀（x1,x2,…,xi）后到达状态i=k的概率。前向算法初始值的设置与Viterbi算法一样：递归计算过程和最终计算如下：（5-29）

（5-30）

（5-31）

与前向算法相对应，给定一个序列X=(x1,x2,…,xL)，令k(i)为在给定状态i=k下后缀（xi+1,xi+2,…,xL）的概率。反向算法初始化如下：递归计算和终止计算如下：（5-32）

（5-33）

（5-35）

（5-34）

利用正向和反向概率，可以计算出P(i=k|X)。由于HMM的阶数为1，当前的状态仅依赖于前一个状态，则根据条件概率的定义，我们得到解（5-36）（5-37）5、HMM模型的参数估计应用中假设有一个HMM模型，其中的状态转换概率和字符释放概率都是已知的。然而在实际中，情况并非如此。 我们所知道的仅仅是一些实例问题是要根据给定的n个字符串重构M，使得M产生这n个字符串具有最大的概率。由于各个字符串是独立产生的，则若使用对数表示，则目标就是寻找一个*，使得其中这里的n个字符串X(1),X(2),…,

X(n)通常被称为“训练序列”。（5-38）（5-39）（5-40）特殊情况：假设已知与字符串序列X(1),X(2),…,X(n)

相对应的状态序列(1),(2),…,(n)，可以计算从状态k到状态l的转换数Fkl和在状态k下释放字符b的次数Ek(b)。则关于最大似然估计值为：为了避免零概率，当处理数量较少的样本时，需要对Fkl和Ek(b)进行修正：rkl、rk(b)为拉普拉斯修正项，通常情况下为1，可以解释为预先假设的均匀分布。但是在某些情况下，这些修正项可能取其他的值，例如已知状态转换或字符释放的信息，或已有的先验知识。在一般情况，不知道状态序列(1),(2),…,(n)，这时，寻找最优参数集在数学上是一个NP-完全问题，可以用Baum-Welch算法或期望最大（EM）算法解决这个问题。具体的求解算法如下：（1）初始化，给中的参数赋予初值；（2）计算从状态k到状态l转换的期望次数，使用与计算P(X,i=k)时相同的参数（见公式5-36），则（5-45）这样，对所有训练序列X(j)（j=1,…,n）的所有位置i（i=1,…,L(j)，L(j)为序列X(j)的长度）进行求和运算，按下式计算期望值Fkl：其中k(j)(i)是针对序列X(j)的正向计算结果，k(j)(i+1)是反向计算结果。接下来计算在状态k释放字符b的期望次数：（5-47）

（5-46）

（3）重新计算的参数值Fkl和Ek(b)，正如在第一种情况所做的一样（参见公式（5-41）和公式（5-42））；（4）反复执行步骤（2）、（3），直到Score(X(1),X(2),…,

X(n)|)的增量小于给定的一个值很小的参数为止。EM算法保证目标函数Score(X(1),X(2),…,

X(n)|)单调增加，并且概率的对数值接近于0，保证算法收敛。需要注意，收敛的是目标函数，并非是的参数。当目标函数变化趋缓时，的参数值可能波动较大，这意味着算法所得到的结果不稳定。Baum算法的主要问题是目标函数存在若干局部极大，算法不能保证找到全局最大点，算法收敛的点可能是局部极大点。克服局部极大缺陷的一种方法是执行算法若干遍，每次给取不同的初始值。如果算法多次计算结果到达同一个极大点，则可以认为该点是全局最大点。6、基于HMM模型的序列比对可以利用HMM将一个序列与一个序列统计特征（profile）进行比对，从而解决多重比对问题。定义一个长度为L的序列统计特征P是一系列的概率集合ei(b)，ei(b)表示在第i（1≤i≤L）个位置上出现字母表中字符b的概率。这样，在给定条件P下序列X=(x1,x2,…,xL)发生的概率为：如果不考虑“空位”，则X与P的比对得分为：这里，p(b)是字符b的背景出现频率。（5-49）定义一个基本HMM模型，有L个“匹配”状态M1,M2,…,ML，它们对应与统计特征的匹配。所有这些状态顺序连接起来，即状态Mj连接到后继Mj+1，如图5.5所示。其中从状态Mj释放字符b的概率为ej(b)。为了在比对中允许插入“空位”的操作，在上述基本模型中加入“插入”状态I1,I2,…,IL，并假设每个插入状态Ij,有一个来自相应匹配状态Mj的连接，有一个到匹配状态Mj+1的连接，还有一个自循环连接。根据“空位”的惩罚原则，给这些状态转换赋予适当的概率。（5-50）

同样，为了允许“删除”操作，可以进一步假如“删除”状态D1,D2,…,DL，这些状态不能释放任何字符。删除状态依然顺序连接，同时增加从Dj到Ij的连接及从Ij到Dj+1的连接。完整的HMM模型如图5.5所示：D1D2D3I2I3I4BeginEndM1M2M3I1

图5.5用于序列多重比对的HMM模型

下面介绍一种Viterbi类似算法，将X=(x1,x2,…,xm)与长度等于L的统计特征P进行比对。对于每一个1≤j≤L和1≤i≤m，定义：（1）vjM(i)代表子序列(x1,x2,…,xi)与HMM模型P的匹配对数得分值，该匹配以状态Mj释放字符xi作为最后操作;

（2）vjI(i)代表子序列(x1,x2,…,xi)与HMM模型P的匹配对数得分值，该匹配以状态Ij释放字符xi作为最后操作;

（3）vjD(i)代表子序列(x1,x2,…,xi)与HMM模型P的匹配对数得分值，该匹配以状态Dj结束（不释放任何字符）。模型P中特殊状态“begin”的初始值为：为了计算vjM(i)、vjI(i)和vjD(i)的值，使用Viterbi算法中的相同技术，但现在的模型有两个特点：（1）模型中的每一个状态最多只有3个引入连接，如上图所示；（2）“删除”状态不释放任何字符。（5-51）“匹配”状态Mj的三个前驱同属于上一层,即j-1层，有“插入”状态Ij的三个前驱属于同一层,