在序列同源性检测中完成模式词典组成的方法和设备

摘要

在本发明的词典组成方面,一种基于计算机的用于处理数据库中多个序列的方法包括以下步骤。首先,该方法包括对多个包括用于组成每个序列的字符在内的序列中的每一个进行评价。然后,生成至少一个表示数据库中序列的至少一个子集的字符模式。该模式具有一个与其相关连的统计学重要性,该模式的统计学重要性由一个用于表示在数据库中模式所支持的序列的最小数量的数值所确定。

说明书

本申请对以A.Floratos和I.Rigoutsos的名义于1998年10月30日提交的临时美国专利申请第60/106,295号要求优先权，并涉及以A.Floratos和I.Rigoutsos的名义与其同时提交的名为“完成序列同源性检测的方法和设备”的美国专利申请，在此将它们引作参考。

本发明一般涉及数据库搜索，及更具体地涉及用于检测查询序列和数据库中与一个给定的应用例如遗传研究相关连的序列之间的序列同源性的方法和设备。

在遗传研究领域内，在将一个新基因测序之后的第一步是鉴定该基因的功能的努力。最普遍和直观的用于完成此目标的方法采用以下生物学事实-如果两个肽延伸在序列水平上表现出足够的相似性(即只须通过少量的插入、缺少和/或氨基酸突变即可从一个得出另一个)，则它们可能在生物学上是相关的。这类方案的例子描述于A.M.Lesk的“计算分子生物学”(计算机科学和技术百科全书，A.Kent和J.G.WilliamS编辑，Marcel Dekker,New York,1994)；R.F.Doolittle的“我们从序列数据库中学到了和将要学到什么”(Computers and DNA,G.Bell和T.Marr编辑，21-31，Addison-Wesley,1990)；C.Caskey,R.Eisenberg,E.Lander和J.Straus的“DNA专利的Hugo声明”(Genome Digest,2：6-9,1995)；以及W.R.Pearson的“蛋白质序列比较和蛋白质进化”(分子生物学中智能化系统的讲座，Cambridge,England,1995)。

在此框架之内，得到关于一个新基因的功能的线索的问题成为氨基酸串中鉴定同源性的问题。一般而言，同源性系指两个或多个序列或串之间的相似性或关系。因此，给定一个查询序列Q(例如新基因)和一组特性很清楚的蛋白质集合D后，需要寻找Q之中相似于D中序列区域的所有区域。

用于实现此任务的第一方案基于一种称为动态规划的技术。此方案描述于S.B.Needleman和C.D.Wunsch的“可用于搜索两个蛋白质的氨基酸序列中相似性的一般方法”(Journal Of MolecularBiology,48：443-453,1970)；以及T.F.Smith和M.S.Waterman的“公共分子子序列鉴定”(Journal Of Molecular Biology,147：195-197,1981)中。不幸的是，此方法的计算量要求很快地使它成为不实际，特别是当如今那样搜索大型数据库时更为如此。一般而言，该问题在于动态规划变体花费很多时间来计算同源性，而这事实上并不重要。

为解决此问题，已经建议了不少算法，它们只集中于发现很局部的相似性。这些算法中最有名的当推FASTA和BLAST。FASTA算法描述于W.R.Pearson和D.J.Lipman的“生物序列比较的改进工具”(Proc.Natl.Acad.Sci85：2444-2448,1988)；以及D.J.Lipman和W.R.Pearson的“快速的和灵敏的蛋白质相似性搜索”(Science,227：1435-1441,1989)中。BLAST算法描述于S.Altschul,W.Gish,W.Miller,E.W.Myers和D.Lipman的“一个基本局部对比搜索工具”(J.Mol.Biology,215：403-410,1990)中。在大多数情况下，可以首先寻找无空位(ungapped)同源性，即只由于突变而不是由于插入或缺少而引起的相似性，从而提高性能。此方案的逻辑依据在于，在两个肽串之间的任何无空位同源性中，有机会至少存在着一对子串，它们的匹配不包含空位。对这些子串的定位(无空位同源性)可用作获取整个(有空位)同源性的第一步。

然而，鉴定查询和数据库序列之间的相似区域只是过程的第一部分(计算量最大的部分)。第二部分(生物学家最感兴趣的部分)是评价这些相似性，即判定它们是否足够以便维持查询和相应的数据库序列之间的推理关系(功能的、结构的或其他)。这类评价通常通过组合生物学信息和统计学推理而完成。通常，相似性量化为一个为每对相关区域计算的分数。此分数的计算涉及空位罚分(对于有空位对比而言)和那些具有任何给定的氨基酸变化为另一种氨基酸的进化概率的合适突变矩阵的使用。这些矩阵的例子是PAM矩阵(见M.O.Dayhoff,R.M.Schwartz和B.C.Orcutt的“蛋白质中进化变化模型”(Atlas of Protein Sequence and Structure,5：345-352,1978))和BLOSUM矩阵(见S.Henikoff和J.G.Henikoff的“来自蛋白质模块中的氨基酸替代矩阵”(Proc.Natl.Acad.Sci89：915-919,1992))。然后通过计算这一分数能够单纯地靠偶然性上升的概率(在某些统计学模型中)来评价此罚分的统计学重要性，例如见S.Karlin,A.Dembo和T.Kawabata的“分子序列中高评分片段的统计学组成”(The Annals of Statistics,2：571-581,1990)；以及S.Karlin和S.AltSchul的“使用一般评分方案评价分子序列的统计学重要性的方法”(Proc.Natl.Acad.Sci87：2264-2268,1990)。取决于使用的统计学模型，此概率可以决定于一些因素，例如：查询序列的长度，基础数据库的大小等。然而，无论使用什么常规统计学模型，总是存在着所谓“灰区”，例如统计学上不重要的分数标示事实上在生物学上的重要相似性。因此不幸的是，对于一个统计学模型能够如何好地接近生物学实际，最后不可避免总有一个限制。

解决统计学重要性与弱相似性的关系的内在困难的选代方法是在推断用作进化远缘同源性的模型的序列描述中使用生物学知识。BLOCKS(见S.Henikoff和J.Henikoff的“数据库搜索中蛋白质模块的自动组装”(Nucleic Acids Research,19：6565-6572,1991))是一个采用模式推理概况的系统，这些概况是通过在PROSITE(见S.Henikoff和J.Henikoff“基于块数据库搜索蛋白质家族分类”(Genomics,Vol.19,pp.97-107,1994))数据库中定义的蛋白质分类而获取的，用于功能地评注新基因。其优点在于此分类是由熟悉已知相关的蛋白质家族的专家编辑的。其结果是，即使弱相似性也能鉴定和使用于评注过程中。另一方面，需要很多关于哪些蛋白质的确是相关的并能随后由模式表示的知识。此外，总有以下危险，即一个蛋白质家族实际上包含比目前认为的更多的成员。将这些其他成员排除于考虑之外时，有可能获得对该家族“过度匹配”的模式，即过于严格而不能将它们外推至未鉴定的家族成员中。

因此，显然需要一种方法和设备，用于通过独一的词典组成技术建立改进的模式词典，该词典组成技术提供改进的序列同源性检测，以及需要其本身不限于只搜索评注的序列的用于序列同源性检测的方法和设备。

本发明提供下面将更详细地描述的改进的模式词典组成技术和改进的序列同源性检测技术而为以上所述和其他需要提供解决方法。

在本发明的序列同源性检测方面，一个基于计算机的在数据库中多个序列和查询序列之间检测同源性的方法包括以下步骤。首先，该方法包括访问与数据库相关连的模式，每个模式表示数据库中的一个或多个序列的至少一部分，其次，查询序列与模式比较以便检测查询序列的一个或多个部分是否与由模式所表示的数据库的序列的一个或多个部分同源。然后为每个检测出与查询序列同源的序列生成一个分数，其中序列分数是基于根据所检测序列的每个同源部分生成的个别分数而得，该序列分数表示查询序列与检测的序列之间的同源性程度。

在本发明的词典组成方面，一个基于计算机的用于处理数据库中多个序列的方法包括以下步骤。首先，该方法包括评价多个包括用于组成每个序列的字符在内的序列中的每一个序列。然后，生成字符的至少一个模式以表示数据库中序列的至少一个子集。该模式具有与其相关连的统计学重要性，模式的统计学重要性由用于表示模式在数据库中支持的序列的最小数决定。

因此，与现有技术显著不同，本发明的方法学的依据是在任意的数据库上完成无指导的模式发现而不要求数据库的任何现有划分。BLOCKS方案假定数据库已被(外部专家)划分为生物学上相关序列的子集。个别地处理每个子集即获得概况。此方案的结果是，BLOCKS无法处理任意数据库，由于并非所有这类数据库都划分为相关的子集。事实上，BLOCKS只用于此处引用的SwissProt数据库，它使用描述于也在此处引用的PROSITE数据库中的蛋白质组。另一方面，本发明优选地使用整个数据库作为其输入及提供一个自动化的方法学，用于决定哪些模式是重要的和哪些不是。

此外，本发明提供一个新的统计学框架，用于评价所发现的模式的统计学重要性。不像现有框架，本发明的方案在其计算中引入存储器的概念。也即，例如当查询序列上的区域A与某数据库序列上的区域B比较时，在评价最终相似性分数时考虑到A与数据库中所有其他序列的相似性。

