DNA序列分析东南大学生物信息学实验室文档格式.docx
- 文档编号:18595408
- 上传时间:2022-12-29
- 格式:DOCX
- 页数:26
- 大小:464.88KB
DNA序列分析东南大学生物信息学实验室文档格式.docx
《DNA序列分析东南大学生物信息学实验室文档格式.docx》由会员分享,可在线阅读,更多相关《DNA序列分析东南大学生物信息学实验室文档格式.docx(26页珍藏版)》请在冰豆网上搜索。
功能序列分析的准确性来自于对“功能序列”和“非功能序列”的辨别能力。
在实际应用中,解决具体识别问题的过程是一个多阶段的过程。
首先收集已知的功能序列和非功能序列实例,并且要求这些序列之间是非相关的。
将这些序列混合在一起,形成两个集合。
一个集合是训练集(trainingset),用于建立完成识别任务的数学模型。
另一个集合是测试集或控制集(controlset),用于检验所建模型的正确性。
用训练集中实例对预测模型进行训练,使之通过学习后具有正确处理和辨别能力。
然后,用模型对测试集中的实例进行“功能”与“非功能”的判断,根据判断结果计算模型识别的准确性。
设Tp是程序正确识别的功能序列数,Tn为正确识别的非功能序列数,Fn是被错误识别为非功能序列的功能序列数,Fp是被错误识别为功能序列的非功能序列数。
以Sn和Sp分别代表识别程序对功能序列的敏感性和识别的特异性,其计算公式如下(Baldietal.,2000):
(5-1)
(5-2)
对于一个实用程序,既要求有较高的敏感性,也要求有较高的特异性。
如果敏感性很高,但特异性比较低,则在实际应用中会产生高比率的假阳性;
相反,如果特异性很高,而敏感性比较低,则会产生高比率的假阴性。
对于敏感性和特异性需要进行权衡,给出综合评价指标。
对于一个识别程序准确性可按下式进行综合评价:
(5-3)
另一个综合评介指标为相关系数,其计算计算公式为:
(5-4)
Baldi等人对生物信息学中分类算法的预测准确性评价进行了深入的研究,详细结果参见文献(Baldietal.,2000)
在检测算法的可行性时,需要从已知的数据中按照不同的方式选择训练集和测试集,反复进行训练-测试,经过多次实验,得到关于该方法可行性的统计结果。
在评判一个识别程序或比较各个识别程序时,测试集的构成非常关键,在不同的测试集上进行测试可能会得到不同的准确性结果,甚至准确性相差很大。
因此,为了公正客观地评价功能序列识别,必须建立标准的功能序列测试集合。
目前已有一些这样的标准集合,如基因转录剪切位点的测试集合、编码区域的测试集合等。
第二节核苷酸关联分析
对于一个给定的基因组,最简单的计算就是统计DNA序列中各类核苷酸出现的频率。
对于随机分布的DNA序列,每种核苷酸的出现是均匀分布的,出现频率各为0.25。
而真实基因组的核苷酸分布则是非均匀的,如表5.1中酵母基因组核苷酸出现频率。
表5.1酵母基因组核苷酸出现频率
核苷酸
频率
A
0.3248693727808
C
0.175********92
G
T
在统计过程中,如果同时计算DNA的正反两条链,则根据碱基配对原则,A和T、C和G的出现频率相同。
如果仅统计一条链,则虽然A和T、C和G的出现频率不同,但是非常接近。
例如,对M.jannaschii基因组DNA的单向链进行统计分析,其结果见表5.2。
由于基因和其它功能区域在正反两条链上出现的可能性通常一样,所以正反两条链在信息的组织结构方面不应该有差别,因而核苷酸出现频率也不应该有偏差。
再根据正反两条链碱基互补的原则,可以推知,单链上A和T、C和G的出现频率相近。
表5.2M.jannaschii单链核苷酸出现频率
0.344
0.155
0.157
0.343
进一步的统计分析可以发现,不同基因组中两个连续核苷酸出现的频率也是不相同的。
4种核苷酸可以组合成16种两联核苷酸。
为了计算每个核苷酸对(i,j)的出现频率,扫描整个基因组DNA序列,累计核苷酸j在核苷酸i之后出现的次数,再除以总核苷酸对的数目。
在单条DNA序列上,总核苷酸对的数目等于序列的长度减1。
对酵母基因组两联核苷酸的统计结果见表5.3,其中核苷酸对出现频率最高的达到0.119,而出现频率最低的只有0.028。
表5.3酵母基因组两联核苷酸频率表
核苷酸对
AA
0.1193400681800
AC
0.0520605330203
AG
.0558********
AT
.0975********
CA
0.0583060967492
CC
0.0325646199051
CG
.028*********
CT
GA
.0557********
GC
0.0348050746970
GG
GT
TA
.0915********
TC
TG
TT
如果令Pij代表两联核苷酸(i,j)的出现频率,Pi代表核苷酸i的出现频率,则Pij’=Pij/(PiPj)的值反应核苷酸i和j的先后关系。
如果Pij’=1,则在两个连续的位置上,核苷酸i和j的出现是相对独立的。
对于酵母基因组,PA=0.3248,PAA=0.1193,PAA’=0.1193/(0.3248*0.3248)=1.131>
1,表明在两个连续位置上“A”的出现不是独立的,而是相关的。
同样,对于相隔一定距离k(k代表核苷酸个数)的两个核苷酸,也可能具有一定的相关性。
假设Pij(k)代表核苷酸j出现在核苷酸i之后第k个位置的频率,则可定义一个反应统计相关性的互信息I(k)(HerzelandGroß
e,1997):
(5-5)
I(k)值得大小实际上反应了距离为k的两个核苷酸之间的相关性的程度。
在进行编码区域识别时,常常需要对三联核苷酸进行统计分析,这实际上是分析密码子的使用偏性。
由于密码子的简并性(degeneracy),每个氨基酸至少对应1种密码子,最多有6种对应的密码子。
在基因中,同义密码子的使用并不是完全一致的(Karlinand1996;
Ghosh2000)。
不同物种、不同生物体的基因密码子使用存在着很大的差异(Granthametal.,1980,1981;
GouyandGautier,1982)。
进一步的研究表明,不同的密码子使用模式的形成可能与基因的GC含量有关。
在一些单细胞生物如Escherichiacoli、Saccharomycescerevisiae中,高表达的基因密码子的使用偏性一般比较大,这主要是由于基因的碱基组成和mRNA翻译时的tRNA选择两大因素造成的(Ikemura1981;
Sharpetal.,1986)。
最近的一些研究表明,基因密码子的使用与基因编码的蛋白的结构和功能有关,基因密码子的使用与基因表达的生理功能有着密切的联系(Heleneetal.,1998;
Richardetal.,2000)。
mRNA中的稀有密码子的使用与蛋白质结构域的连接区和规则二级结构单元的连接区有关,翻译速率在连接区会降低(Thanarajand1996);
同时,在表达具有不同二级结构的蛋白质时,mRNA区段的翻译速率有所不同(ThanarajandArgos1996)。
这些均说明蛋白质折叠方式与mRNA序列之间存在一定的相关性。
最近的研究表明,蛋白的三级结构与密码子使用概率有密切的关系,通过对密码子的聚类分析,可以很清晰地将具有不同三级结构蛋白质的编码基因分成不同的类,而具有相似三级结构蛋白的编码基因则大致聚在同一类中,从而证明基因密码子的使用偏性与蛋白质三级结构具有密切的相关性。
不同基因的密码子的使用与所编码的蛋白的三级结构的相关性的研究,可以用来对未知基因的密码子使用偏性情况进行预测。
同时,这一相关性还可用来对未知空间结构的蛋白质进行空间结构的预测,从而预测蛋白的功能,这在后基因组时代有着非常重要的研究意义(顾万君等,2001)。
进一步的研究发现,在不同物种中,类型相同的基因具有相近的同义密码子使用偏性,对于同一类型的基因由物种引起的同义密码子使用偏性的差异较小(周童等,2001)。
第三节功能位点分析
所谓功能位点(functionalsite)就是与特定功能相关的位点,是生物分子序列上的一个功能单元,或者是生物分子序列上一个较短的片段。
功能位点又称为功能序列(functionalsequence)、序列模式(motif)、信号(signal)等。
核酸序列中的功能位点包括转录因子结合位点、转录剪切位点、翻译起始位点等。
在蛋白质序列分析中,常使用序列模式这个名词,蛋白质的序列模式往往与蛋白质结构域或者作用部位有关。
基因组序列中若干个相邻的功能位点组合形成功能区域(functionalregion)。
不同的功能位点可能相互覆盖或相交。
密码子是蛋白质编码序列中的功能单元,密码子与tRNA相互作用。
与单个氨基酸相对应的密码子是非重叠的。
其他的功能单元(如蛋白质结合位点)比较复杂。
对于未知功能的序列进行功能信号检测是DNA序列分析的一个重要任务。
1、利用共有序列(consensus)搜索功能位点
共有序列又称一致性片段,它是描述核酸序列中功能位点的最常用方法(SolovyevandKolchanov,1997)。
共有序列是关于功能位点特征的描述,它描述了功能位点每个位置上核苷酸进化的保守性,而这种保守性是与功能相关的。
目前国际上有许多数据库包含各种DNA或RNA功能位点的共有序列,如基因转录调控元件。
利用共有序列进行功能位点分析牵涉到两个方面的问题,一是如何构造共有序列,二是如何利用共有序列在给定的核酸序列上搜索寻找功能位点,并计算所找到的功能位点的可靠性。
共有序列具有以下几个方面的特征:
(1)共有序列中既有保守的位置,在这些位置上仅允许出现特定类型的核苷酸,也有可变的位置;
(2)任何位置上的核苷酸可以用表5.4中15种类型之一来表示:
表5.4核苷酸表示符号
核苷酸类型编号
符号
内容
1
2
3
4
5
R
A,G
6
Y
T,C
7
W
A,T
8
M
A,C
9
K
G,T
10
S
G,C
11
D
A,G,T
12
V
A,G,C
13
B
C,T,C
14
H
T,C,A
15
N
A,G,T,C
(3)在分析一个共有序列与一个DNA片段的关系时,必须评估其统计的显著性;
(4)由于共有序列一般很小,所以不用各种核苷酸出现的频率来描述共有序列,而只用各种核苷酸出现的次数来描述共有序列。
设一个共有序列的长度为L,用Nl来刻画该共有序列的核苷酸组成特征,Nl是属于对应功能位点第l类核苷酸的个数。
显然有:
(5-6)
设该共有序列有M个保守位置,用Ml来描述这些保守位置,Ml是保守位点第l类核苷酸的个数,显然有:
(5-7)
注意,在公式(5-6)中,加和过程不包括“15”类,因为如果某个位置上核苷酸的类型号为“15”(类型为N,代表任意核苷酸),则该位置是非保守的。
共有序列与一个长度为L的序列片段之间同源变化仅可能出现在L-M个位置上,被第l类核苷酸所占据的可变位置数等于Nl-Ml。
可利用共有序列分析DNA序列中的功能位点或搜索序列数据库。
当然在这之前首先要构造共有序列,其前提条件是获得多个实际同类功能位点的序列多重比对。
具体的构造过程如下(Rogozinetal.,1997):
(1)初始化共有序列为一系列可变位置,以“N”代表;
(2)在可变位置寻找出现次数最多的核苷酸,并将该位置转化为保守位置;
(3)对当前所得到的共有序列进行特异性检查,若通过检查,转(5),否则转(4);
(4)形成与当前共有序列一致的位点子集,转
(2);
(5)从原位点集合中删除与当前共有序列一致的位点,若还有剩余位点,则转
(1),构造另外的共有序列。
[1][2][3][4][2][3][4][2]
TTATGTTATGTNNNNtTATGtTATGTNNNCtACGctACGc
ATATAATATA非特异tACGCTacgc非特异tTGTctTGTc
TACGCTACGCtTGTCtTGTCtCCActCCAc
TTGTCTTGTCtCCACtCCACTNSNC
TCCACTCCACTNNNC
NNNNNTNNNN
[3][5][5]
TNSNCConsensus1:
Consensus2:
特异TNSNCNTATN
剩余位点:
TTATG
ATATA
图5.1共有序列构造示例
图5.1是一个构造共有序列过程的示例。
在给定的序列中搜索与共有序列一致(或者匹配)的序列片段是一个很简单的过程,只要检查序列片段在对应位置的成分是否与共有序列的描述相符合。
需要注意的是,在很长一段DNA序列中或者在众多的数据库序列中找到一条符合条件的序列片段可能并不一定是我们期望的结果,可能是由于随机因素的影响使一个片段被选中。
共有序列越短、被搜索的序列越多、越长,则随机影响越大。
必须对查找结果进行统计的显著性分析,只有当显著性大于某个给定的阈值时,才能认为是真正找到了与共有序列相一致的序列片段。
共有序列是关于序列特征的一种定性描述。
对于DNA序列,它能够说明序列每个位置可能出现的碱基类型,但是不能准确地说明各位置上不同类型碱基出现的可能性大小。
因此需要定量的序列特征描述方式。
2、用感知矩阵分析功能位点
这种方法用权系数描述功能位点各位置上每种核苷酸的相对重要性。
感知矩阵(或加权矩阵)是根据一系列功能位点的多重对比排列结果而建立的,其大小为4n,4代表碱基的种类数目,n代表功能位点的长度。
图5.2即是一个感知矩阵M,其中每行代表不同的核苷酸,而每列对应于功能位点上的特定位置。
矩阵的每一个元素M(a,j)的值代表第a种核苷酸在功能位点第j个位置上出现的得分,a{A,T,G,C}。
22
-3
19
-1
-5
-19
-9
16
图5.2感知矩阵示例
令感知矩阵描述的位点长度为n,对于给定一个序列s=a1a2…an,根据对应位置上核苷酸的类型,取感知矩阵中对应的权值,加和以后得到该序列的得分
(5-8)
设S=ATTGCA,则
Ws=1+6+14-5+8+19=43
假设T和T'
分别是功能位点阈值和非功能位点阈值。
在检测功能位点信号时,如果WsT,则S是功能位点;
如果WsT'
,则S是非功能位点。
设A代表功能位点的训练集合,将该集合分解为两个部分,其中A+代表功能位点集合,A-代表非功能位点集合。
感知矩阵M的构造算法如下:
(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)若训练集合中的所有实例都处理过,则循环结束,转(7),否则继续执行循环体,直到处理完所有实例;
(7)如果M稳定,则结束;
否则转
(2)。
上述算法反复调整感知矩阵M的元素值,直到M矩阵能够正确识别训练集中的所有功能位点和非功能位点。
对于最终得到的感知矩阵,要求其具有敏感性和特异性,每一列上的元素值应该尽可能地有明显的差别,以便反应功能位点各个位置上的特点。
与感知矩阵类似,如果令矩阵每一个元素M(a,j)的值代表第a种核苷酸在功能位点第j个位置上出现的概率,则M是一个概率矩阵,与第三章中的特征统计图相对应。
假设各个位置上出现的碱基是相互独立的,即任何两个位置上的碱基是不相关的,那么对于给定一个序列s=a1a2…an,可以计算出功能位点序列为s的概率:
(5-9)
如果分别统计功能位点和非功能位点,通过计算可以形成两个矩阵M和M’,进一步计算可以判断一个给定的序列究竟属于功能位点,还是属于非功能位点。
给定一个序列s=a1a2…an,定义似然比LR(M,M’,s):
(5-10)
在进行功能位点检测时,计算LR(M,M’,s),并与给定的阈值L比较,如果LR(M,M’,s)>
L,则序列s可能是一个功能位点。
为了能够在M(a,j)或M’(a,j)非常小的情况下正确计算,将似然比的定义改变为对数似然比,其计算公式如下:
(5-11)
概率矩阵M和M’的每个元素是一个0和1之间的正数。
如果令一个4n新矩阵U的元素(a,j)的值为log2(M(a,j)/M’(a,j)),则矩阵U的每个元素值可能是正值,也可能是负值。
实际上,矩阵U就是前面介绍的感知矩阵。
第四节隐马尔柯夫模型
1、马尔柯夫链(Markovchain)(Feller,1968;
张新生等,2000)
考虑一个具有多个状态的系统S,S={s1,s2,…,s|s|},令S0、S1、…、St为一系列在各个时刻系统状态的变量,即状态链。
对于每个1到|S|的整数,它们分别与状态链中的一个状态相联系,并且在任何时刻,这条链都处于一个特殊的状态。
当且仅当对于任何t有
(5-12)
则St形成一条马尔柯夫链。
简单地说,就是系统未来的状态仅依赖于当前状态。
St称为在时刻t系统链的状态。
一条马尔柯夫链完全决定于初始分布P(S0)和转换概率Pt=P(St+1|St)。
我们仅讨论确定模型的马尔柯夫链,即状态转换概率不随时间变化的马尔柯夫链。
令状态转换矩阵为F=(fij),fij代表从状态si移动到状态sj的概率。
生物序列可以被描述为一个随机过程的输出,其中对于一个给定的核酸在位置p出现的概率依赖于已占据前面k个位置的核酸,这样一种表示称为k阶马尔柯夫模型。
一个序列具有不同的统计性质(如二目频率或三目周期性),不同的功能区域(如编码区域、非编码区域)对应于不同的马尔柯夫模型。
作为一个例子,下面介绍马尔柯夫链在识别CpG岛中的应用。
CpG岛是一类长度在几百bp的特殊DNA序列,其中CG核苷酸对出现的频率非常高。
CpG岛在基因组中有重要的生物学意义,而识别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,
(5-13)
整个序列X的发生概率为
(5-14)
为了处理方便,添加两个特殊的字符‘B’(begin)和‘E’(end),使得x0=‘B’,xL+1=‘E’,则上述公式简化为:
(5-15)
令fst+为CpG岛内的字符转换概率,fst-为CpG岛外的字符转换概率(见表5.5),则X的对数似然得分为
(5-16)
上述计算值越大,则X越可能是CpG岛。
表5.5CpG岛内部和外部的转换概率(Shamir,1999)
+
ACGT
-
0.1800.2740.4260.120
0.300
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- DNA 序列 分析 东南大学 生物 信息学 实验室