真核生物5‘端外显子和翻译起始位点的预测

北京大学校长基金论文集（2003年）真核生物翻译起始位点预测1真核生物5’端外显子和翻译起始位点的预测PredictionofFirstExonsandTranslationInitiationSitesinEukaryoticGenomes力学与工程科学系00级胡钢清GangqingHuDept.ofMechanicsandEngineeringScience摘要真核生物翻译起始位点（TIS，translationinitiationsites）的正确预测对于基因的正确注释有着重大的意义。本文试图通过发展和结合一些现有的算法来较好地预测出真核生物翻译起始位点。本文算法的得出主要基于三种方法自然结合。我们首先基于简单的权重矩阵方法（WMM，WeightMatrixMethod）【1】，发展了多状态信号的熵距离判别分析方法,由此可以得到一些待定ATG信号，这些信号包含了大部分的翻译起始位点。接下来，我们采用WMM找出与该ATG对应的GT信号并将它们之间的序列提出（这里ATG信号以及GT信号的定义详见正文）。通过引入模糊词汇的概念，我们发现翻译起始位点前存在某些模糊词汇，这样可以通过判断ATG信号前是否含有某个“模糊词汇”进一步选出待定ATG信号。最终我们结合了佘振苏教授提出的以统计语言学为特色的多元熵距离（multivariateentropydistance，MED）方法【2】对这些ATG信号对应的阅读框在20维空间中进行多中心聚类形成了预测翻译起始位点的完整算法并达到了较好的结果。1引言2001年2月，人类基因组计划的相关组织和机构公布了人类基因组测序草图【3】，人类基因组计划的提前完成，基因组的研究热点已转向揭示基因信息结构的复杂性与遗传语言的根本规律。其中，基因预测算法的研究也成为对基因组序列进行统计分析的重要目标。所谓基因预测，就是在对DNA序列编码潜能（codingpotentials）提出某种模式（pattern）描述的基础上，对一未知的DNA序列上完整的基因结构进行注释。对于原核生物，由于其基因结构较为简单，在基因组的DNA链上表现为一个编码蛋白质的基因对应为一段连续的开放阅读框（openreadingframe，ORF），因此，基因预测的问题也相对简单，本人所在的佘振苏教授研究小组已经在这方面取得了很好的结果【4】。但是对于真核生物，它们的基因结构远比原核生物的复杂，许多基因是断裂基因，间断成外显子（exon）和内含子（intron），并且exon在序列中长度比例极小【5】。生物体通过对剪接位点（SpliceSite）的剪接等过程最终构成蛋白质。现有很多真核生物基因结构的预测算法就是针对SpliceSite的预测，即对外显子的预测，著名的有Genescan，GeneMark等。本质上可以将外显子分为4类：包含起始编码子（startcodon）和终止编码子（terminationcodon）的singleexon，包含起始编码子和供体位点（donorsite）北京大学校长基金论文集（2003年）真核生物翻译起始位点预测2的firstexon，包含受体位点（acceptorcite）和供体位点的internalexon，以及包含受体位点和终止编码子的terminalexon【6】。现有基因预测软件中对于firstexon、terminalexon以及singleexon，它们的识别精度通常都大大低于internalexon的识别精度【7】。Firstexon的预测实际上包含了对翻译起始位点的预测。本文重点在于翻译起始位点的预测算法设计，同时也给出了5’端外显子（即firstexon）的预测结果。一般认为，对基因结构预测所用到的信息有三部分组成：序列上的信号（signal），内容（content）统计以及与已知基因的相似度（similarity）【8】。我们用到的信息只是包含了前两部分，对ATG信号的多状态熵距离判别分析以及找出模糊词汇是属于信号的范畴，而对阅读框运用MED方法实际上关系到了内容统计。在预测时，我们采用的算法是简单的线性判别。2材料我们采用了Kulp&Reese（1995）、ALLSEQ【9】以及HMR195【7】三部分数据集，其中前两者的部分序列作为学习集，最后者的部分数据作为测试集。如引言所提到，真核生物的基因按是否含有内含子可以分为多外显子基因以及单外显子基因。实际上，前者在数目上占绝大多数，因而本文关心的是对多外显子基因中的翻译起始位点以及5’端外显子的研究与预测。Kulp&Reese（1995）提供的数据集是通过运用BLASTP【10】删除一些相似和冗余的序列后得到的。数据库包含了人类的单外显子基因186个，多外显子基因304个，其中多外显子基因中共包含1798个外显子。我们直接采用了其中的304个多外显子基因作为学习集的一部分。ALLSEQ数据集经过了一些严格的标准【11】筛选、整理出来的。它包含的是一些脊椎动物的多外显子基因，有570个，其中有外显子数目为2649。我们采用ALLSEQ数据集构成学习集的另一部分。HMR195数据集是Rogic等于2001年整理出的，其中包含了人、小鼠和大鼠共195条基因的注释。数据集中人、小鼠、大鼠的比例为103：82：10，包含152个多外显子基因和43个单外显子基因。我们采用其中152个多外显子基因作为测试集。这样我们一共得到了874个基因作为学习集，152个基因作为测试集。3方法3．1ATG信号的多状态熵距离判别分析方法真核序列中，基因一般是由编码甲硫氨酸（Met）的密码子ATG开始，我们把这样的ATG称为翻译起始位点，而将DNA序列上出现的核苷酸三联码“ATG”称为ATG信号。为了描述ATG信号，我们取出由ATG信号前mbp，后nbp组成的核苷酸序列L。用qij（j＝1,…,m+n，i=0,…,3分别表示核苷酸A、C、G、T）表示出现在该序列上第j个位置上的核苷酸i，qij是一个4×(m+n)的矩阵，其特征如图1所示。北京大学校长基金论文集（2003年）真核生物翻译起始位点预测3图1用矩阵q来反映序列的核苷酸的位置以及内容我们可以简单地将ATG信号分为两类，一类是翻译起始位点，剩下的ATG信号归为另一类。为了刻画这两类信号的区别，我们分别统计出他们对应于L的权重矩阵0M、1M（它们都是4×（m＋n）的矩阵）,进而依次定义出该序列的相对熵kLH以及熵距离LE如下：310log(k12mnijkLijkjiijqHqM其中＝，），01LEH-HLL然而通过引入多状态熵距离判别分析，我们将对ATG信号有更细致的描写。为了描叙处于不同状态的ATG信号（“状态”这个词将在下文阐明），我们要用到阅读框这个概念：在同一DNA序列上，如果用TISp、ATGp来分别表示翻译起始位点以及某个ATG信号的位置，则可以定义出该ATG信号的阅读框ATGATGmod(,3)TISfrpp（表示ATGTISpp对3的求模运算），可以看出翻译起始位点对应ATG信号的阅读框为0。阅读框概念的引入是有生物学意义的，生物体通过对剪接位点剪接等操作形成mRNA，在mRNA的核苷酸序列中只有与翻译起始位点有相同阅读框的3个核苷酸才可能编码氨基酸。在第一个供体位点前，ATG信号在DNA序列上出现的位置决定了该ATG可能具有的某个属性：N：出现在翻译起始位点之前，或出现在翻译起始位点之后且阅读框与翻译起始位点对应的阅读框不同。S：出现在翻译起始位点位点上。C：出现在翻译起始位点之后且阅读框与翻译起始位点对应的阅读框相同。我们猜想，每一个ATG信号并不是独立的，它与邻近的其它ATG信号存在某些关联，尤其对于具有同一阅读框的ATG信号这种关联将更为强烈。这样的话，考察同一阅读框的3个相邻ATG信号之间的关联性，可以归纳出中间那个ATG信号可能具有的6种状态：1，当ATG信号与翻译起始位点有相同阅读框时，有5种：NNN_0（NNN_0,表示这3个ATG的属性依次为：N，N，N。后面的数字表示这3个ATG的阅读框，后类推），NNS_0，NSC_0，SCC_0，CCC_0。2，当ATG信号不与TIS同阅读框时归为第6种：NNN_1,2。为了刻画具有NSC_0状态的ATG信号与具有其它状态信号之间的区别，我们取出与其有相同阅读框的前后各一个ATG信号，将这三个ATG信号各自对应的取其前mt后nt核苷酸序列Lt（t=1,2,3）拼接起来构成一条序列总长为l＝北京大学校长基金论文集（2003年）真核生物翻译起始位点预测431()tttmn的序列S（如图2）。图2将三条序列拼接起来，构成一条新的序列S（这里假定（m，n）分别为（4，7）、（2，9）以及（5，6））在统计出不同状态对应于S权重矩阵1,...,6iMi)（＝后，类似于将ATG信号分为两类时相对熵的定义，我们给出将ATG信号分为6个状态时的相对熵的定义310log(lijksijkjiijqHqM其中k=1,2,...,6)为了刻画具有状态NSC_0的ATG信号与其它状态的ATG信号的差别，我们同样可以引入熵距离的概念，不过这里的熵距离是一个5维向量：13456(,,,,)ssssssEEEEEE2(kkssEH-Hs其中k=1,3,4,5,6）如果认为非编码区内的核苷酸之间是相互独立的以及结合生物学上每个氨基酸是由3个核苷酸编码而成这个事实，我们可以认为将序列Lt划分为pt个片段Ltj（其中311tjpttiLS，如图3），并将每个片段统计出的熵距离组成一个向量1(,...,)pEEE是有意义的。这样S就被化为n＝31ttp个片断，第m个片断的相对熵、熵距离以及S的熵距离依次为：213,120log(,SlijksmijkjliijqHqllM这里表示第m个片段在上的起始与结束位置）134562,,,,,,,,(,,,,),,kksmsmsmsmsmsmsmsmEEEEEEEH-Hs,m其中k=1,3,4,5,6113366,1,,1,,1,(,...,,,...,,...,,...,)sssnssnssnEEEEEEE北京大学校长基金论文集（2003年）真核生物翻译起始位点预测5图3将序列Lt划分为各个片断可以看出，每一ATG信号都可以在5×n维的熵距离空间中对应一个向量sE。这样对于给定的一个ATG信号，我们构造出其sE，计算该sE与各种状态ATG信号对应sE中心的距离lk，取dk=σk*lk(k=1,2,…6,0σk为参数)，取最小值对应的状态为该ATG最可能的状态，进而判断出该ATG的属性。表1给出了在相同条件下，采用多状态熵距离判别分析方法（各种参数如表2）与WMM【1】对测试集的ATG信号进行分析的对照结果。从表中可以发现：（1）对于在第一个内含子前面的序列，在要求较高精度（93%）地识别出TIS的条件下，采用多状态熵距离判别分析方法总比WMM方法能减少不少假阳性。（2）多状态熵距离判别分析方法从WMM方法中减少假阳性ATG信号数目的能力(减数/WMM数目)是随分析范围的从5’端往后延伸而降低的。即与WMM相比，保守地估计多状态熵距离判别分析方法只在分析第一个受体位点前的ATG信号占优势。而这对我们设计可以算法就已经够用了，如下文可以看出我们只需要学习第一个受体位点前信息。表1多状态熵距离判别分析方法与WMM的对照结果（这里把判断为TIS的不具有NSC_0状态的ATG信号成为假阳性的ATG信号）分析到第一个供体位点之前分析到第一个受体位点之前TIS=ATG信号假阳性的ATG信号减少数目减数/WMM数目TIS=ATG信号假阳性的ATG信号减少数目减数/WMM数目多状态熵距离数目WMM数目多状态熵距离数目WMM数目134832848160.01913414891403-86-0.0613810131040270.02613717681764-4-0.00143123513651300.09514221232157320.0151451300