使用此处描述的加强统计学模型可以检测重要局部相似性，而它们在现有方案中将无法检测到。这允许本发明的系统在完成相似性搜索中能够比使用现有技术的系统具有更高灵敏度。

还有，本发明提供一个自动化方法，可以利用基础数据库D中可用的部分评注信息。此方法学允许用户利用看上去不重要的更细的相似性。例如，当一个模式与查询序列区域A匹配时，可以观察到也与该模式匹配的所有数据库区域。如果所有(或更多)这些数据库以相同方式评注，则此评注可以转送至查询区域A。可以证明，用以上方式部分地评注查询序列，对于总的序列评注是有用的。

本发明还提供一个详细的方法学，用于将数据库聚类为高度同源序列的组。在遗传数据处理应用中，此方法学允许对多结构域蛋白质进行正确处理。

还知道，此处描述的本发明概念可以在一个网络例如因特网上以客户-服务器关系实施。这允许用户在一个远程客户设备上输入一个查询序列，后者通过网络传送至服务器并在服务器中处理。该服务器然后通过网络向用户的客户设备返还同源性搜索的结果。

下面结合附图的本发明实施例的详细描述将使本发明的这些和其他目的、特征和优点更为明显。

图1是根据本发明实施例的序列同源性检测系统的框图；

图2是本发明的序列同源性检测系统的示例性硬件实施框图；

图3是本发明的序列同源性检测系统的基于网络的实施框图；

图4是用于阐述根据本发明一个实施例的搜索引擎方法学的高级流程图；

图5阐述根据本发明一个实施例的给定查询序列的模式匹配过程例子；

图6阐述根据本发明一个实施例的为一个具体查询序列生成的散列表例子；

图7阐述根据本发明一个实施例的给定查询序列的链接过程的例子；

图8阐述按照根据本发明一个实施例的给定查询序列的给分过程所生成的有向加权图例子；

图9是用于阐述本发明的搜索引擎方法学的匹配和链接阶段的实施例的流程图；

图10是用于阐述本发明的搜索引擎方法学的评分阶段的实施例的流程图；

图11阐述在SP34中具有给定主干结构的模式分布及与具有相同主干的随机分布的比较；

图12至15是用于阐述根据本发明一个实施例的词典组成方法学的流程图；

图16是用于阐述根据本发明一个实施例的数据库清除过程的流程图；及

图17至30阐述与本发明相关连的实验结果。

下面将在示例性遗传数据处理应用的范围内解释本发明。然而应该理解，本发明不限于这类具体应用。相反，本发明更一般地应用于从任意数据库中建立模式词典(在恰当地翻译数据库记录为等效的序列表示之后)以及相对于数据库中的数据为任何给定查询记录完成无限制同源性搜索。

最初参照图1，显示一个根据本发明一个实施例的序列同源性检测系统的框图。示例性系统100包括一个搜索引擎模块110，一个模式词典120，一个词典组成模块130和一个源数据库140。如下面将要解释的，搜索引擎110从一个用户接收一个查询序列并完成对模式词典120的搜索以便试图找出存在词典中的模式，该模式用于表示以某种方式相似于查询序列的数据库140中序列。在查询之前，词典组成模块130从数据库140中组成模式词典120。此词典组成过程系指信息收集或挖掘(mining)。搜索引擎110将某些或全部查询结果(例如来自数据库的同源性序列)返还给用户。

图2是序列同源性检测系统100的示例性硬件实施的框图。如图所示，系统100可以根据一个处理器210，一个存储器220和I/O设备230实施。应该知道，此处所用“处理器”一词试图包括任何处理设备，例如包括一个CPU(中央处理单元)的设备。此处所用“存储器”一词试图包括与一个处理器或CPU相关连的存储器，例如RAM，ROM，固定存储设备(例如硬盘驱动器)，可装卸存储设备(例如软盘)，闪烁存储器等。此外，此处所用“输入/输出设备”或“I/O设备”一词试图包括例如一个或多个用于进行查询和/或输入数据至处理单元的输入设备，例如键盘，以及/或一个或多个用于呈现查询结果和/或其他与处理单元相关连的结果的输出设备，例如CRT显示和/或打印机。还应该理解，“处理器”一词可以意味着不止一个处理设备，及与一个处理设备相关连的不同元素可以由其他处理设备共享。因此，如此处所描述的，软件部件包括用于完成本发明方法学的指令或代码，它们可能存于一个或多个相关连的存储设备中(例如ROM，固定或可装卸存储器中)，以及当准备好以供使用时，可以部分地或全部装载(例如装载入ROM中)并由一个CPU执行。

图3是本发明的序列同源性检测系统的基于网络实施的框图。如图所示，一个客户计算机系统310通过网络330例如因特网与一个服务器计算机系统320通信。然而，该网络也可是一个私有网络和/或局域网。根据图3的实施，图1中所示系统100的元素的全部或一部分位于服务器330上并由它执行。一个在他的客户计算机系统(例如个人计算机、膝上和/或某些其他类型个人处理设备)上远程地操作的用户通过在计算机系统中运行的应用软件(例如网络浏览软件和/或与搜索引擎相关连的图形用户界面)输入一个查询序列。该查询以常规方式通过网络330，并由服务器320处理。服务器320接收该查询并执行按照一个存储的模式词典的本发明搜索引擎方法学。该词典可能由本发明的词典组成模块根据一个源数据库组成。服务器通过网络返还部分或全部查询结果(例如来自数据库的同源性序列)至客户。应该理解，服务器可表示不止一个计算机系统。也即，图1中一个或多个元素可能位于它们自己的计算机系统(例如带有其自己的处理器、存储器和/或I/O设备)上并由其执行。

给出一个本发明序列同源性检测系统的元素的一般描述和不同示例性硬件实施后，现在详细地阐述本发明的不同方法学。

下面在一个与序列同源性检测系统100相关连的示例性实施例中描述与搜索引擎模块110和词典组成模块130相关连的相应方法学。应该理解，与搜索引擎模块相关连的本发明方法学也可以与其他已知模式词典一起使用。相似地，与词典组成模块相关连的本发明方法学可用于建立与其他已知搜索引擎一起使用的模式词典。

为便于参考，详细描述的其余部分将分为以下段落：(Ⅰ)定义；(Ⅱ)搜索引擎；(Ⅲ)词典组成；以及(Ⅳ)实验结果。

Ⅰ．定义

以下段落提供某些评注，以供下面描述本发明不同方面之用。

∑系指用于构成序列的字符集合。在生物学设置中(此处所优选的)，我们所关心的序列是蛋白质，该集合是20种氨基酸的集合。术语蛋白质/序列将要在下面交换地使用，对于术语字符/氨基酸的词也是如此。

D系指蛋白质的基础数据库，即在其上建立模式集合(模式词典或生物学词典)的数据库。我们将要在整个说明书中一直使用的一个示例性数据库如下(包含三个序列)：

D={s₁,s₂,s₃}，其中

s₁=ARQSTLUMNPQ

s₂=FDSALQFTGMRA

s₃=RKMFPQDDSLA

Ⅱ系指模式集合，即此处系指生物学词典或模式词典120。获取Ⅱ的正确方法将要在下面名为“词典组成”的段落中阐述。模式是用于描述肽家族的正规表达式。由一个单个模式表示的多肽家族预料会包含相关的(结构上、功能上、进化上)氨基酸的延伸。更具体地，给出氨基酸的字母∑后，我们定义Ⅱ中的一个模式P为以下形式的正规表达式：

∑(∑U{‘·’})^*∑

其中‘．’(系指“不必关心字符”)系指可由任意残基占据的位置。作为一个正规表达式，每个模式P定义一个包含可通过将一个∑中的任意残基替代每个不必关心字符而从P中获得的所有串的多肽语言。还有，Ⅱ中的每个P与D中的至少K_min个序列相匹配。K_min是一个整数，其计算在下面“词典组成”段落中阐述。在下面的例子中，我们将要假定一个具体值。与模式P匹配的数据库序列的区域记录于模式P的偏移量表L_D(P)中。这是一个表，包含所有对(j,k)以使模式P与在偏移量k处的数据库第j个序列相匹配。

对于以上引入的数据库例子，并假定K_min=2，则模式集合为P={A.Q.T,M.PQ}。此集合中的两个模式出现在下面输入序列(匹配位置用下划线表示)中：

    A.Q.T                  M.PQ

S₁：ARQSTLUMNPQ     S₁：ARQSTLUMNPQ

S₂：FDSALQFTGMRA    S₃：RKMFPQDDSLA

这两个模式的偏移量表如下：

L_D(A.Q.T)={(1,1),(2,4)}

L_D(M.PQ)={(1,8),(3,3)}

应该知道每个括号内的第一项是序列号及第二项是偏移量。对应于序列中任何字符的偏移量是该字符离序列始端的距离。例如，(2,4)表示序列是S₂及模式A.Q.T开始位置离序列S₂始端的距离是四个字符。

Q系指查询蛋白质。本发明搜索引擎的一个目标是鉴定D中数据库序列与用户可能提供的任何查询序列Q之间的序列同源性。作为例子，我们将要使用查询Q=JLANQFTLMDPQDLA。此序列具有一些与数据库序列同源的区域。下面我们将要显示它们(相似区域的一对再次用下划线显示)：

Q：JLANQFTLMDPQDLA      Q：JLANQFTLMDPQDLA

s₁：ARQSTLUMNPQ        s₂：FDSALQFTGMRA

