《chapter4DNA序列分析》由会员分享，可在线阅读，更多相关《chapter4DNA序列分析（178页珍藏版）》请在金锄头文库上搜索。
1、第第4章章 DNA序列分析序列分析 DNA序列分析基因序列基因表达调控信息 寻找基因牵涉到两个方面的工作寻找基因牵涉到两个方面的工作 ：识别与基因相关的特殊序列信号预测基因的编码区域结合两个方面的结果确定基因的位置和结构 基因表达调控信息隐藏在基因的上游区域，在组成上具有基因表达调控信息隐藏在基因的上游区域，在组成上具有一定的特征，可以通过序列分析识别这些特征。一定的特征，可以通过序列分析识别这些特征。 第第1节节 DNA序列分析步骤和分析结果评价序列分析步骤和分析结果评价在DNA序列中，除了基因之外，还包含许多其它信息，这些信息大部分与核酸的结构特征相关联，通常决定了DNA与蛋白质或者DNA
2、与RNA的相互作用。存放这些信息的DNA片段称为功能位点如启动子（Promoter）、基因终止序列（Terminator sequence）、剪切位点（Splice site）等。发现重复元素数据库搜索分析功能位点序列组成统计分析综合分析一个基本的DNA序列分析方案功能序列分析的准确性来自于对“功能序列”和“非功能序列”的辨别能力。两个集合： 训练集（training set）用于建立完成识别任务的数学模型。测试集或控制集（control set）用于检验所建模型的正确性。用训练集中实例对预测模型进行训练，使之通过学习后具有正确处理和辨别能力。然后，用模型对测试集中的实例进行“功能”与“非功能
3、”的判断，根据判断结果计算模识别的准确性。收集已知的功能序列和非功能序列实例收集已知的功能序列和非功能序列实例（这些序列之间是非相关的（这些序列之间是非相关的 ）训练集训练集（training set）测试集或控制集测试集或控制集（control set）建立完成识别任务的模型建立完成识别任务的模型检验所建模型的正确性检验所建模型的正确性对预测模型进行训练，对预测模型进行训练，使之通过学习后具有使之通过学习后具有正确处理和辨别能力。正确处理和辨别能力。进行进行“功能功能”与与“非功能非功能”的的判断，根据判断结果计算判断，根据判断结果计算模识别的准确性。模识别的准确性。识别识别“功能序列功能序
4、列”和和“非功能序列非功能序列”的过程的过程 Sn 敏感性敏感性Sp特异性特异性Tp是正确识别的功能序列数，Tn为正确识别的非功能序列数，Fn是被错误识别为非功能序列的功能序列数，Fp是被错误识别为功能序列的非功能序列数。敏感性和特异性的权衡对于一个实用程序，既要求有较高的敏感性，也要求有较高的特异性。如果敏感性很高，但特异性比较低，则在实际应用中会产生高比率的假阳性；相反，如果特异性很高，而敏感性比较低，则会产生高比率的假阴性。对于敏感性和特异性需要进行权衡，给出综合评价指标。 对于一个识别程序准确性可按下式进行综合评价：另一个综合评介指标为相关系数，其计算计算公式为：选择训练集和测试集在检
5、测算法的可行性时，需要从已知的数据中按照不同的方式选择训练集和测试集 测试集的构成非常关键在不同的测试集上进行测试可能会得到不同的准确性结果，甚至准确性相差很大。 建立标准的功能序列测试集合。如基因转录剪切位点的测试集合、编码区域的测试集合等。 第第2节节 核苷酸关联分析核苷酸关联分析对于一个给定的基因组，最简单的计算就是统计DNA序列中各类核苷酸出现的频率。对于随机分布的DNA序列，每种核苷酸的出现是均匀分布的出现频率各为0.25。而真实基因组的核苷酸分布则是非均匀的核苷酸核苷酸 频率频率 A0.3248693727808 C0.1751306272192 G0.1751306272192
6、T0.3248693727808 酵母基因组核苷酸出现频率酵母基因组核苷酸出现频率在统计过程中，如果同时计算DNA的正反两条链，则根据碱基配对原则，A和T、C和G的出现频率相同。如果仅统计一条链，则虽然A和T、C和G的出现频率不同，但是非常接近。核苷酸核苷酸 频率频率 A0.344C0.155G0.157T0.343单链核苷酸出现频率单链核苷酸出现频率 基因和其它功能区域在正反两条链上出现的可能性通常一样 核苷酸出现频率也不应该有偏差 正反两条链在信息的组织结构方面不应该有差别 单链上A和T、C和G的出现频率相近。正反两条链碱基互补的原则 单链上A和T、C和G的出现频率相近的解释两联核苷酸频率
7、 不同基因组中两个连续核苷酸出现的频率也是不相同的4种核苷酸可以组合成16种两联核苷酸酵母基因组两联核苷酸频率表酵母基因组两联核苷酸频率表对酵母基因组对酵母基因组两联核苷酸的两联核苷酸的统计结果统计结果其中核苷酸对其中核苷酸对出现频率最高出现频率最高的达到的达到0.119而出现频率最而出现频率最低的只有低的只有0.028令: Pij 代表两联核苷酸（i，j）的出现频率 Pi 代表核苷酸i的出现频率 则 ： Pij= Pij/(PiPj) 的值反应核苷酸i和j的关联关系 如果Pij=1，则在两个连续的位置上，核苷酸i和j的出现是相对独立的。关联性分析关联性分析 对于酵母基因组 PA=0.3248
8、 PAA=0.1193 PAA =0.1193/（0.3248*0.3248） =1.131 1 表明在两个连续位置上“A”的出现不是独立的，而是相关的。关联性分析关联性分析 同样，对于相隔一定距离k（k代表核苷酸个数）的两个核苷酸，也可能具有一定的相关性。假设Pij(k)代表核苷酸j出现在核苷酸i之后第k个位置的频率，则可定义一个反应统计相关性的互信息I(k) I(k)值得大小实际上反应了距离为k的两个核苷酸之间的相关性的程度三联核苷酸基因密码子 在进行编码区域识别时，常常需要对三联核苷酸进行统计分析，这实际上是分析密码子的使用偏性。由于密码子的简并性（degeneracy），每个氨基酸至少
9、对应1种密码子，最多有6种对应的密码子。在基因中，同义密码子的使用并不是完全一致的。不同物种、不同生物体的基因密码子使用存在着很大的差异基因密码子的使用与基因编码的蛋白的结构和功能有关，与基因表达的生理功能有着密切的联系 蛋白的三级结构与密码子使用概率有密切的关系通过对密码子的聚类分析，可以很清晰地将具有不同三级结构蛋白质的编码基因分成不同的类，而具有相似三级结构蛋白的编码基因则大致聚在同一类中，从而证明基因密码子的使用偏性与蛋白质三级结构具有密切的相关性。 在不同物种中，类型相同的基因具有相近的同义密码子使用偏性对于同一类型的基因由物种引起的同义密码子使用偏性的差异较小 针对酵母第一染色体的
10、分析结果针对酵母第一染色体的分析结果第第3节节 功能位点分析功能位点分析功能位点（functional site）与特定功能相关的位点，是生物分子序列上的一个功能单元，或者是生物分子序列上一个较短的片段。功能位点又称为功能序列（functional sequence）、序列模式（motif）、信号（signal）等。核酸序列中的功能位点包括转录因子结合位点、转录剪切位点、翻译起始位点等。在蛋白质序列分析中，常使用序列模式这个名词，蛋白质的序列模式往往与蛋白质结构域或者作用部位有关。功能位点示意功能位点示意基因组序列中若干个相邻的功能位点组合形成功能区域（functional region）。功
11、能位点分析的任务发现功能位点特征识别功能位点1、利用共有序列搜索功能位点共有序列（consensus）又称一致性片段共有序列是关于功能位点特征的描述，它描述了功能位点每个位置上核苷酸进化的保守性 例如: NTATN 利用共有序列进行功能位点分析牵涉到两个方面的问题，如何构造共有序列如何利用共有序列在给定的核酸序列上搜索寻找功能位点，并计算所找到的功能位点的可靠性共有序列具有以下几个方面的特征：（1）共有序列中既有保守的位置，也有可变的位置；（2）任何位置上的核苷酸可以用15种类型之一来表示：核苷酸表示符号核苷酸表示符号符符 号号含含 义说 明明GG腺腺嘌嘌呤呤AA鸟嘌嘌呤呤TT胸腺胸腺嘧啶CC
12、胞胞嘧啶RG or A嘌嘌呤呤YT or C嘧啶MA or C氨基氨基KG or T羧基基SG or C强强氢键(3个个氢键)WA or T弱弱氢键(2个个氢键)HA or C or T非非GBG or T or C非非AVG or C or A非非T(非非U)DG or A or T非非CNG or A or T or C任意碱基任意碱基共有序列构造过程: (1)初始化共有序列为一系列可变位置，以“N”代表；(2)在可变位置寻找出现次数最多的核苷酸，并将该位置转化为保守位置；(3)对当前所得到的共有序列进行特异性检查，若通过检查，转（5），否则转（4）；(4)形成与当前共有序列一致的位点子集
13、，转（2）；(5)从原位点集合中删除与当前共有序列一致的位点，若还有剩余位点，则转（1），构造另外的共有序列。TTATGATATATACGCTTGTC TCCAC TTATGATATATACGCTTGTC TCCAC TNNNN tTATG tACGC tTGTC tCCAC tTATG tACGC tTGTC tCCAC TNNNC 1 2 3 4 2 3 NNNNNTNN NN非特异非特异 TNNNC非特异非特异 tACGc tTGTc tCCAc 4 2 tACGc tTGTc tCCAc 3 TNSNC 特异特异 5 Consensus1： TNSNC 剩余位点：剩余位点： TTATG
14、 ATATA 5 Consensus2： NTATN TN NSC在给定的序列中搜索与共有序列一致的序列片段数据库搜索共有序列表示方法的缺点：是关于序列特征的一种定性描述，对于DNA序列，它能够说明序列每个位置可能出现的碱基类型，但是不能准确地说明各位置上不同类型碱基出现的可能性大小。2、用感知矩阵分析功能位点用权系数描述功能位点各位置上每种核苷酸的相对重要性感知矩阵（或加权矩阵）根据一系列功能位点的多重对比排列结果而建立的其大小为4n 4代表碱基的种类数目，n代表功能位点的长度 矩阵的每一个元素M(a,j)的值代表第a种核苷酸在功能位点第j个位置上出现的得分,a A,T,G,C。123456
15、A18227-319T26142-10G3110-50-19C5-916880感知矩阵示例感知矩阵示例对于一个序列s=a1a2an，根据对应位置上核苷酸的类型，取感知矩阵中对应的权值，加和以后得到该序列的得分设S=ATTGCA，则 Ws = 1+6+14-5+8+19=43 T功能位点阈值T非功能位点阈值如果Ws T，则S是功能位点；如果WsT，则S是非功能位点。感知矩阵感知矩阵M的构造算法的构造算法令A+代表功能位点集合 A-代表非功能位点集合过程如下： （1）初始化M为零矩阵； （2）执行过程（3）-（6）的循环； （3）逐步取训练集合中的每个实例Si，如果Si A+，转 过程（4）；如果
16、Si A-，转过程（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矩阵能够正确识别训练集中的所有功能位点和非功能位点。对于最终得到的感知矩阵，要求其具有敏感性和特异性，每一列上的元素值应该尽可能地有明显的差别，以便反应功能位点各个位置上的特点。与感知矩阵类似，如
17、果令矩阵每一个元素M(a,j)的值代表第a种核苷酸在功能位点第j个位置上出现的概率，则M是一个概率矩阵。假设各个位置上出现的碱基是相互独立的，即任何两个位置上的碱基是不相关的，那么对于给定一个序列s=a1a2an，可以计算出功能位点序列为s的概率：如果分别统计功能位点和非功能位点，通过计算可以形成两个矩阵M和M，进一步计算可以判断一个给定的序列究竟属于功能位点，还是属于非功能位点。给定一个序列s=a1a2an，定义似然比LR(M,M,s)：在进行功能位点检测时，计算LR(M,M,s)，并与给定的阈值L比较，如果LR(M,M,s)L，则序列s可能是一个功能位点。概率矩阵M和M的每个元素是一个0和
18、1之间的正数。如果令一个4n新矩阵U的元素(a,j)的值为log2(M(a,j)/M(a,j)则矩阵U的每个元素值可能是正值，也可能是负值。 实际上，矩阵U就是感知矩阵。第第4节节 隐马尔柯夫模型隐马尔柯夫模型1、马尔柯夫链（Markov chain） 考虑一个具有多个状态的系统S，S=s1,s2,s|s|,令S0、S1、St为一系列在各个时刻系统状态的变量，即状态链。对于每个1到|S|的整数，它们分别与状态链中的一个状态相联系，并且在任何时刻，这条链都处于一个特殊的状态。当且仅当对于任何t有 则St形成一条马尔柯夫链。简单地说，就是系统未来的状态仅依赖于当前状态。St称为在时刻t系统链的状态
19、。一条马尔柯夫链完全决定于初始分布P(S0)和转换概率Pt=P(St+1|St)。 令状态转换矩阵为F =(fij) fij代表从状态si移动到状态sj的概率。生物序列可以被描述为一个随机过程的输出，其中对于一个给定的核酸在位置p出现的概率依赖于已占据前面k个位置的核酸，这样一种表示称为k阶马尔柯夫模型。ATCGTAGCAT.一个序列具有不同的统计性质（如二目频率或三目周期性）不同的功能区域（如编码区域、非编码区域）对应于不同的马尔柯夫模型。马尔柯夫链在识别CpG岛中的应用CpG岛是一类长度在几百bp的特殊DNA序列，其中CG核苷酸对出现的频率非常高。 ACGCGCGTACGCGAATCpG岛
20、在基因组中有重要的生物学意义，而识别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 岛内的字符转换概率 f
21、st-为CpG 岛外的字符转换概率则X的对数似然得分为上述计算值越大，则X越可能是CpG岛。 CpG岛内部和外部的转换概率岛内部和外部的转换概率 另外一个待解决的问题是： 给定DNA序列，确定CpG岛的位置。直接的方法：对窗口内的子序列计算得分Score(Xk)，具有正值的Xk 就是可能的CpG岛子序列起始位置为k+1 ，长度为l问题：事先不知道CpG岛的长度但是假设CpG岛的长度为l如果l比较大，而真实的CpG岛又比较小，则上述概率计算值不足以证实CpG岛；如果l取值比较小，则难以找出整个CpG岛。这是该算法的最大不足之处，需要考虑其他的算法。HMM2、隐马尔柯夫模型（HMM）功能位点的正则
22、表达式来表示相当于一致性序列这里的正则表达式描述了一个功能位点的构成规律，或者说描述了功能位点各个位置上核苷酸的组成。TGCCAGG ?ACACATC问题：对于每个位置，仅仅说明可能的取值，而没有说明各种取值出现的可能性大小例如，用这样的方法无法区分下面两条序列究竟哪一个更可能属于功能位点：TGCC-AGGACACATC第一个序列中，假设所有位置上都是取已知出现次数最少的字符而对于第二个序列，所有位置上都是取已知出现次数最多的字符。显然，第一个序列几乎肯定不是功能位点，而第二个序列几乎可以肯定是功能位点，但是用正则表达式表达却无法将两种极端的情况分开。隐马尔柯夫模型可以用于生物序列分析，该模型
23、在生物信息分析方面有重要的应用。一阶隐马尔柯夫模型包括有限数目的系统状态、离散的字母表、状态转换矩阵和字符释放概率。 一个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为
24、终止状态。例如，对于前面给出的两个序列例如，对于前面给出的两个序列ACACATC和和TGCTAGG，它们的得分分别为：，它们的得分分别为：P（ACACATC）= 0.8 1.0 0.8 1.0 0.8 0.6 0.4 0.6 1.0 1.0 0.8 1.0 0.8= 4.7 10-2P(TGCTAGG) = 0.2 1.0 0.2 1.0 0.2 0.6 0.2 0.6 1.0 1.0 0.2 1.0 0.2= 0.0023 10-2从上述计算结果可以看出，两个序列差别非常大从上述计算结果可以看出，两个序列差别非常大 一个功能位点的HMM模型是通过对一系列的功能位点实例进行机器学习而形成的用这
25、样的模型可以定量的计算一个序列片段是功能位点的可能性计算方法是从模型的第一个状态出发，根据序列的核苷酸组成，将相应的状态值与状态转换值连乘，结束于最后一个状态 一个检测CpG岛的HMM模型 有8个状态，状态名称和释放的字符为： 状态： A+ C+ G+ T+ A- C- G- T- 释放字符： A C G T A C G T 其中，带有“+”号的状态表示在CpG岛内部，用“-”号标记的状态代表CpG岛外部。假设字符处于CpG岛内的概率是p 处于CpG岛外的概率是q可以得到状态转换概率CpG岛岛HMM模型中的状态转换概率模型中的状态转换概率 解码问题：解码问题：给定一个隐马尔柯夫模型M =（A，
26、S，）和一个字符序列X，在M中为X寻找一条最优路径*，在路径中的每一个状态都选择释放一个字符，要求使得P(X|*)最大，记为：在处理 CpG岛问题中，最优路径可以帮助我们寻找CpG岛所在的位置。如果找到最优路径*，则这条路径穿过的“+”状态将对应于CpG岛。3、Viterbi算法求解HMM模型的最优路径基本思想：动态规划算法给定一个字符序列X=(x1,x2,xL) ，以vk(i)代表序列前缀（x1,x2,xi）终止于状态k（kS,1iL）的最可能路径的概率。求解过程如下：（1）初始化（2）对于每个i=0,L-1及每个lS，按下式进行递归计算：（3）最后，计算序列X终止于状态“end”最可能的路
27、径概率，即P(X|*)的值在正向的递归计算过程中，保持向前推进的反向指针，这样，在正向计算完成后，根据反向指针重构最优路径*。算法的时间复杂度为O(L|S|2),空间复杂度为O(L|S|)。在概率的计算过程中，需要使用大量的乘法运算，在有限计算精度的情况下，会产生误差。如果使用对数值，可以解决这个问题。因此，以vk(i)代表序列前缀（x1,x2,xi）终止于状态k（kS,1iL）的最可能路径的对数得分值，则初值按如下方式设置递归计算及最终得分计算改为（5-26）4、前向概率和反向概率给定一个隐马尔可夫模型M =（A，S，） 一个字符序列X=(x1,x2,xL)要求计算模型M产生X的概率P(X|
28、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模型中可能的路径非常多，穷举每条路径显然是不合适的。下面介绍解决该问题的前向算法（for
29、ward algorithm）与反向算法(backward algorithm)。算法的根本任务是对于每个1iL及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）
30、 利用正向和反向概率，可以计算出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） 特殊情况：假设已知与字符串序
31、列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）算法解决这个问题。具体的
32、求解算法如下： （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）反复执行步骤
33、（2）、（3），直到Score(X(1), X(2), X(n)|)的增量小于给定的一个值很小的参数为止。EM算法保证目标函数Score(X(1), X(2), X(n)|)单调增加，并且概率的对数值接近于0，保证算法收敛。需要注意，收敛的是目标函数，并非是的参数。当目标函数变化趋缓时，的参数值可能波动较大，这意味着算法所得到的结果不稳定。Baum算法的主要问题是目标函数存在若干局部极大，算法不能保证找到全局最大点，算法收敛的点可能是局部极大点。克服局部极大缺陷的一种方法是执行算法若干遍，每次给取不同的初始值。如果算法多次计算结果到达同一个极大点，则可以认为该点是全局最大点。6、基于HMM模型
34、的序列比对可以利用HMM将一个序列与一个序列统计特征（profile）进行比对，从而解决多重比对问题。定义一个长度为L的序列统计特征P是一系列的概率集合ei(b)，ei(b)表示在第i（1iL）个位置上出现字母表中字符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)。为了在比
35、对中允许插入“空位”的操作，在上述基本模型中加入“插入”状态I1,I2,IL，并假设每个插入状态Ij,有一个来自相应匹配状态Mj的连接，有一个到匹配状态Mj+1的连接，还有一个自循环连接。根据“空位”的惩罚原则，给这些状态转换赋予适当的概率。（5-50） 同样，为了允许“删除”操作，可以进一步假如“删除”状态D1,D2,DL，这些状态不能释放任何字符。删除状态依然顺序连接，同时增加从Dj到Ij的连接及从Ij到Dj+1的连接。完整的HMM模型如图5.5所示：D1D2D3I I2 2I I3 3I I4 4BeginBeginEndEndM M1 1M M2 2M M3 3I I1 1 图图5.5
36、 用于序列多重比对的用于序列多重比对的HMM模型模型 下面介绍一种Viterbi类似算法，将X =(x1,x2,xm)与长度等于L的统计特征P进行比对。对于每一个1jL和1im，定义： （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”的初始值为：为了
37、计算vjM(i)、vjI(i)和vjD(i) 的值，使用Viterbi算法中的相同技术，但现在的模型有两个特点：（1）模型中的每一个状态最多只有3个引入连接，如上图所示；（2）“删除”状态不释放任何字符。（5-51）“匹配”状态Mj的三个前驱同属于上一层,即j-1层，有“插入”状态Ij的三个前驱属于同一层,即j层，于是有（5-52）（5-53）“删除”状态Dj的三个前驱同属于上一层,即j-1层，并且由于Dj不释放任何字符，所以有最优得分计算公式如下： （5-54） （5-55） 前面我们介绍了如何将一个序列与一个以HMM模型代表的统计特征进行比对，那么怎样建立代表统计特征的HMM模型呢？即如何
38、确定模型中的各个状态、各状态的转换概率、各状态的字符释放概率？给定序列X =(x1,x2,xm)，定义前向概率为：反向概率为：前向概率的计算过程如下： （1）初始化 （2）递归计算（5-56）（5-57）反向概率的计算过程如下： （1）初始化 （2）递归计算（5-58）（5-59）利用HMM模型可以得到关于多重比对的近似结果。假设有n个序列S(1), S(2), S(n),如果代表统计特征的HMM模型已知，则将S(i)与比对，并将所得到的比对融入多重比对。如果还不知道代表统计特征的HMM模型，则按下述方法根据给定的序列构造：为选择合适的长度L，初始化的状态转换概率和各个状态的字符释放概率，利用
39、Baum算法训练模型，然后按照前面的方法根据进行多重序列的比对。可以进一步拓展上述多重序列比对的方法，形成在给定的多个序列中识别相似序列模式（pattern）的算法。 7、Gibbs采样 假设有n个序列S(1), S(2), S(n)及整数w，要求寻找这些序列的共同模式，即对于每一个序列S(i) ，找出一个最大长度等于w的子序列，使得这n个子序列的相似度最大。令a(1), a(2),a(n)分别为各子序列的起始下标，并令cij为字母表中第j个字符在这n个子序列第i位出现的次数，qij代表第j个字符在序列模式第i位出现的频率，pj代表第j个字符在所有序列中出现的频率。我们的目的是使得下面的得分最
40、大：（5-60） 具体的计算过程如下： （1）初始化：随机选择a(1), a(2),a(n)； （2）随机选择1zn，对于除S(z)之外的序列计算cij、qij和pj; （3）在序列S(z)中寻找最佳的子序列，并更新a(z),这里需要将S(z)与当前的模式统计特征进行局部比对； （4）重复执行第（2）步和第（3）步，直到得分值的变化小于预先设定的值为 应用Gibbs采样算法在许多情况下可以得到较好的效果。但是该算法也存在一些问题，如，算法可能会陷入局部极小，不能保证收敛到最佳的模式。又如，选择不同的长度值w，所得到的结果可能相差非常大。8、其他方法功能位点的其他分析方法有神经网络方法及相邻核苷
41、酸分析方法等。有一点需要说明，前面我们介绍的是单个功能位点的识别分析，其识别结果往往不太好。实际上在核酸序列中，各种功能位点共存，协同作用，形成功能域。因此，在进行功能位点识别时，最好同时考虑多个位点，考虑各种位点之间的关联，对它们进行同时识别，可能会得到更好的识别结果。这就是所谓功能域的识别。第第5节节 基因识别基因识别基因识别是生物信息学领域里的一个重基因识别是生物信息学领域里的一个重要研究内容要研究内容 基因识别问题，在近几年受到广泛的重基因识别问题，在近几年受到广泛的重视视 当人类基因组研究进入一个系统测序阶段当人类基因组研究进入一个系统测序阶段时，急需可靠自动的基因组序列翻译解释时，
42、急需可靠自动的基因组序列翻译解释技术，以处理大量已测定的但未知功能或技术，以处理大量已测定的但未知功能或未经注释的未经注释的DNA序列序列 原核基因识别原核基因识别重点在于识别编码区域重点在于识别编码区域非翻译区域（非翻译区域（untranslated regions, UTR） 编码区域两端的编码区域两端的DNA，有一部分被转录，有一部分被转录，但是不被翻译，这一部分称为非翻译区域但是不被翻译，这一部分称为非翻译区域 5UTR-基因上游区域的非翻译区域基因上游区域的非翻译区域 3UTR-基因下游区域的非翻译区域基因下游区域的非翻译区域对于任何给定的核酸序列（单链DNA或mRNA），根据密码子
43、的起始位置，可以按照三种方式进行解释。例如，序列ATTCGATCGCAA这三种阅读顺序称为阅读框（reading frames）CAA A ATTCGATCGATTCGATCGCAAATTCGATCGCA（1）（3）（2）一个开放阅读框（一个开放阅读框（ORF,open reading frame）是一个没有终止编码的密码子序是一个没有终止编码的密码子序列。列。原核基因识别任务的重点是识别开放阅读原核基因识别任务的重点是识别开放阅读框，或者说识别长的编码区域。框，或者说识别长的编码区域。基于基因密码子特性的识别方法基于基因密码子特性的识别方法辨别编码区域与非编码区域的一种方法辨别编码区域与非编
44、码区域的一种方法是检查终止密码子的出现频率是检查终止密码子的出现频率 终止密码子出现的期望次数为：终止密码子出现的期望次数为： 每每21个（个（ 64/3）密码子出现一次终止密）密码子出现一次终止密码子码子 基本思想：基本思想：如果能够找到一个比较长的序列，其相应如果能够找到一个比较长的序列，其相应的密码子序列不含终止密码子，则这段序的密码子序列不含终止密码子，则这段序列可能就是编码区域。列可能就是编码区域。基本算法：基本算法：扫描给定的扫描给定的DNA序列，在三个不同的阅读序列，在三个不同的阅读框中寻找较长的框中寻找较长的ORF。遇到终止密码子以遇到终止密码子以后，回头寻找起始密码子。后，回
45、头寻找起始密码子。这种算法过于简单，不适合于处理短的这种算法过于简单，不适合于处理短的ORF或者交叠的或者交叠的ORF。识别编码区域的另一种方法是分析各种识别编码区域的另一种方法是分析各种密码子出现的频率密码子出现的频率 将一个随机均匀分布的将一个随机均匀分布的DNA序列翻译成氨基酸序列翻译成氨基酸序列，则在氨基酸序列中上述序列，则在氨基酸序列中上述3种氨基酸出现的种氨基酸出现的比例应该为比例应该为6:4:1例如，亮氨酸、丙氨酸、色氨酸分别有例如，亮氨酸、丙氨酸、色氨酸分别有6个、个、4个和个和1个密码子个密码子但是在真实的氨基酸序列中，上述比例并不但是在真实的氨基酸序列中，上述比例并不正确正
46、确这说明这说明DNA的编码区域并非随机的编码区域并非随机假设在一条假设在一条DNA序列中已经找到所有的序列中已经找到所有的ORF，那么可以利用密码子频率进一步那么可以利用密码子频率进一步区分编码区分编码ORF和非编码和非编码ORF马尔柯夫链模型马尔柯夫链模型利用这种方法，可以计算一个利用这种方法，可以计算一个ORF成为成为编码区域的可能性。编码区域的可能性。一个简单的统计模型一个简单的统计模型假设相继的密码子是独立的，不存在前后依假设相继的密码子是独立的，不存在前后依赖关系。赖关系。 令令fabc代表密码子代表密码子abc在编码区域出现的频率在编码区域出现的频率给定序列给定序列a1,b1,c1
47、, a2,b2,c2, an+1,bn+1从密码子从密码子a1b1c1开始的阅读框，其开始的阅读框，其n个密码子个密码子的出现概率为的出现概率为第二种和第三种阅读框第二种和第三种阅读框n个密码子出现的概率个密码子出现的概率分别为分别为第第i个阅读框成为编码阅读框的概率个阅读框成为编码阅读框的概率计算：计算：算法：算法：在序列上移动长度为在序列上移动长度为n的窗口，计算的窗口，计算Pi根据根据Pi的值识别编码的阅读框的值识别编码的阅读框基于编码区域碱基组成特征的识别方法基于编码区域碱基组成特征的识别方法编码序列与非编码序列在碱基组成上有编码序列与非编码序列在碱基组成上有区别区别单个碱基的组成比例
48、单个碱基的组成比例多个碱基的组成多个碱基的组成通过统计分析识别编码序列通过统计分析识别编码序列分析实例分析实例2、真核基因识别问题、真核基因识别问题 真核基因远比原核基因复杂：真核基因远比原核基因复杂：一方面，真核基因的编码区域是非连续一方面，真核基因的编码区域是非连续的，编码区域被分割为若干个小片段。的，编码区域被分割为若干个小片段。另一方面，真核基因具有更加丰富的基另一方面，真核基因具有更加丰富的基因调控信息，这些信息主要分布在基因因调控信息，这些信息主要分布在基因上游区域。上游区域。 基因识别基本思路基因识别基本思路 找出基因两端的功能区域找出基因两端的功能区域:转录启动区转录启动区 终
49、止区终止区 在启动区下游位置寻找翻译起始密码子在启动区下游位置寻找翻译起始密码子 识别转录剪切位点识别转录剪切位点剪切给体位点剪切给体位点剪切接受体位点剪切接受体位点各种不同的方法有不同的适应面，而不各种不同的方法有不同的适应面，而不同的方法有时可以结合起来以提高基因同的方法有时可以结合起来以提高基因识别的准确率。识别的准确率。关键问题是如何提高一个识别算法的敏关键问题是如何提高一个识别算法的敏感性（感性（sensitivity，Sn）和特异性和特异性（specificity，Sp）。）。 3、基因识别的主要方法、基因识别的主要方法两大类识别方法：两大类识别方法：从头算方法（或基于统计的方法）
50、从头算方法（或基于统计的方法）根据蛋白质编码基因的一般性质和特征进行识别，根据蛋白质编码基因的一般性质和特征进行识别，通过统计值区分外显子、内含子及基因间区域通过统计值区分外显子、内含子及基因间区域 基于同源序列比较的方法基于同源序列比较的方法利用数据库中现有与基因有关的信息（如利用数据库中现有与基因有关的信息（如EST序序列、蛋白质序列），通过同源比较，帮助发现新列、蛋白质序列），通过同源比较，帮助发现新基因。基因。最理想的方法是综合两大类方法的优点，最理想的方法是综合两大类方法的优点，开发混合算法。开发混合算法。基因识别方法有基因识别方法有 ： （1）基于规则的系统）基于规则的系统 （2）
51、 语义学方法语义学方法 （3） 线性辨别分析（线性辨别分析（LDA） （4） 决策树决策树 （5） 动态规划动态规划 （6） 隐马尔柯夫模型隐马尔柯夫模型 （7） 剪切对比排列剪切对比排列 （spliced alignment）4、编码区域识别、编码区域识别两类方法两类方法 ：基于特征信号的识别基于特征信号的识别 内部外显子内部外显子剪切位点剪切位点5端的外显子一定在核心启动子的下游端的外显子一定在核心启动子的下游3端的外显子的下游包含多聚信号和终端的外显子的下游包含多聚信号和终止编码止编码 基于统计度量的方法基于统计度量的方法 根据密码子使用倾向根据密码子使用倾向双联密码统计度量等双联密码统
52、计度量等在一个基因中，第在一个基因中，第i个（个（i=1，64）密码子相对使用密码子相对使用倾向倾向RSCUi的定义如下：的定义如下：Obsi是该基因中第是该基因中第i个密码子实际出现的次数个密码子实际出现的次数Expi是对应密码子期望的出现次数是对应密码子期望的出现次数 aai是统计的第是统计的第i个密码子出现的次数个密码子出现的次数 syni是所有与第是所有与第i个密码子同义密码子出现的次数个密码子同义密码子出现的次数RSCU大于大于1表示相应密码子出现的次数比期望次数表示相应密码子出现的次数比期望次数高，而小于高，而小于1则表示出现次数相对较少。则表示出现次数相对较少。（5-66） （5
53、-65） 密码子使用倾向密码子使用倾向设一段DNA序列为S，从S的第i位到第j位的双联密码统计度量IF6（i，j）定义为： fk是从第k位开始的双联密码的频率Fk是该双联密码随机出现的频率（5-67）双联密码统计度量双联密码统计度量通过相似搜索发现编码区域或者外显子通过相似搜索发现编码区域或者外显子 EST（Expressed Sequence Tags）cDNA 蛋白质序列蛋白质序列目前大多数预测程序都将数据库相似性目前大多数预测程序都将数据库相似性搜索的信息结合进基因预测过程搜索的信息结合进基因预测过程同时考虑序列特征信号和统计度量同时考虑序列特征信号和统计度量GRAIL用人工神经网络识别
54、编码区域用人工神经网络识别编码区域输入是一系列反映功能位点信号特征和序列编码统计特征的参数输入是一系列反映功能位点信号特征和序列编码统计特征的参数输出就是对一段输出就是对一段DNA序列是否是编码区域的判别结果序列是否是编码区域的判别结果神经网络具有非线性映射能力，能够发现输入和输出之间的高阶相关性神经网络具有非线性映射能力，能够发现输入和输出之间的高阶相关性5、构建基因模型基因识别最终任务是建立完整的基因结构模型一个理想的基因识别程序应该能够发现完整的基因结构 （,e1, i1, , in-1, en , ） ATG-外显子1内含子外显子外显子n-UAG基因剪切位点基因剪切位点 剪切给体（剪切
55、给体（donor）位点位点- “gt” 接受体（接受体（acceptor）位点位点- “ag”基因的可变剪切基因的可变剪切gene A基因可变剪切示意基因可变剪切示意构建基因模型方法构建基因模型方法 剪切位点形成外显子和内含子的边界剪切位点形成外显子和内含子的边界 搜集候选外显子搜集候选外显子 候候选基因基因候选基因是一条候选基因是一条非相交非相交的外显子和内含的外显子和内含子的链，表示为子的链，表示为 （i0, e1, i1, , en, in） 其中其中ij代表内含子（代表内含子（0 j n） el代表外显子（代表外显子（1 l n） i0和和in并非真实的内含子，它们分别代表并非真实的内
56、含子，它们分别代表基因两侧的非编码序列基因两侧的非编码序列候选基因位于给定的候选基因位于给定的DNA序列，并满足下列一序列，并满足下列一致性条件：致性条件： （1）所有外显子加起来的长度是）所有外显子加起来的长度是3的整数倍；的整数倍； （2）在各个外显子内部（除最后一个外显子的）在各个外显子内部（除最后一个外显子的最后一个密码子），没有终止编码；最后一个密码子），没有终止编码； （3）第一个内含子）第一个内含子-外显子边界（外显子边界（i0, e1）是翻是翻译起始编码，而最后一个外显子译起始编码，而最后一个外显子-内含子边界内含子边界（en, in）是终止编码。是终止编码。位点图位点图（分层
57、标注剪切位点）（分层标注剪切位点）另设两个特殊的顶点，即起点（另设两个特殊的顶点，即起点（source）和终点（和终点（sink）。从起点到终点的任何一条路径代表一个可能的基因结构。从起点到终点的任何一条路径代表一个可能的基因结构。例如例如： 位点图上的路径位点图上的路径候选基因所对应的道路图中的路径候选基因所对应的道路图中的路径求最优路径求最优路径每一条弧附加一个权值每一条弧附加一个权值外显子、内含子度量每个节点附加权值每个节点附加权值剪切位点度量 综合评价综合评价 6、用于基因识别的HMM模型隐马尔柯夫模型隐马尔柯夫模型HMM是一条状态不可见是一条状态不可见的马尔柯夫链，其当前状态的输出是
58、可见的马尔柯夫链，其当前状态的输出是可见的。的。每个状态按照一定的概率分布随机地从字母每个状态按照一定的概率分布随机地从字母表中取出字符并释放。表中取出字符并释放。 扩展的隐藏马尔柯夫模型（扩展的隐藏马尔柯夫模型（GHMMs ）对对HMM进一步抽象，产生更一般的马尔柯夫进一步抽象，产生更一般的马尔柯夫模型，以分析复杂的脊椎动物基因模型，以分析复杂的脊椎动物基因 。(1) 信号传感器模型信号传感器模型将剪切位点、起始编码区域或者终止编码区域看成是DNA序列上的功能位点或者信号位点，用HMM来进行分析 内含子区域内含子区域 外显子区域外显子区域 保守位点保守位点 根据对比排列，根据对比排列，形成具
59、有形成具有19状态状态的的HMM模型。模型。 对前一节所介绍的对前一节所介绍的HMM模型进行修改，模型进行修改，可以处理双联核苷酸的问题，即将可以处理双联核苷酸的问题，即将4种种概率分布扩展为概率分布扩展为16种。种。假设一段序列为假设一段序列为ACTGTC，则，则 P(ACTGTC)=p1(A) p2(C A) p3(T C) p4(G T) p5(T G) p6(C T) 其中其中p1是状态是状态1对于对于4种核苷酸的概率，种核苷酸的概率，p2(x y)状态状态2的条件概率。的条件概率。(2) 编码区模型由于密码子的长度为3，因此密码子模型的最后一个状态应该至少为2阶。对于2阶的状态，具有
60、64种概率分布，可根据已知编码区域进行统计计算而得到64种分布。例如：p(ACA)=c(CAA)/c(CAA)+c(CAC)+c(CAG)+c(CAT) p(CCA)=c(CAC)/c(CAA)+c(CAC)+c(CAG)+c(CAT) p(GCA)=c(CAG)/c(CAA)+c(CAC)+c(CAG)+c(CAT) p(TCA)=c(CAT)/c(CAA)+c(CAC)+c(CAG)+c(CAT) 其中，c(xyz)是密码子xyz的计数。这样的模型可以检测无结束编码的区域，因为对应于三个结束编码TAA、TAG和TGA的p(ATA)、p(GTA)和p(ATG)自动为0。(3) 组合模型将上述
61、模型扩展，使之可以识别具有多个外显子的基因。 改进后的模型见下图、基于剪切比对的基因识别方法基本思想是：利用数据库中的同源信息进行基因识别，包括DNA、RNA和蛋白质数据库。 其方法是：首先通过分析所有可能的剪切接受体位点和剪切给体位点，构建一组候选的外显子。然后进一步分析候选外显子，探查所有可能的外显子组合，寻找一个与已知目标蛋白质或其他表达序列最匹配的组合 一种半自动的综合方法识别基因过程：（1）选择所有长度大于50bp并介于保守的剪切接受位点和给体位点之间的ORF，作为候选的外显子；预选预选（2）对于候选的外显子计算其6目编码度量值，并从大到小将它们排列起来；减小搜索范围减小搜索范围（3
62、）对照蛋白质序列数据库进行搜索，寻找相似体。搜索，筛选搜索，筛选、基因识别程序介绍表表5.7 基因识别程序及访问地址基因识别程序及访问地址（HP主页；主页；ESE-mail服务器；服务器；WSweb服务器；服务器；CL客户客户/服务器协议；服务器协议；EX有可执行代码；有可执行代码；SC有源代码）有源代码）表表5.8 各程序的性能比较（敏感性各程序的性能比较（敏感性(1)被预测出的真实编码核酸的被预测出的真实编码核酸的%；敏感性敏感性(2)被正确识别出的编码外显子的被正确识别出的编码外显子的%；特异性特异性(1)预测出的编码核酸为真实编码核酸的预测出的编码核酸为真实编码核酸的%；特异性特异性(
63、2)预测出外显子为真实外显子的预测出外显子为真实外显子的%）基因识别方法存在的问题和局限性：基因识别方法存在的问题和局限性： （1）关于基因的定义不明确）关于基因的定义不明确 统一定一定义 （2）目前的方法仅仅识别蛋白质编码基因）目前的方法仅仅识别蛋白质编码基因 转录信号信号 （3）现有的许多方法仅检测单个基因）现有的许多方法仅检测单个基因部分基因、多重基因部分基因、多重基因 （4）基于同源分析的方法是保守的）基于同源分析的方法是保守的不可能发现新的基因不可能发现新的基因 （5） 忽视关于基因结构的生物学知识忽视关于基因结构的生物学知识 基因表达的真实分子机制基因表达的真实分子机制 目标基因的
64、分析及应用目标基因的分析及应用 目标基因的分析及应用数据公布，形成文章cDNA文库的构建随机挑取克隆5或3端测序序列聚类和拼接功能分类及代谢分析EST数据注释分析平台的构建文献检索与数据收集cDNAcDNA序列序列序列序列基因组序列基因组序列基因组序列基因组序列蛋白质序列蛋白质序列翻译翻译Codon biasGC Content酶切位点酶切位点引物设计引物设计编码区预测编码区预测基因结构分析基因结构分析选择性剪切选择性剪切SNP序列比对序列比对功能注释功能注释KEGGGO系统发育树系统发育树蛋白质理化性质蛋白质理化性质二级结构预测二级结构预测结构域分析结构域分析重要信号位点分析重要信号位点分析
65、三级结构预测三级结构预测核苷酸序列分析核苷酸序列分析基因编码区组分分析基因编码区组分分析GC含量含量/Codon bias引物设计引物设计限制性核酸内切酶位点预测限制性核酸内切酶位点预测基因编码区结构分析基因编码区结构分析基因结构分析基因结构分析选择性剪切分析选择性剪切分析/SNP分析分析基因调控区域分析基因调控区域分析蛋白质序列分析蛋白质序列分析蛋白质一级序列蛋白质一级序列蛋白质理化性质分析蛋白质理化性质分析蛋白质二级结构蛋白质二级结构蛋白质二级结构预测蛋白质二级结构预测蛋白质序列信号位点分析蛋白质序列信号位点分析蛋白质超二级结构蛋白质超二级结构蛋白质结构域分析蛋白质结构域分析蛋白质三级结构
66、蛋白质三级结构蛋白质三维结构模拟蛋白质三维结构模拟序列比对序列比对序列比对注释序列比对注释多序列比对多序列比对系统发育分析系统发育分析系统发育分析系统发育分析核苷酸序列分析核苷酸序列分析内容内容分析软件分析软件基因编码区分析基因编码区分析Spidey 选择性剪切分析选择性剪切分析根据根据Spidey分析结果分析结果蛋白质序列分析蛋白质序列分析实验内容实验内容分析软件分析软件蛋白质理化性质分析蛋白质理化性质分析ProtParam/ProtScale蛋白质二级结构分析蛋白质二级结构分析TMpred蛋白质三维结构分析蛋白质三维结构分析SWISS-MODEL /SWISS-PdbView系统发育分析系
67、统发育分析实验内容实验内容分析软件分析软件系统发育分析系统发育分析PHYLIP/TreeView实例实例 核苷酸序列分析核苷酸序列分析核苷酸序列分析开放读码框（open reading frame, ORF)的识别基因结构分析内含子/外显子剪切位点识别选择性剪切分析CpG 岛的识别核心启动子/转录因子结合位点/转录启始位点的识别转录终止信号的预测GC含量/密码子偏好性分析核苷酸序列分析开放读码框的识别开放读码框的识别开放读码框（open reading frame, ORF) 是一段起始密码子和终止密码子之间的碱基序列ORF 是潜在的蛋白质编码区基因组DNA序列cDNA，mRNA，EST核苷酸
68、序列分析ORFORF基因开放阅读框基因开放阅读框/ /基因结构分析识别工具基因结构分析识别工具Getorfhttp:/saturn.nchc.gov.tw:9091/Pise/getorf.htmlWeb/LinuxPlotorfhttp:/saturn.nchc.gov.tw:9091/Pise/plotorf.htmlWeb/LinuxORF Finder http:/www.ncbi.nlm.nih.gov/gorf/gorf.html WebBestORFhttp:/ http:/opal.biology.gatech.edu/GeneMark/WebGene Finderhttp:/
69、rulai.cshl.org/tools/genefinder/（Dr. Michael Zhang ）WebFGENESHhttp:/ LinuxFgeneSB/ FgeneSVhttp:/ http:/compbio.ornl.gov/generation/WebGeneBuilder http:/r.it/webgene/genebuilder.html WebFGENESH+ /+http:/ Web/LinuxGenomeScan http:/genes.mit.edu/genomescan.html WebGeneWise http:/www.sanger.ac.uk/Softwa
70、re/Wise2/ WebGRAILhttp:/grail.lsd.ornl.gov/grailexp/Web/Linux/WindowsBCM Gene Finderhttp:/searchlauncher.bcm.tmc.edu/seq-search/gene-search.htmlWeb核苷酸序列分析ORFORF内含子内含子/外显子剪切位点识别外显子剪切位点识别对基因组序列的读码框区域进行预测内含子5端供体位点(donor splice site): GT内含子3端受体位点(acceptor splice site): AG预测工具：GENSCAN，GENEMARKNetGene2, S
71、plice View核苷酸序列分析Gene StructureGene Structure内含子内含子/外显子剪切位点识别外显子剪切位点识别如何分析mRNA/cDNA的外显子组成？RNASPL与相应的基因组序列比对，分析比对片段的分布位置预测工具：Spidey，SIM4，BLAT，BLAST，FASTA 核苷酸序列分析Gene StructureGene Structure基因开放阅读框基因开放阅读框/基因结构分析工具基因结构分析工具对基因组序列的读码框区域进行预测对基因组序列的读码框区域进行预测NNSplicehttp:/www.fruitfly.org/seq_tools/splice.h
72、tmlWeb Splice Viewhttp:/r.it/webgene/wwwspliceview.htmlWeb NetGene2http:/www.cbs.dtu.dk/services/NetGene2/WebSPL/SPLM/RNASPL/FSPLICEhttp:/ Web/LinuxSpideyhttp:/www.ncbi.nih.gov/spideyWebPROT_MAPhttp:/ Sim4http:/gamay.univ-perp.fr/analyse_seq/sim4http:/globin.cse.psu.edu/Web/LinuxBLAThttp:/www.cse.uc
73、sc.edu/kent/src/unzipped/blat/LinuxBLASTftp:/ftp.ncbi.nlm.nih.gov/BLAST/ExecutablesWeb/Windows/LinuxFASTAftp:/ftp.virginia.edu/pub/fasta/win32_fasta/fasta34t21b5d.zipWeb/Windows/Linux核苷酸序列分析Gene StructureGene Structure选择性剪切选择性剪切（Alternative splicing）分分析析选择性剪接是调控基因表达的重要机制了解不同物种、细胞、发育阶段、环境压力下基因的调控表达机制
74、分析方法：查询选择性剪切相关的网站多序列比对核苷酸序列分析Gene StructureGene Structure查询选择性剪切相关的网站查询选择性剪切相关的网站http:/www.ebi.ac.uk/asd/index.html0综合综合http:/splicenest.molgen.mpg.de/综合综合http:/cgsigma.cshl.org/new_alt_exon_db2/综合综合http:/166.111.30.65/AsMamDB/哺乳动物哺乳动物http:/www.tigr.org/tigr-scripts/tgi/splnotes.pl?species=humanhttp
75、:/prosplicer.mbc.nctu.edu.tw/http:/www.bit.uq.edu.au/altExtron人人http:/www.cse.ucsec.edu/kent/intronerator/altsplice.html线虫线虫http:/ StructureGene Structure 从已知基因的功能推测剪切机制从已知基因的功能推测剪切机制http:/prosplicer.mbc.nctu.edu.tw/查询查询NOX1核苷酸序列分析Gene StructureGene Structure 基于序列比对分析选择性剪切基于序列比对分析选择性剪切在序列上高度相似的在序列上高
76、度相似的mRNA/cDNA/EST序列序列相匹配的基因组序列相匹配的基因组序列序列比对序列比对对分布位置进行分析对分布位置进行分析cDNA/mRNA/EST 序列比对序列比对收集序列收集序列核苷酸序列分析Gene StructureGene Structure 评判的标准：评判的标准：来自来自Unigene的高质量数据的高质量数据Exon至少有至少有3条条ESTs覆盖覆盖Exon周围有周围有GT-AG信号信号Blast比对比对Score值值100相似度相似度95%S.Gupta et al., Genome wide identification and classification of a
77、lternative splicing based on EST data, 2004, 20(16): 2579-2585基因周围调控序列分析基因周围调控序列分析CpG岛位于真核生物基因转录起始位点上游，GC含50% ，长度200bp转录起始位点(Transcription start site, TSS)PY2CAPY5核心启动子(Core promoter element)TATA box，Pribnow box 上游启动子元件(Upstream promoter element)CAAT box，GC box，SP1，Otc转录终止信号AAUAAA，UUUUUU操纵子、终止子、增强子、
78、沉默子核苷酸序列分析Gene StructureGene Structure 启动子数据库启动子数据库TransFac http:/www.gene- EPD http:/www.epd.isb-sib.ch/ TRRD http:/wwwmgs.bionet.nsc.ru/mgs/gnw/trrd Jasparhttp:/jaspar.cgb.ki.se/cgi-bin/jaspar_db.plZhang Labhttp:/rulai.cshl.org/software/index1.htm DBTSShttp:/dbtss.hgc.jp/index.htmlMIRAGEhttp:/www.
79、ifti.org/ Bacillus subtilis http:/dbtbs.hgc.jp/ Drosophila melanogaster http:/www-biology.ucsd.edu/labs/Kadonaga/DCPD.html E. coli http:/arep.med.harvard.edu/ecoli_matrices/ Human http:/zlab.bu.edu/mfrith/HPD.html PlantProm http:/ =plantprom&group=data&subgroup=plantprom Planthttp:/www.dna.affrc.go.
80、jp/PLACE http:/oberon.fvms.ugent.be:8080/PlantCARE/index.html Saccharomyces cerevisiae http:/cgsigma.cshl.org/jian/ 核苷酸序列分析Gene StructureGene Structure CpG Island 分析分析CpG Island http:/www2.ebi.ac.uk/cpg/ WebCpG finderhttp:/ Network Promoter Predictionhttp:/www.fruitfly.org/seq_tools/promoter.htmlWeb
81、Softberry: BPROM, TSSP, TSSG, TSSWhttp:/ =index&group=programs&subgroup=promoterWebMatInspectorhttp:/www.gene-regulation.de/WebTRANSPLORERhttp:/www.biobase.de/pages/products/transplorer.htmlWebRSAThttp:/rsat.ulb.ac.be/rsat/WebSIGSCANhttp:/bimas.dcrt.nih.gov/molbio/signal/Web/WindowsConsInspector ftp
82、:/ariane.gsf.de/pub/win95_NT/ConsInspector.exeWindows/LinuxCisterhttp:/zlab.bu.edu/mfrith/cister.shtmlWeb 核苷酸序列分析Gene StructureGene Structure 转录终止信号预测转录终止信号预测Hcpolyahttp:/r.it/webgene/wwwHC_polya.htmlWeb POLYAHhttp:/ =polyah&group=programs&subgroup=promoterWebpolyadq http:/rulai.cshl.org/tools/polya
83、dq/polyadq_form.htmlWeb核苷酸序列分析Gene StructureGene Structure 编码区综合分析编码区综合分析核苷酸序列分析Gene StructureGene StructureCpG岛分析No调控序列所在位置Cister结果：881-896CCAAT908-923CCAAT转录终止信号GetOrfGenScan735-773964-10201054-11461112-11561341-16251054-1490(1054-1145,1268-1490)CCTAGTCCAGACGCCATGGGT比对分析（Blastx，Blastn，Blastp）Blast
84、x结果： gammaglobin：1054-11461266-1493http:/zlab.bu.edu/mfrith/HPD.htmlHuman gene 5HSA004013: -10001000GeneBuilderhttp:/r.it/webgene/genebuilder.html核苷酸序列分析ToolsTools核苷酸序列综合分析软件GeneBuilderhttp:/r.it/webgene/genebuilder.htmlDNA Toolhttp:/www.crc.dk/dnatools/downloads/setup/dt6_setup.exeSEQtoolshttp:/www
85、.seqtools.dk/DNAssisthttp:/www.dnassist.org/dnassist20.zipGeneTool http:/ DNAmanhttp:/ Striderhttp:/ pDRAW32http:/ NTIhttp:/ Sequence assembly Sequence manipulation Homology comparison Multiple alignmentGene structure analysis Primer/Oligo analysis Restriction analysis Codons analysis核苷酸序列分析ToolsToo
86、ls实例 一步骤一：基因结构分析使用工具：Spidey步骤二：选择性剪切分析使用Spidey的分析结果数据：序列Seq1来自拼接结果序列Seq2Seq6来自Seq1与nr数据库BLASTN比对结果（Blosum62, E value=0.001）比对得到来自人类的complete cds 9条，去除冗余的序列后剩下5条Accession number of Seq2Seq6: AF127763，AF166326，AF166327, AF166328, BC075014核苷酸序列分析 practicepracticeSpideyNCBI开发的在线预测程序http:/www.ncbi.nih.go
87、v/spidey基于BLAST和Dot View局部联配的算法 核苷酸序列分析 practicepractice序列在线提交形式：界面中有两个窗口：上方窗口用于输入基因组序列（直接粘贴序列或用Genbank ID/AC号）下方窗口用于输入cDNA/mRNA序列（直接粘贴序列或用Genbank ID/AC号）可同时输入多条cDNA/mRNA序列与同一条基因组序列进行分析主要选项主要选项/参数参数输入基因组序列输入基因组序列Z83819输入输入lesson7.seq文档中的文档中的5条序列条序列判断用于分析的序列间的差异，判断用于分析的序列间的差异，并调整比对参数并调整比对参数设定内含子的长度，设
88、定内含子的长度，默认长度：内部内含子默认长度：内部内含子为为35kb, 末端内含子为末端内含子为100kb比对阈值比对阈值选择物种选择物种输出格式输出格式核苷酸序列分析 practicepractice输出结果输出结果第一条序列为基因组序列第一条序列为基因组序列浅蓝色为浅蓝色为mRNA/cDNA序列序列橘黄色为外显子橘黄色为外显子外显子对应于外显子对应于基因组上的基因组上的起始起始/ /结束位置结束位置外显子对应于外显子对应于mRNA/cDNA上的上的起始起始/结束位置结束位置外显子外显子长度长度一致性一致性百分比百分比错配和错配和gap序列联配结果序列联配结果外显子外显子序号序号核苷酸序列分析 practicepracticeNox基因AF127763，AF166326，AF166327和AF127763 核苷酸序列分析 practicepractice