因此，搜索引擎鉴定例如以上所示的相似性。如果两个相同长度的区域一个放置于另一个之上时它们之中有数个匹配字符排列在一起，则这两个区域是相似的。下面将要精确地给出相似性的正确评注；现在，有依据说它涉及每个可能的字符对的分数使用。每个这类分数是匹配程度的量度，用于鉴定在生物学上两个字符排列在一起的概率多大。

给定一个模式P，P的“主干”定义为字母{1,0}之上的串，该串是通过在P中将P的每个残基转换为字符‘1’和每个“不必关心字符”为字符‘0’而获得的。例如，模式P=“A..DFE”的主干是串“100111”。主干将模式集合划分为等效的类别，每个类别包含所有共享相同主干的模式。

另一个根据本发明可以采用的概念是模式“密度”。一般而言，该密度描述G(P)(其中G(P)系指多肽语言，它包含可通过将∑中的一个任意残基替代每个“不必关心字符”而从P中获得的所有串)的两个成员之间最小同源性数量，并由两个整数L和W(L≤W)确定：如果P的每个以氨基酸开始和结束并具有至少为W的长度的子串包含L或更多残基，则模式P具有一个<L,W>密度。在每个这类模式中，在模式长度上的残基数量的比例至少是L/W。整数L和W是我们的优选方法的参数，它们的数值控制所完成搜索中允许的相似性数量。这些参数也详细地描述于1998年2月13日提交的系列号为09/023,756的美国专利申请中，该申请提到“TEIRESIAS”算法，它对1997年6月12日提交的美国临时申请系列第60/049,461号要求优先权，后者的内容此处引作参考。注意到，按照定义，一个<L,W>模式具有至少L个残基。

还有，给出一个模式P和一个序列S，S中任何属于G(P)的子串称为P的匹配位置。P的偏移量表包含P的所有匹配位置的第一字符的偏移量。

给出以上定义后，现在我们可以一般地描述根据本发明的用于改进序列同源性检测的优选方案，例如与系统100(图1)相关连。序列同源性检测包括两个不同阶段：信息收集和搜索。

首先，在完成任何搜索之前，挖掘基础数据库D。此挖掘过程也指信息收集或词典组成。在此步骤中，收集所有重要的<L,W>模式并将每个这类模式P与其偏移量表L_D(P)关联起来(用于判定一个模式是否重要的具体准则将要在搜索引擎段落中详细描述)。

第二步是实际搜索。给出一个查询序列Q后，我们鉴定所有由Q匹配的模式P(在过程第一阶段中收集的那些模式)。对于每一个这类P，我们将与P匹配的Q的区域及所有数据库序列中也与P匹配的相应区域(这些区域可以容易地通过偏移量表L_D(P)访问)配成对。最后，配成对的区域在两个方向内延伸和对比并使用一个(用户规定的)突变矩阵评分，然后将最高分数的匹配连同所含对比内容一起报告。

值得指出，此处信息收集阶段是一个对D的一次性计算。所得结果存储于一个文件(图1的模式词典120)中，并且每次在数据库D上执行搜索过程时都使用它。

使用模式来描述相关多肽的动力在于生物学事实。具体地，已知存在不少基本元素(或者是结构特性，例如α螺旋、β链、环等，或者是较大功能单元，例如基序、模块和结构域)，它们是构成蛋白质的结构单元。用于鉴别物种的进化所用关键机理之一是在蛋白质序列中氨基酸位置的突变。然而功能上/结构上的重要区域更为抵制这类突变。因此可以合理地认为能够通过发现以下内容鉴定这类生物学上相关的多肽：(a)在它们一级结构内的保守位置；及(b)增大的可重用性。在我们的术语中，这些性质对应于具有预料不到的高支持值的特性。

然而，重要的是重申，此处描述的本发明搜索引擎方法学可以与其他已知模式词典一起使用。相似地，本发明词典组成方法学可以用于建立模式词典以供其他已知搜索引擎使用。

应该知道，这两种方法学将要分别在下面的段落Ⅱ和Ⅲ中描述。虽然词典组成方法是在搜索方法之前应用的，但为便于阐述，我们按照相反顺序讨论这些过程，从搜索引擎方法开始，后随之以词典组成方法。

Ⅱ．搜索引擎

现在参照图4，显示用于阐述根据本发明一个实施例的搜索引擎方法学的高级流程图。此方法学可能由图1中搜索引擎110采用。搜索引擎的操作可以分成两个不同的阶段：(ⅰ)模式匹配加上链接(块402)；和(ⅱ)评分(块406)。

第一阶段相对于查询序列Q检查Ⅱ中的每个模式P(回忆П系指以上提到的模式词典120)，从而隔离出与Q匹配的所有模式。下面我们描述一个用于完成此“检查匹配”过程的具体算法；然而可以使用任何匹配算法。注意到图4中阶段1的“复杂度检查”(块404)。在某些情况下，有可能一个模式P与查询Q匹配但不希望考虑此匹配。这一例子是所谓“低复杂度”模式。这类模式有时由于生物学序列的特性所引起。低复杂度模式差不多都由相同氨基酸例如模式“A.A..AAA.A.A”所组成，并且由于某些蛋白质具有长的重复氨基酸区域而出现。当然对于同源性检测目的这类模式被考虑为不重要，以及如果把由这些模式引起的任何匹配都加以忽略则更好些。允许系统用户在搜索引擎内将一个“复杂度检查”部件设置为它的“ON”或“OFF”状态，从而让系统用户来决定这么做或不这么做。现在记住当此部件设置为“ON”时，即使P中某些模式与查询蛋白质Q匹配，它们也将被忽略。下面我们将要描述当复杂度检查是“ON”时忽略与Q匹配的模式P的精确条件。

继续看阶段I，每个与Q匹配的模式P生成一个Q与也和P匹配的所有数据库区域之间的局部同源性。这些后者区域可以容易地通过P的偏移量表L_D(P)访问。假定在偏移量i处P与Q匹配，L_D(P)中每个区域(j,k)产生片段(i,j,k,l)，其中l是模式P的长度。这将在下面描述。最后，当匹配过程继续时兼容的片段链接在一起而组成较长片段(下面将要描述兼容片段的评注和链接操作)。在阶段1的末尾，我们具有一个集合R，它包含数据库D中所有与Ⅱ中至少一个模式P匹配的序列以使P也与Q匹配。R中每个序列S伴随用于描述Q和S之间模式引起的同源性的片段。

考虑我们以上引进的例子。查询序列Q=JLANQFTLMDPQDLA与Ⅱ中的两个模式P1=“A.Q.T”和P2=“M.PQ”都匹配。给出P1在偏移量3处与Q匹配以及P2在偏移量9处与Q匹配，这两个匹配产生以下四个片段：

(3,1,1,5)(3,2,4,5)-自L_D(P1)中

(9,1,8,4)(9,3,3,4)-自L_D(P2)中

以及该集合R是：

R={s₁-(3,1,1,5),(9,1,8,4)

   s₂-(3,2,4,5)

   s₃-(9,3,3,4)}

其中R中每个序列S_i带有一个片段表。注意到在此具体例子中没有可能链接。

图4中阐述的搜索引擎方法学的第二阶段对R中每个序列S赋予一个分数。存在着不少方案用于计算一个给定S_j的这个分数。当然每一次开始时都计算由S_j携带的所有片段的分数。每个片段获得一个分数(这些分数称为片段分数)。基于突变矩阵M来完成评分。突变矩阵是20×20实数阵列。这一矩阵的第(i,j)个条目标示在进化中第i种氨基酸已经变为第j种氨基酸的概率。此处为此目的，有依据假定M是来自一个∑×∑-＞R的函数，当其输入为两个氨基酸A1，A2时，它返还一个实数。由于可以使用许多突变矩阵，用户可以使用选项以便选择具体使用的矩阵M。

例如，假定我们使用一元的矩阵M，即对于所有氨基酸都用M(A,A)=1，及对于所有不同氨基酸A，B都用M(A,B)=0。考虑以上集合R中的第一个序列，也即携带片段(3,1,1,5)和(9,1,8,4)的序列S₁。我们显示如何为这两个片段中的第一个评分(另一个以及R套中所有片段都相似地评分)。设想我们已经对比(一个在另一个之下)片段所含长度为5的两个蛋白质区域，即在Q的偏移量3处和在S₁的偏移量1处开始的区域：

ANQFTL    (自Q)

ARQSTL    (自s₁)

然后在所有对比的列上将值M(X,Y)求和而计算该片段的分数，其中X，Y是在任何给定列中对比的两个氨基酸。对于以上片段，此分数是：

M(A,A)+M(N,R)+M(Q,Q)+M(F,S)+M(T,T)+M(L,L)

=1+0+1+0+1+1=4。

以上描述的为片段评分的方案是基本评分方案。也即，系统用户可以设置一系列选项以便修改计算片段分数的方式。例如，如果系统参数“extend”(它是一个整数及将要在下面描述)已经设置为大于零的值，评分时不但考虑由该片段描述的蛋白质区域，而且也考虑两个区域左侧和右侧的“extend”氨基酸的区域(评分过程完全和以上描述的一样进行，只是现在考虑较长区域)。此外，如果已经设置“gapped_alignment”选项，则在对比延伸的区域(即基本片段的左侧和右侧的那些区域)时，我们还使用空位以便使对比分数最大。

在以上过程的末尾(不取决于所使用的不同评分方法)，为集合R中的每个片段计算片段分数。然后在评分阶段的最后步骤中利用这些片段分数，即将Q与R中所有序列S_j之间的相似性量化。此量化过程通过对R中每个S_j赋予一个分数而完成；此分数称为S的“序列分数”(以将它区别于片段分数)。理想情况下，序列S_j的序列分数越高，则此S_j与Q越相似。

在给S_j评分过程中，我们只考虑S_j所携带片段的片段分数。但还是有数个选项。在最简单的情况下，S_j的序列分数定义为S_j片段的所有片段分数之间的最大值。下面描述第二个更复杂的方案。此处首先为待评分的序列S_j建立一个有向图。此图的顶点是所有S_j所携带的片段。每个顶点被赋予对应于该顶点片段的片段分数。如果

i＜=i＇和k＜=k＇，

即如果由两个片段(区域Q[i..i+l-1]和Q[i＇..i＇+l＇-1])所描述的两个查询区域的相对顺序与由两个片段(区域S_j[k..k+l-1]和S_j[k＇..k＇+l＇-1])所描述的两个查询区域的相对顺序相同，则自片段(i,j,k,l)作一条边至片段(i＇,j＇,k＇,l＇)。与这些顶点一样，也对每条边赋予一个分数，以表示查询中区域的位移(即i＇-i之差)相对于S_j上区域的位移(即k＇-k之差)是如何正规。位移之间差别越大(即｜(i＇-i)-(k＇-k)｜之值)，则该边的分数越小。在建立该图之后，我们能够应用任何最长路径算法来鉴定具有最高分数的路径(路径分数定义为路径中所有顶点和边的分数之和)。此分数然后成为S_j的序列分数。

以上已经描述了计算片段和序列两者的分数的数种方式。一般而言，能够使用任何其他“生物学上合理”的评分方案以作为替代。

现在参照图5，6和7，将要描述搜索引擎方法学400的模式匹配、链接和评分过程的更具体例子。再次，在由搜索引擎110实施的搜索阶段中，提供查询蛋白质Q给该系统，并且鉴定相似于Q的数据库序列S∈D及报告回至用户。该搜索阶段利用一套通过挖掘输入数据库D而获得的模式集合П。为示例目的，有依据假定П是一组在以上“定义”段落中描述的形式的模式<L,W>的集合。每个模式P∈П由其偏移量表L_D(P)所伴随及在D中具有至少为K_min的支持值。数值L、W和K_min是我们优选方法的参数，以及下面将要在“词典组成”段落中描述设置它们的方式。

当一个查询序列Q提供给系统时，要做的第一件事情是找出所有由Q匹配的P∈П。通过使用一种散列技术可以很快地完成此点，该技术出现在D.Gusfield的“串、树和序列上的算法：计算机科学和计算生物学”(剑桥大学出版社，62-63,1997)。更具体地，对于在Q之内的每个位置，我们生成W个散列值，每个在该位置处起始长度为2、3、…、(W+1)的子串都有一个值。对于每个这类子串，相应的散列值只决定于子串的第一个和最后的字符，以及决定于这两个字符中间的残基数量。

图5提供给定查询序列的过程例子。在此例中，显示了为在序列Q的位置6处起始的W=4的子串生成的散列值。子串s使用的散列值是：

H(s)=((av(first_char)-av(＇A＇))+(av(last_char)-

av(＇A＇))^*26)^*W+gap

其中av(c)是字符c的ASCII值，first_char和last_char分别是s的第一和最后字符以及gap是s的第一和最后字符中间残基数量。注意到，由于<L,W>密度的限制，gap始终小于W。

对应于具体值h的散列条目包含查询序列Q的所有偏移量p，以使在p处起始的子串(长度最多为W+1)的散列值为h。图6给出为一个具体查询序列生成的散列表例子。图6中显示为序列Q=AFGHIKLPNMKAMGH生成的散列表的快照。为标记表的条目，我们不用实际数字散列值而使用一个描述所有具有具体散列值的串的模式。每个散列条目指向一个偏移量表。表中每个偏移量标志着Q中与有关的散列条目相关的子串的开始。

为检查一个模式P∈П是否由Q匹配，我们使用一个计数器阵列C[1..｜Q｜]，其长度等于Q的长度。初始地，阵列的每个条目设置为0。从P中偏移量为1处开始，我们找出在P内所有对应于一个残基的偏移量j，并排除对应于最后残基的偏移量。对于每个这类j，使F为在j处起始并就只包含两个残基的最短子串。使OL标志Q中由对应于F的散列表条目所指向的偏移量表。如果OL不空，则对于每个偏移量p∈OL，计数器C[p-j+1]按1增量。如果模式P就只包含n个残基，则在此过程末尾，只当Q在偏移量j处与P匹配时，计数器C[i]才具有值(n-1)。以上描述的匹配技术的一个优点在于它通常需要的时间与查询序列Q的长度几乎成直线关系，并且只决定于模式P中残基数量。

一旦发现模式P∈П由在偏移量i处开始的Q的子串所匹配，我们需要将该Q的子串及也与P匹配的所有数据库区域关联起来。只须扫描只包含这些区域的偏移量表L_D(P)，即可容易地完成此点。更具体地，每个条目(j,k)∈LD(P)标示在第j个数据库序列S_j的偏移量k处起始的子串是G(P)的一个元素。查询序列Q与数据库序列S_j之间的局部相似性然后登记为一个四联组(i,j,k,l)，称为一个片段，它与S_j相关连。数值l=｜P｜是局部相似性的长度。

有时，两个不同的与Q和数据库序列S_j两者都匹配的模式P和P’对应于Q和S_j之间的相同的局部相似性。这一情况的例子阐述于图7中。在这类情况下，个别的对应于两个模式的片段必须链接成一个。具体地，两个与S_j相关连的片段(i,j,k,l)和(i＇,j＇,k＇,l＇)只当满足以下情况才称为兼容，即：

k＜=k＇和k+1+w len＞k＇和k＇-k=i＇-i

其中w_len是由用户规定的整数参数；w_len允许链接不相交的片段，只需一个片段在另一个末尾之后不超过w_len个位置就开始即可。将(i,j,k,l)和(i＇,j＇,k＇,l＇)链接在一起所得片段为：

(i,j,k,max(l,k＇-k+l＇))

作为找出由Q和S_j两者匹配的模式P∈П的结果，每当一个新片段与一个数据库序列S_j相关连时，即将兼容片段链接。如果已有片段与S_j相关连，及与新到的片段兼容，则将新的和已有的片段的相关对去除，并由它们的链接结果替代。

已经鉴定Q与数据库序列之间的所有局部相似性之后，剩下的任务是评价这些相似性。通过对每个与至少一个片段相关连的数据库序列S_j赋予一个分数(使用一个用户规定的评分矩阵)，即可完成此事。有数个选项可用于评分功能。此处给出本发明原理之后，熟悉技术的人知道其他评分方式。

如上所述，一个方案是个别地为S_j的每个片段评分并将这些分数中的最高分赋予S_j。可用以下两个方法中的任何一个为片段(i,j,k,l)评分：

不允许空位：在此情况下，根据由片段所含无空位对比来计算分数，即查询的区域Q[i,i+l-1]和序列的区域S_j[k,k+l-1]的对比。此外，用户可以有选项来设置可变的＂extend＂以便“围绕”片段延伸对比。如果此值大于0，则根据区域Q[i-extend,i+l-1+extend]和S_j[k-extend,k+l-1+extend]的无空位对比来计算分数。

允许空位：此选项只当extend＞0时方才可用，并允许在对比区域内存在空位，从而获得围绕该片段的区域的更精细评分。

如上所述，也提供其他评分选项，考虑到与当今正被评分的数据库序列S_j相关连的片段相对顺序。在如上所述的个别地为每个片段评分之后，一个方案建立一个有向加权图，如图8所示。此图的顶点是与S_j相关连的片段，并当：

i＜=i＇和k＜=k＇时，

在片段(i,j,k,l)和(i＇,j＇,k＇,l＇)之间存在一条有向线。每个顶点被赋予一个等于相应的片段的分数的权值，而每条边E根据以下内容加权：(a)两个片段如何靠近，即(i＇-i-1)的值；及(b)两个片段之间的位移如何正规，即(i＇-i)与(k＇-k)差别多大。此图内的一条路径的分数是该路径的所有顶点和边的权值之和。然后计算具有最大分数的路径并将该分数赋予S_j。

现在参照图9和10，相应的流程图总结性地阐述由本发明搜索引擎模块所完成的两个阶段的实施例。图9阐述匹配和链接阶段的实施例900，而图10阐述评分阶段的实施例1000。

在图9中，假定数据库D中每个序列S_j具有片段SegL(S_j)的一个关联表。初始地，所有这些表都是空的。还有，集合R初始地是空的。当由图9流程图所描述的计算进行时，R由序列S_j所占用。当这一序列插入R时，它带有其片段表SegL(S_j)。

因此，对于Ⅱ中每个模式P(块902)，搜索引擎完成以下操作。在步骤904，搜索引擎判定P是否与Q匹配。如果不是，则进至词典中下一个P。如果是，则在步骤906搜索引擎判定用户是否已允许复杂度检查部分进行工作。如已允许，则在步骤908引擎判定P与Q的匹配是否为低复杂度匹配(下面将更详细地阐述)。如果是，则引擎进至词典中下一个P。如果不是，则对于P与Q匹配的所有偏移量i(块910)和对于L_D(P)中的所有(j,k)(块912)，引擎完成以下操作。在步骤914，它将片段(i,j,k,｜P｜)与SegL(S_j)中任何兼容片段链接起来。然后引擎将结果加入SegL(S_j)中。

在步骤916，引擎判定S_j是否在R中。如果是，则引擎回至步骤914。如果不是，则引擎将S_j和SegL(S_j)加入R中。为所有P与Q匹配的偏移量i和L_D(P)中的所有(j,k)完成步骤914至916。对模式词典中每个P重复整个过程(步骤904至916)。

已经完成匹配和链接过程后，搜索引擎可以执行图10中的评分操作。引擎在步骤1006中为R中所有序列S_j(块1002)和S_j中所有片段s(块1004)计算s的片段分数。然后引擎在步骤1010中为R中所有序列S_j(块1008)计算S_j的序列分数。最后，在步骤1012中，引擎报告R中得到最高分数的S_j以及由它们相应的序列分数所含的局部对比。

再参照图4，如前所述，搜索引擎模块110可能包括一个复杂度检查部分(例如图9的步骤906)。复杂度检查部分负责将由于低复杂度区域而生成的局部同源性去除，首先，分两个阶段实行低复杂度检查，这些都是在词典建立阶段(“词典组成”段落)以及搜索阶段(此段落)两个期间。

在词典建立阶段期间，低复杂度区域以两种方式处理。首先，当寻找输入数据库中的模式时，我们去除了(即从输入中去除)所有由L或相同氨基酸的更多连续出现所组成的蛋白质区域(L是我们在词典建立阶段期间设置的一个整数；这里完全可以假定它具有某个固定值)。这考虑如以下用下划线显示的低复杂度区域(其中那些点标示在所阐述串的左侧和右侧还有更多氨基酸)：

……ASDFHRTYIUSFFFFFFFFFFFFFFFFFFAKJRBVCJ……

然而，这只是低复杂度区域的一个情况。还存在很多。例如考虑以下区域的下划线部分：

GFWRETIOJIFPAPAPAPAPAPAPAPAPAPAPAPAJSHDGF……

为处理这种类型的区域(例如具有一般化重复组成部分者)，我们也去除给定模式P的所有重叠出现部分。换言之，如果模式P在偏移量k₁和k₂处(其中k₂＞k₁)与数据库序列S_j匹配及k₂-k₁小于P的长度，则不在模式P的偏移量表中放置这两个偏移量中任何一个。例如，在以上所示区域中，模式“P.P.PA”具有长度为6，它在偏移量12和14处出现(在其他位置之间)，也即在重叠位置出现，因为14-12=2而2＜6。

现在在搜索引擎阶段期间，我们已有两种方式扑获和去除低复杂度同源性。第一个是以上所给例子的一般化。简而言之，我们去除所有不是“语言上丰富”的模式，即它们表现出一个具体氨基酸的过分表示。为此目的，我们允许用户设置参数V的值(0与1之间的实数)。只当与查询序列Q匹配的模式P的可变性(variability)v(P)不大于值V时，才可以进一步加以考虑。具体地，对于每个模式，我们定义其可变性为：

即使通过以上描述的可变性测试，还要进行第二层检查。此第二层试图扑获低复杂度的更细概念。为了解它如何工作，考虑以下例子。我们假定查询蛋白质是以下简单串：

Q=FRGDSAAABBBBAABBSJIEKL

以及考虑模式P=“A…B..AB”。该模式在偏移量7处与查询匹配，如下所示：

A…B..AB

FRGDSAAABBBBAABBSJIEKL

匹配区域及其最接近的周围是低复杂度区域(它只由‘A’和‘B’组成)。然而该模式P的可变性只有0.5。为处理此字符的低复杂度区域，我们允许用户定义整数“margin”和“min_m”(其中min_m＜=2^*margin)以及一个百分比perc。然后检查在所考虑的模式(此处为模式“A…B..AB”)的实际匹配位置的左侧margin个字符内和右侧margin个字符内(此处查询的偏移量7)的近似匹配。如果将模式P放置于给定偏移量处，至少有模式中的正规字符中perc％个字符与查询的基础字符匹配，则认为模式P在该查询的偏移量处近似地匹配。例如，如果perc=75％，则模式“A…B..AB”在偏移量6和8处近似地与Q匹配，如下所示：

A…B..AB

FRGDSAAABBBBAABBSJIEKL(在偏移量6处)

A…B..AB

FRGDSAAABBBBAABBSJIEKL(在偏移量8处)

其中在这些偏移量中的每一个偏移量处，75％的模式正规字符(即4个中的3个)与相应的查询字符匹配。已经给定参数的margin，min_m和perc后，现在可以说在此检查层次上何时查询与数据库区域之间获得的局部同源性是低复杂度的。考虑一个在偏移量X处与查询Q匹配及在偏移量Y处与数据库序列S匹配的模式P。如果以下情况之一是真实的，则该匹配考虑为低复杂度：(ⅰ)该模式在X的左侧和右侧2^*margin个字符的至少min_m个字符内近似地与Q匹配；或(ⅱ)该模式在Y的左侧和右侧2^*margin个字符的至少min_m个字符内近似地与序列S匹配。

Ⅲ．词典组成

如前所述，在优选实施例中，词典组成方法学在搜索引擎从用户接收一个查询序列之前完成。这是因为，再参照图1，搜索引擎模块110最好利用由词典组成模块130所组成的模式词典120。词典组成模块130实施本发明数据库处理方法学以便组成模式词典(或生物学词典)，这将要在下面描述。然而，如前所述，模式词典120也可由此处所阐述的以外的其他搜索引擎使用。也即，现有搜索引擎可能利用从根据本发明的源数据库中挖掘出的模式。当然，根据优选实施例，假定根据此处阐述的本发明方法学所组成的模式词典将要由也在此处阐述的本发明搜索引擎使用。

模式词典组成阶段(也称为信息收集阶段)确定所考虑的数据库D中找到的所有重要的<L,W>模式的集合П。实质上这是一个数据挖掘过程，其中利用D以便试图找出D的序列之间的隐含关系。此想法集中于那些没有预料的关系，并且依靠这些性质假定它们具有生物学关系。为此目的，将要通过模式在D内的支持值阐述模式的重要性。更具体地，将要定义数量K_min(最小支持值)，以使每个其支持值至少为K_min的模式可以说明为统计学上重要的。所有这些模式(也有个别例外，它们不容忍最小支持值的要求)将要包括于集合П内，作为搜索阶段的输入。

回忆K_min的概念首先在上面的“定义”段落中引入。“密度”的概念也已引入。回忆密度描述G(P)的任何两个成员之间的同源性的最小数量(其中G(P)系指多肽语言，它包含所有可以通过用来自∑的任意残基替代“不必关心字符”而从P获得的串)，并且由两个整数L和W定义(L＜=W)：如果模式P中每个以氨基酸起始和结束并具有至少为W的长度的子串包含L个或更多残基，则该模式具有密度<L,W>。再次，这些参数描述于以上包括的于1998年2月13日提交的美国专利申请系列号09/023,756中，它提出“TEITRESIAS”算法。

虽然本发明优选方法学在组成模式词典Ⅱ中利用参数L和W，但应该知道可以采用其他已知用于确定一组序列的任何两个成员之间的同源性的最小数量的技术。

设置参数L、W涉及对一些有时矛盾和互连的因素的考虑。例如，比例L/W描述在搜索阶段中D中查询序列与蛋白质之间允许的同源性数量。小的L/W允许检测弱相似性。由于数个值对(L,W)会产生相同比例L/W，那么什么应该是L和W的正确设置值呢?选用大量的L通常会使信息收集阶段的运行时间过长(除非L/W接近1)。此外，选择大的L值将忽略只有数个氨基酸的弱模式，而它们往往处于感兴趣的模式之中(即通常被现有的相似性搜索工具所忽略)。另一方面，选择太小的L(例如2或3)可能无用，由于在此情况下在输入数据库D中的具有L+i个残基的<L,W>模式的分布与具有D的氨基酸组成的随机数据库中的相应分布并无明显区别。在大多数通常情况下，应该知道可以完全任意地选择L、W和K_min的值。然而，为实际上保证所发现的模式是远在统计学噪音水平之上，我们使用一个统计学框架(即一种设置以上所述参数的方法)来增强模式发现过程。

为使以上观点更明显，考虑图11，它将已知为SwissProt Rel.或SP34的测试数据库中的模式分布(见A.Bairoch和R.Apweiler的“SWISS-PROT蛋白质序列数据库及其在1998年的补充TrEMBL”，(Nucleic Acids Res26：38-42,1998))与相应的随机分布进行比较。图11阐述一种具有SP34中给定主干结构的模式的分布(该分布由符号“o”表示)以及与相同主干的随机分布(该分布由符号“+”表示)互相比较。主干的概念首次在以上“定义”段落中引入。曲线上的一点(X,Y)表示正好有(给定主干结构的)Y个模式而这些模式中的每一个具有支持值X，即它正好由X个不同的数据库序列所匹配。此处显示的结果是使用SP34的一个“清除”版本(数据库的清除将要在下面阐述)而获得的。对于SwissProt，我们计算正好具有L个残基的模式<L,W>的支持值(对于图11中所示L，W的值)。然后将结果列表以便为每个可能的主干建立一行：对应于一个给定主干B的行的第i列标示在SwissProt中具有支持值i的模式(该主干结构的)的数量。随机分布是对于SwissProt中N=2000个随机地弄乱的版本使用完全相同的方法所获得的(图13阐述弄乱过程，用于获得每一个弄乱的版本)。在此情况下，给定主干B的行是通过将对应于所有2000个表中的B的行平均而得。其结果是，第i列给出具有主干B的模式的平均数的足够正确的评价，这些模式出现在具有SwissProt的残基组成的随机数据库内的正好i个序列中。在图11中，我们画出所选主干的SwissProt结果与相同主干的平均分布的曲线。虽然呈现的结果涉及具体主干，但如使用其他主干将不会有质的变化。

注意到我们是使用2000个采样点(输入数据库的随机弄乱的版本)。这只是为了阐述的目的。原则上实际中的采样点可以任意地设置。一般而言，当这些点的数量增大时，所获得的评价能够更正确地收敛于它们的真实值。给出对于待计算的评价的任何希望的信任度，就能够使用标准统计学理论以确定应该使用多少采样点。

如图11所示，只当L变为5或更大时，我们开始区别SwissProt中与随机数据库的组成偏倚(compositional bias)(以模式表示)。一般而言，L的值取决于基础数据库D的大小：数据库越大，则此值应该越高。所示对SwissProt的结果是使用L=6的值而获得的。对于W，我们选择值15，以使比例L/W(即最小允许同源性)为40％。

已经设置L和W值后，剩下来的事是确定最小支持值K_min应该为何值。由于每个较大模式包含至少一个正好具有L个氨基酸的子模式，我们只集中于正好具有这么许多L个残基的模式。一个方案是选择K_min以使在K_min中或更多不同序列中出现一个模式的概率为小。但是对图11(d)更细致地观看后，发现此方案可能太严格。具体地，考虑一个支持水平K=15。随机分布标示可以只靠偶然性预期一个和两个具有支持值K的模式之间。因此，根据以前所述准则，在SwissProt内具有支持值15的模式被认为是不重要的。然而这两个分布在该支持水平上具有给人印象深刻的差别。具体地，虽然K=15处的随机分布平均值大约为1.5，而SwissProt内有大约180个具有支持值15的模式。

因此，似乎如果考虑到在隔离中的模式的概率，其结果是去除许多根据以上分布是在噪音水平之上的模式。此观察提示我们应该使用不同的重要性准则。

现在参照图12至15，我们提出一个阐述用于确定重要性准则的优选方案的流程图。也即，我们提供一种方法学，用于计算K_min的值。给定K_min值后，通过在其中包括源数据库D中具有至少该值K_min作为支持值的所有模式，从而组成模式词典Π。因此，应该知道图1的词典组成模块130可以完成图12至15中阐述的过程。

通常在我们的方案中，我们不是寻找个别模式，而是一起考虑一个具体主干结构的所有模式。更具体地，对于任何给定主干B和一个基础数据库D，使N_B,K为：

N_B,K=具有主干B并在D内具有支持值K的模式数量。

还有，使X_B,K为对应于N_B,K的随机变量(在D的所有弄乱的版本的空间内定义的)。然后最小支持值K_min是满足以下不等式的第一个数值K：

其中threshold是一个用户规定的概率，隐含出自以上不等式的最小支持值K_min的信任度水平。较小阈值导致K_min的较大值及对于事实上所选模式的较大重要性。

因此，作为用于确定重要性准则K_min的过程的输入，我们有源数据库D，整数参数L和K，一个标示物种数量的整数N和其值在0与1之间的实数threshold。当然，作为过程输出，我们得到整数K_min以使D中具有支持值K_min或更多的所有模式都是统计学重要的，从而包括于模式词典中以供接收用户查询时搜索之用。

以下流程图的解释使用不同标志，它们中某些已在上面引入过。然而，为简明起见，应用以下定义：给定任何模式P，P的主干B(P)定义为当用‘1’替代P的每个正规字符及用‘0’替代每个“不必关心字符”时所得{1,0}上的串，例如如果P=A..F.G..R，则B(P)=100101001。如果B是一个任意主干及P是一个模式以使B(P)=B，则我们说P就是B模式。N_B,K即称为在D之中具有支持值K的B模式的数量，Xⁱ_B,K是第i个随机数据库中具有支持值K的B模式的数量。m_B,K是所有Xⁱ_B,K的平均值及S_B,K是所有Xⁱ_B,K的方差。应该知道，由于我们没有随机变量X_B,K的分布的分析描述，我们采用标准采样技术。因此，对于给定数据库D，我们能够为随机变量X_B,K的平均值m_B,K和方差S_B,K两者实验地计算正确的评价。

初始地参照图12，总过程1200开始时在D上执行TEIRESIAS算法(例如如上所述及在以上所包括的于1998年2月13日提交的美国专利申请系列号09/023,756中所描述者)并在步骤1202中计算N_B,K。虽然优选TEIRESIAS算法，但应该理解可以使用其他常规技术计算N_B,K。然后对于i=1至N(块1204)完成以下步骤。

在步骤1206，生成一个随机数据库R_D_i。此步骤将要在图13中进一步解释。如过程1300中所解释，R_D_i(块1302)是通过对D中每个序列S(块1304)计算S中字符的随机排列(块1306)而生成的。S中字符的随机排列称为S＇。S＇加至R_D_i中(步骤1308)。重复该过程直至处理完D中所有序列S(块1310)。因此，R_D_i包括所有随机排列S＇。回至图12，在步骤1208中，在R_Di上运行TEIRESIAS以便计算Xⁱ_B,K。步骤1206和1208是用于所有i的，也即使用直至i=N(块1210)。

然后，对于每个B,K(块1212)我们使用Xⁱ_B,K以计算m_B,K和s_B,K。此步骤将要进一步在图14中解释，如过程1400中所示，首先将s_B,K设置为0(步骤1402)。然后，对于i=1至N(块1404)，作为s_B,K与Xⁱ_B,K之和，计算s_B,K(块1406)。对于所有i重复该过程(块1408)，最后将s_B,K除以N以计算s_B,K的平均值(块1410)。然后在步骤1412至1420中计算方差m_B,K。首先在步骤1412中将m_B,K设置为0。然后，对于i=1至N(块1414)，作为m_B,K与(Xⁱ_B,K-s_B,K)²之和，计算m_B,K(块1416)。对于所有i重复该过程(块1418)，最后将m_B,K除以N以计算方差m_B,K(块1410)。

回至图12，在步骤1216中，使用m_B,K和s_B,K计算p_B,K。此步骤将要进一步在图15中解释。如过程1500中所示，在步骤1502中定义一个实数C，以使：$$ >>>N>>B>,>K>\; >=>>(>>m>>B>,>K>\; >+>1>.>96>>>S>>B>,>K>\; >>>N>>>(>1>+>>>>1>.>96>>>2>N>>>>)>>>>)>>+>C>>(>>>S>>B>,>K>\; >>1>+>>>>1>.>96>>>2>N>>>>>>)>>>$$

其中N表示物种数或试验次数的具体值，例如2000。因此在步骤1504中，p_B,K作为1/C²的值计算出来。应该理解，p_B,K是概率Pr[X_B,K＞N_B,K]的上限。总而言之，我们使用X_B,K的物种平均值和方差为手头N_B,K的值计算C。我们知道常数C与统计学中众所周知的Chebychev不等式相关连。注意到使用95％的可信度水平来计算常数C，然而这不是要求。也即，任何其他值都可使用。

回至图12，为每个B，K重复步骤1214(图14)和1216(图15)。然后在步骤1220中确定K_min为最小的K，以使max_B{p_B,X}≤threshold。在下一段落(SwissProt Rel.34)中呈现的测试情况下，threshold的值如此选择以使K_min=15，即在此支持水平上，只能偶然地遇到1.5个给定主干结构的模式。这里有一个折衷：我们愿意允许少量模式导出的局部同源性(偶然性的结果，如以上的1.5个模式)以便能扑获在SwissProt内出现的相同支持水平上由其他模式隐含的多得多的统计学重要相似性。

在下一段落中提供某些实验结果之前，我们首先解释在完成本发明词典组成方法学之前清除数据库的概念。此过程阐述于图16中，及也可由图1中词典组成模块130实施。数个数据库包含几组高度同源的序列(例如血红蛋白α链蛋白质)这些组不但通过引入大量模式而使模式发现过程慢下来，而且它们也能虚假地抬高一个模式的重要性。当模式在一个高度同源的序列家族中多次出现而很少在其外部出现时，就会出现此种情况。

为处理这些问题，可在模式发现过程开始之前“清除”数据库D。如图16所示，清除过程1600涉及将高度相似的蛋白质鉴定和合并在一起(步骤1602)。如果在最佳地对比两个序列后，较短序列的位置中有X％(例如50％)与较长序列中的位置完全相同，则将这两个序列放置于相同的组中。所得组称为冗余组。将要在其上完成信息收集过程的集合D＇由以下组成：(a)D中那些被发现为与其他蛋白质不是足够地同源的序列；及(b)冗余组的每一组中最长的序列(步骤1604)。最后，使用TEIRESIAS算法分别地处理每组冗余组(步骤1606)，收集模式直至该组中所有序列都与这些模式中至少一个匹配。此方案保证即使包含多结构域蛋白质的组也能通过为每一结构域生成至少一个模式来正确地加以处理。值得指出处理冗余组所得模式通常都很密(残基的数量比“不必关心字符”的数量多得多)和长。这是组序列的高度同源性的后果。对于这些模式，我们允许在搜索阶段中近似的匹配。

Ⅳ．实验结果

在此段落中我们讨论与本发明优选实施例相关连的实验结果。也即，根据作为试验数据库的SwissProt Rel.34，实现以上所述的本发明的词典组成(信息收集)和搜索引擎方法学两者，从而生成以下结果。信息收集阶段中发现的模式的定量和定性描述在下面第一子段(B)中给出，其中分析这些模式为SwissProt完成的覆盖及评注它们之中最经常出现者。在下面第二子段(B)中，我们在数个查询序列上呈现搜索阶段结果。

A．信息收集

SwissProt的处理从如前段中所述清除它开始。此过程的结果详细地解释于图17中。在SwissProt上的清除过程生成9,165个高度相似序列的冗余组。清除过的数据库(信息收集阶段将要使用的那个数据库)由下法组成：从原始输入中去除高度相似序列，然后通过在其中增加每个冗余组中的最长序列而增强所得的集合。

清除过的数据库可用后，TEIRESIAS在其上应该做的全部事情是设置参数L，W和K_min的值。如前所述，我们使用设置值L=6和W=15。此外，在此处所报告的结果中，在计算方差时，我们选择threshold值为10^-11和可信度水平为95％。为这些设置值计算的K_min值成为15。在清除过的数据库上使用以上规定的L，W和K_min的值而运行TEIRESIAS，生成了一个具有534,185个模式的集合П(模式词典)。

挖掘清除过的数据库只是信息收集阶段的第一步。还需要在9,165个冗余组上应用发现过程。再次，我们使用TEIRESIAS以便处理每组收集足够的<6,15>模式的这类组从而保证组中每个序列与至少一个模式匹配。这些模式然后加至集合П以便组成将要由搜索阶段使用的模式的最后集合П。图18提供关于在整个SwissProt Rel.34上由这些模式完成的覆盖的信息。由一个模式覆盖的数据库区域正好是那些与该模式匹配的子串。注意到，对于密和长的模式(大部分来自冗余组的处理结果)，我们允许近似的匹配，其中模式的“大部分”(具体地，模式的残基的80％)由一个区域匹配。值得指出，大部分未覆盖的序列是片段。更具体地，只有231个具有大于50的大小。图19给出П中模式的以下特性的分布：(ⅰ)SwissProt Rel.34模式的长度；(ⅱ)氨基酸或残基的数量。

如图18中例子所示，已经达到搜索阶段所想达到的成功关键目标(也即SwissProt的好覆盖)。余下需要回答的问题是所发现的模式是否为生物学上相关的。为努力解决此问题，我们分析了这些模式中最经常出现的模式。所得评注示于图20。自此分析可以明显地看出(至少对于那些检查过的模式)，模式发现过程鉴定出在生物学上重要的序列特征。

图20阐述具有最高支持的100个模式。只要可能，一个分类中的模式彼此已经对比。小写斜体字是为方便而用，它们用作以下方括号表达式的占位符号：a：[STGDAR],b：[STGDK],c：[STGDKY],d：[STGK],e：[GASMDL],f：[GISETV],g：[LIVMFY],h：[LIVMF],i：[LIVMA],j：[LIVMC],k：[LIVMF],l：[ILVMF],m：[QKCS],n：[KRQA],o：[IVTNF],p：[QKCASN],q：[QKIAGN],r：[RKAHQN],s：[KRQNE],t：[KRQMN],u：[LFYIMS],v：[AGSPE]。一个括号标示一个可由括号中残基中任何一个所占用的位置。

应该知道并不是所有发现的模式表现出这种清楚的功能特异性。它们之中有一些对应于那些为在功能上评注一个蛋白质传统地被认为是不重要的区域(例如环、卷曲螺旋、跨膜区)。然而，有些时候甚至这类弱相似性也能为蛋白质区域的特征化提供有用提示。我们已经实施两个允许使用此潜能的机理。首先，向用户提供所有由查询序列匹配的模式的表。在大多数情况下，一个有专长的用户能够鉴定哪些模式具有生物学重要性。选择具体模式可以使评分精细化，从而只集中于由此模式覆盖的数据库区域。其次，当基础数据库包括不同数据库序列区域的评注时，此评注与模式一起用于提取有用信息。下一子段中给出这两个机理的使用例子。

B．搜索

为阐述搜索阶段(并阐述如何使用它)，我们选择两个查询序列。第一个是很好地研究和评注过的核心组蛋白3蛋白质(SwissProt：H31_HUMAN)，而第二个是尚未特征化的来自詹氏甲烷球菌(Methanococcus Jannaschii)的ORF(SwissProtID：YZ28_METJA)。

H31_HUMAN

由于核心组蛋白在细胞内包装DNA的中心作用，它们是广泛研究的对象。这些小蛋白质富含带正电荷的氨基酸，它们有助于与带负电荷的DNA双螺旋结合，见J.D.Watson,N.H.Hopkins,.W.Roberts,J.Steitz和A.M.Weiner的“基因的分子生物学”，The Benjamin/Cummings Publishing Company，第四版，1987。四个核心组蛋白(H2A、H2B、H3和H4)结合在一起，成为一个八聚体结构(相似于园柱形楔形体)，它为146个bps长的DNA片段提供基底以供围绕，从而在细胞染色质内组成核小体复合物。

SwissProt Rel.34数据库包含33个序列，它们被评注为组蛋白3，其中有在人类中找到的核心组蛋白3蛋白质H31_HUMAN。使用我们的同源性检测工具将此序列搜索的最高评分结果列于图21中。挨着每个序列给出序列与H31_HUMAN之间的最高评分局部对比的相似性分数。图21中提到的分数是使用PAM 130矩阵得到的(见M.O.Dayhoff,R.M.Schwartz和B.C.Orcutt的“蛋白质中进化变化的模型”，Atlas of Protein Sequence and Structure,5：345-352,1978)，以及将其最高评分片段的分数赋予数据库中每个匹配的序列。

SwissProt Rel.34中的所有33个组蛋白3被正确地鉴定为与H31_HUMAN同源。此外，发现数个其他模式(YB21_CAEEL,CENA_HUMAN,CSE4_YEAST,YL82_CAEEL,CENABOVIN,YMH3_CAEEL)具有与H31_HUMAN的广泛局部相似性。查看这些蛋白质的评注，可以发现它们是已知的相似组蛋白-3的蛋白质。还有一个最后注解，在SwissProt Rel.34中出现的H3_NARPS(一个已知组蛋白3)只是一个片段，这是为什么在结果表中它得分最低的原因。

图22给出为查询序列H31_HUMAN生成的对比的所选视图(最高和最低评分两者)。在图22中显示H31_HUMAN与一个高度相似的(H3_YEAST)和中等相似的(CENA_HUMAN)蛋白质的局部对比。对于每个序列，报告有不少局部对比。在每个这类相似性中，相关的查询(＂Query＂)和数据库序列(＂Seq＂)区域一个在另一个之下地列出，在它们之间是所得一致的区域。我们使用‘+’来标示化学上相似的氨基酸。

YZ28_METJA

在某种重要性上讲，H31_HUMAN是一个容易的测试情况，因为数据库包含数个与其高度同源的序列。一个要问的有趣问题是当出现“边界线”序列时，也即不知其同源性的序列时，我们的方法学如何发挥作用。在解决此问题的努力中，系统中出现一个尚未评注的序列YZ28_METJA，一个具有詹氏甲烷球菌的基因组中1272个残基的开放读取帧。

图23中阐述与查询序列YZ28_METJA一起提出时由我们的系统产生的最高评分对比。所用突变矩阵是PAM130。

为了YZ28_METJA的功能评注，以上提到的结果不是十分出色，因为数据库命中的内容涉及十分分散的蛋白质：前两个(NTNO_HUMAN,NTNO_BOVIN)是依赖于钠的降肾上腺素输送器，而最后一个(KAPL_APLCA)是一个激酶。

带着这些问题，我们进而检查YZ28_METJA与数据库序列之间的更密切的相似性。对于此分析，个别地仔细检查每个与YZ28_METJA匹配的模式。应该理解，本发明的搜索阶段允许用户选择任何由手头查询序列匹配的模式并集中于单独由具体模式引起的局部对比而不管所用其他模式。此特征应用于由YZ28_METJA匹配的模式中的每一个。意图在于发现任何这种模式是否专用于一个具体蛋白质家族，从而给出有关YZ28_METJA的功能性的线索。

其结果是，存在三个专用于激酶家族的模式(即模式"Y..S..I...DLK","NIL......IKL"和"I.H.DLK......D")。图24描述它们之中的第一个所产生的最高评分对比之间的几个，即模式“Y..S..I...DLK”所引起的查询序列YZ28_METJA的最高评分局部对比。所用突变矩阵是PAM130。图25包含一个所有包含该特定模式的数据库序列的完整列表。图26和27给出其余两个模式的相应的列表。图28提供(a)由YZ28_METJA匹配的所有模式的分布和(b)这三个激酶专用模式所覆盖的区域的图形表示。

模式“Y..S..I...DLK”在SwissProt内生成24个命中。所有这些蛋白质(NABA_RAT除外，它是一个钠/胆汁酸双传输器)都评注为蛋白质激酶(它们中的两个，KD82_SCHPO和KKKI_YEAST被表征为推定的/可能的激酶)，其中大多数属于丝氨酸/苏氨酸激酶家族或显示与该家族的相似性。此外，“Y..S..I...DLK”不但属于这些蛋白质的激酶结构域，而且实际上包含该结构域的活态位置(氨基酸D)。

在图25中显示包含模式“Y..S..I...DLK”的SwissProt Rel.34序列。所有它们都评注为蛋白质激酶或可能的/推定的蛋白质激酶(差不多全部丝氨酸/苏氨酸变种)。唯一的例外是评注为钠/胆汁酸双传输器的蛋白质NABA_RAT。

图26显示三个模式中的第二个"NIL......IKL"所得相似结果。在此情况下，数据库命中数量是34及它们中所有都是(除两个Yeast和Mycoplasma Hominis中未评注的ORF以外)已知的(或可能的)蛋白质激酶。再次，丝氨酸/苏氨酸激酶是大多数。

最后，第三个模式"I.H.DLK......D"生成30个SwissProt Rel.34命中，所有它们都是已知的或推定的蛋白质激酶。这显示于图27中。此外，如同三个模式中的第一个的情况，模式"I.H.DLK......D"包括激酶结构域的活态位置。

感兴趣的是注意到所有这三个以上提到的模式都是下列一般模式的特定实例(或一部分)：

[LIVMFYC].[HY].D[LIVMFY]K..N[LIVMFYCT][LIVMFYCT][LIVMFYCT]

其中标志[XYZ]标示一个可由残基X,Y,Z中任何一个占用的位置。此更一般的模式是PROSITE数据库条目，其访问号为PS00108，也即丝氨酸/苏氨酸蛋白质激酶活态位置的标记。注意到此PROSITE标记太专门而无法在由以上看过的三个模式所覆盖的YZ28_METJA区域中检出激酶催化位置。这种情况(在人工智能语言中称为训练组的过度表示)是一个典型的由整个宇宙的有限子集训练的学习系统：始终存在着以下危险：正面例子组(在此情况下是由PROSITE使用的已知的丝氨酸/苏氨酸激酶专用集合)是偏倚的，其结果是学习的特征(此处是激酶标记)在解释观察时不是一般地足够，以致无法有效地外推所考虑家族的新实例(即存在假负面)。为解决此问题，可使用尽可能大的训练组，这也是我们在此处提出的方案的难点。

如上所述，图28提供(a)由YZ28_METJA匹配的所有模式的分布和(b)这三个激酶专用模式所覆盖的区域的图形表示。

在图28(a)中，存在410个由YZ28_METJA匹配的模式(在那些由信息收集阶段所发现的模式中)。如果一个模式在一个残基位置之前(或在该处)起始并在该位置之后(或在该处)结束，则该模式“覆盖”该位置。该图表显示，对于每个残基位置(x轴)有多少模式覆盖该位置。如图28(b)所示，在本文中讨论的三个激酶模式在偏移量35(模式"Y..S..I..DLK")处，112(模式"NIL......IKL")处和1052(模式"I.H.DLK......D")处与序列匹配。此处阐述相对于图28(a)中模式分布的尖峰的偏移量。

使用现有评注

在由YZ28_METJA匹配的410个模式中，只有三个以上分析的模式表现出这种清楚的功能特异性。这并不意味着其余407个没有用。这种从两个序列之间得出一个局部相似性的生物学推理并不始终具有功能特性。有时候同源性标示结构的保留，而另外一些时候它可能对应于所比较序列的总功能中支持作用的功能单元(例如DNA连接结构域)。在发现这类较弱相似性的努力中，我们提供了一种方法，用于发现基础数据库中可用的评注。在下面给出的描述中，我们假定SwissProt评注格式。

SwissProt数据库与其序列区域的大部分序列评注相关连(FT线，见A.Bairoch和R.Apweiler的“SWISS-PROT蛋白质序列数据库及其在1998年内的补充TrEMBL”，Nucleic Acids Res26：38-42,1998)。一个典型区域描述看起来像：

FT DOMAIN 528 779 PROTEINKINASE

其中关键词“FT”标示这是一个区域描述线及其余的线通过给出其起始和结束位置(从相关的数据库序列的残基528直至并包括残基779)及其评注(一个蛋白质激酶结构域)而描述区域。

当存在一个模式P时，我们可以使用(如早已提到的)偏移量表L_D(P)以便在数据库中找到所有与P匹配的序列。假定S是如此的序列，其中在S内有一个与P匹配的子串起始于偏移量j处。如果P竞会落于S的一个评注的区域内(或全部或部分)，我们能够将此区域与P关联起来。为每个与P匹配的序列S完成此过程，结果得到一组与P相关连的区域RS_D(P)。图29给出我们的系统为以上所述三个激酶模式中的一个产生的输出的一部分的例子。也即，图29阐述使用SwissProt评注对个别模式的分析：与模式"I.H.DLK......D"匹配的某些数据库序列。对于每一个这类序列，其ID和DE线被报告(见A.Bairoch和R.Apweiler的“SWISS-PROT蛋白质序列数据库及其在1998年内的补充TrEMBL”，Nucleic Acids Res26：38-42,1998)，给出序列的SwissProt名称及其功能性的简短描述。接着是序列内匹配起始处的偏移量。最后有对于所有具有与由该模式覆盖的区域的交点的评注区域的FT线。

现在给出一个由查询序列Q的子序列A匹配的模式P，问题是如何使用RS_D(P)以便将A特征化。可以使用不少方法。例如，如果RS_D(P)足够大及其大多数成员在它们的功能性方面都一致，则可以推理为A很可能具有相同功能性。另一个考虑的是模式P的相对长度及由FT线描述的区域。例如，如果一个模式P具有15个残基的范围而一个包含P的评注的序列区域具有300种氨基酸的长度，则可能不希望传输该区域的评注至P。总之，希望最终用户应用他/她的专长来决定如何最好地利用由该系统提供的信息。

图30阐述两种使用集合RS_D(P)以便评注YZ28_METJA的区域的方法，从而扩展图28(b)中所画图形。也即，图30显示根据由这些片段匹配的模式的评注来将YZ28_METJA的不同片段特征化。通过利用对于也与这些模式匹配的数据库序列的不同区域而言是可用的信息来评注这些模式。再次相对于整个YZ28_METJA之上的模式分布的尖峰来显示这些片段。第一个方法(图30(b))在下列情况下赋予一个评注X(例如X=跨膜区区域)给模式P：(ⅰ)RS_D(P)的大小至少是15；(ⅱ)RS_D(P)中区域的大部分(80％)评注为X；和(ⅲ)RS_D(P)的每个评注为X的区域的至少50％由P覆盖。第二个方法(图30(c))通过允许由模式覆盖的评注区域的百分比变为30％或更多而分享以上要求(ⅰ)和(ⅱ)而不要(ⅲ)。

性能

一个查询序列的同源性搜索的运行时间决定于：(ⅰ)所用模式П集合的大小；和(ⅱ)Q和数据库序列之间局部相似性的实际数量(由与Q匹配的模式所推理而锝)。对于此处使用的SwissProt Rel.34的情况，在具有256MB存储器的Pentium 266兆赫计算机上搜索其大小大约为一千个残基的查询蛋白质的典型搜索大约需要4-6秒。应该提到，以上所述运行时间是在将所有程序数据(模式和它们的偏移量表)存于存储器中的情况下获得的。对于SwissProt，此数据占用大约200MB。

根据本发明的不同方面，我们提供了一种用于完成序列相似性搜索的方法学，该搜索基于发现蛋白质的基础数据库D上的模式及使用这些模式鉴定查询序列和手头数据库的蛋白质之间的同源性。我们描述了一种使用统计学参量以精确地定义一组模式的方法及讨论了如何将存储器引入统计学计算从而在鉴定重要的同源性时提供更大灵敏度。最后，通过使用SwissProt Rel.34数据库作为实验基础来显示该方法学的用途以及显示如何使用该系统来评注查询序列。在此范围内，我们还讨论了利用所发现的模式以及基础数据库的评注的潜力来将查询与数据库序列之间的甚至弱的相似性加以特征化。

有利的是，本发明的序列同源性检测系统的一个方面与基于现有技术模式的同源性检测工具(例如BLOCKS)明显不同点在于所用模式集合的完整性。从大量训练集合中以无监督方式学习模式以使所有蛋白质都在基础数据库D中。这里并不需要任何将序列”应该“考虑为相同家族成员的偏倚建立先验假定。其结果是，所发现模式更为灵敏。此外，通过一起考虑不同功能性的序列，我们能够发现跨越家族边界的弱相似性(例如描述跨膜区区域的模式)。这类相似性虽然不足以供功能评注的推理用，但能够提供关于被检查的查询序列的不同部分的作用的有用信息。

本发明系统的另一个优点是为完成同源性搜索所需时间。当基因数据库的大小很快地增加时，在每次搜索中使用一些模式而不是扫描整个数据库因而加快速度，这能够成为一个有利因素(尤其是希望进行家中测试而不使用公共服务器的用户)。

虽然本发明的阐述性实施例已经参照附图加以描述，但应该理解，本发明并不就只限于这些这些实施例，熟悉技术的人可在不背离本发明的范围或实质的情况下实现不同其他改变和修改。