CN103793626B

CN103793626B - 碱基序列比对系统及方法

Info

Publication number: CN103793626B
Application number: CN201310367008.7A
Authority: CN
Inventors: 朴旻胥; 朴商贤; 吕润九
Original assignee: IND ACADEMIC COOP; Samsung SDS Co Ltd
Current assignee: IND ACADEMIC COOP; Samsung SDS Co Ltd
Priority date: 2012-10-29
Filing date: 2013-08-21
Publication date: 2017-03-01
Anticipated expiration: 2033-08-21
Also published as: KR101480897B1; US20140121986A1; KR20140056560A; CN103793626A; WO2014069767A1

Abstract

本发明公开一种碱基序列比对系统及方法。根据本发明一个实施例的碱基序列比对系统，用于将包含第一序列及第二序列的一对碱基序列比对到参考序列，包括：种子序列生成单元，从所述第一序列及所述第二序列中分别生成一个以上的片段，并由此构成第一种子序列集合及第二种子序列集合；映射值计算单元，将所述参考序列划分为多个区间，并按所述多个区间分别计算包含于所述第一种子序列集合中的种子序列在对应区间内的第一映射值以及包含于所述第二种子序列集合中的种子序列在对应区间内的第二映射值；比对单元，选择计算出的所述第一映射值及所述第二映射值均为基准值以上的第一区间，并在所述第一区间内搜寻所述第一序列及所述第二序列的映射位置。

Description

碱基序列比对系统及方法

技术领域

本发明的实施例涉及一种用于分析基因组的碱基序列的技术。

背景技术

测序仪用于从原始碱基序列中生成作为长度较短的碱基序列的短片段（read），此时将有一对短片段配对（pair）生成。这样形成配对的短片段是在原始DNA中的预定距离内生成，且根据测序仪的种类的不同，形成为在参考序列中具有反向互补（reversecomplement）方向或相同方向。此时生成的两个短片段之间的距离（insert size）及各短片段的长度为事先根据碱基序列分析目的进行设定，且在相同实验中生成的短片段都具有类似的值。在这些成对的短片段中把先生成的称为5'短片段，后生成的称为3'短片段，而且将5'短片段与3'短片段的方向为反向互补关系的称为双末端短片段（paired-end read），反之5'短片段与3'短片段具有相同方向的称为配对短片段（mate-pair read）。

在对这种双末端短片段或配对短片段进行比对（alignment）时要同时考虑以下三个条件。

（1）各短片段与参考序列之间的碱基序列同源性（homology）

（2）两个短片段的比对方向

（3）两个短片段的比对位置之间的距离

现有技术中的比对算法构成为根据条件（1）将两个短片段分别比对到参考序列上之后，在两个短片段的比对位置中选择满足上述条件（2）、（3）的位置。然而如果这样进行双末端短片段或配对短片段的比对，则为了首先获得满足上述条件（1）的各短片段的比对位置，对参考序列中不满足上述条件（2）、（3）的位置也都要进行搜寻，因此存在不必要的计算量过多的问题。

发明内容

本发明实施例的目的在于提供一种能够确保映射（mapping）准确度的同时通过改善映射过程的复杂度而提高处理速度的针对一对短片段的比对方案。

根据本发明一个实施例的碱基序列比对系统，用于将包含第一序列及第二序列的一对碱基序列比对到参考序列，包括：种子序列生成单元，从所述第一序列及所述第二序列中分别生成一个以上的片段，并由此构成第一种子序列集合及第二种子序列集合；映射值计算单元，将所述参考序列划分为多个区间，并按所述多个区间分别计算包含于所述第一种子序列集合中的种子序列在对应区间内的映射值即第一映射值、以及包含于所述第二种子序列集合中的种子序列在对应区间内的映射值即第二映射值；比对单元，选择计算出的所述第一映射值及所述第二映射值均为基准值以上的第一区间，并在所述第一区间内搜寻所述第一序列及所述第二序列的映射位置。

根据本发明另一实施例的碱基序列比对系统，用于将包含第一序列及第二序列的一对碱基序列比对到参考序列，包括：误差估计单元，分别计算所述第一序列及所述第二序列的最小误差估计值；比对单元，从所述第一序列与所述第二序列中选择计算出的所述最小误差估计值较小的序列，并计算该序列的相对所述参考序列的比对位置，且在以计算出的所述比对位置为基准而设定的可映射范围内执行针对剩余序列的全局比对。

根据本发明一个实施例的碱基序列比对方法，用于在碱基序列比对系统中将包含第一序列及第二序列的一对碱基序列比对到参考序列，包括如下步骤：在种子序列生成单元中，从所述第一序列及所述第二序列中分别生成一个以上的片段，并由此构成第一种子序列集合及第二种子序列集合；在映射值计算单元中，将所述参考序列划分为多个区间，并按所述多个区间分别计算包含于所述第一种子序列集合中的种子序列在对应区间内的映射值即第一映射值、以及包含于所述第二种子序列集合中的种子序列在对应区间内的映射值即第二映射值；在比对单元中，选择计算出的所述第一映射值及所述第二映射值均为基准值以上的第一区间，并在所述第一区间内搜寻所述第一序列及所述第二序列的映射位置。

根据本发明另一实施例的碱基序列比对方法，用于在碱基序列比对系统中将包含第一序列及第二序列的一对碱基序列比对到参考序列，包括如下步骤：在误差估计单元中，分别计算所述第一序列及所述第二序列的最小误差估计值；在比对单元中，从所述第一序列与所述第二序列中选择计算出的所述最小误差估计值较小的序列，并计算该序列的相对所述参考序列的比对位置；以及在所述比对单元中，在以计算出的所述比对位置为基准而设定的可映射范围内执行针对剩余序列的全局比对。

根据本发明的各实施例，在将双末端短片段或配对短片段比对到参考序列时，预先选择具有形成配对的可能性的区间，并在对应区间内执行针对所述双末端短片段或配对短片段的比对，从而与现有方法相比可以显著减少计算量。并且，还可以提供一种在对双末端短片段或配对短片段进行比对时，不仅在特定碱基被置换的情况下可以进行比对，而且在特定碱基被插入或删除而存在缺口（gap）状不一致的情况下也可以进行比对的比对算法。

附图说明

图1为用于说明根据本发明一个实施例的碱基序列比对方法100的图。

图2为用于举例说明根据本发明一个实施例的碱基序列比对方法100的步骤104中的最小误差估计值（MEB）计算过程的图。

图3为用于详细说明根据本发明一个实施例的碱基序列比对方法100中的比对步骤114的顺序图。

图4为用于详细说明根据本发明一个实施例的碱基序列比对方法100中的有效配对搜寻过程的顺序图。

图5为示出根据本发明一个实施例的碱基序列比对系统500的模块图。

图6为示出根据本发明另一实施例的碱基序列比对系统600的模块图。

符号说明：

500、600：碱基序列比对系统 502：种子序列生成单元

504：映射值计算单元 506：比对单元

602：误差估计单元 604：比对单元

具体实施方式

以下参照附图说明本发明的具体实施方式。然而这仅仅是示例，本发明并不局限于此。

在对本发明进行说明时，如果遇到对有关本发明的公知技术的具体说明有可能不必要地干扰本发明的主旨的情况，则省略其详细说明。并且，后述的术语均为考虑本发明中的功能而进行定义的，其可能因使用者、运用人员的意图或习惯等而有所不同。因此，要以整个说明书的内容为基础对其进行定义。

本发明的技术思想由权利要求书确定，以下的实施例只是为了将本发明的技术思想有效地传递给本发明所属技术领域中具有普通知识的人员而采用的一种手段。

在对本发明的实施例进行具体说明之前，首先对本发明中使用的术语进行如下说明。

首先，“短片段序列（read sequence）”（或者简称为“短片段（read）”）是指基因组测序仪（genome sequencer）中输出的短碱基序列数据。短片段的长度因基因组测序仪的种类而不同，通常构成为35～500bp(base pair)范围的多种长度，在DNA碱基的情况下，通常用A、C、G、T等四个字母表示。

在本发明的实施例中，基因组测序仪输出一对配对（pair）的短片段。此时，将所述一对短片段中的第一个短片段称为5’短片段而将第二个短片段称为3’短片段，所述5’短片段与3’短片段的方向可形成为反向互补（reverse complement）关系（双末端短片段），或者形成为相同的方向（配对短片段）。例如，对于双末端短片段而言，如果5’短片段为正向（forward）短片段，则3’短片段将是反向互补（reverse complement）短片段，与之相反，如果5’短片段为反向互补短片段，则3’短片段将是正向短片段。并且，对于配对短片段而言，如果5’短片段为正向短片段，则3’短片段也将是正向短片段，与之相反，如果5’短片段为反向互补方向的短片段，则3’短片段也将是反向互补方向的短片段。

“参考序列（reference sequence）”指可对利用所述短片段形成整个碱基序列提供参考的碱基序列。在碱基序列分析中，通过将基因组测序仪所输出的大量短片段参照参考序列进行映射而完成整个碱基序列。在本发明中，所述参考序列既可以是碱基序列分析时预先设定的序列（例如人类的整个碱基序列等），或者也可以将基因组测序仪中产生的碱基序列使用为参考序列。

“碱基（base）”为构成参考序列及短片段的最小单位。如前所述，构成DNA的碱基可由A、C、G、T等四个字母表示的碱基构成，将这些分别称为碱基。换言之，对于DNA而言，可用四种碱基表示，短片段也是如此。

“片段序列（fragment sequence）”（或者简称为“片段”）指一种序列，该序列成为为了短片段的映射而比较短片段与参考序列时的单位。从理论上讲，为了将短片段映射于参考序列，需要把整个短片段从参考序列的最前端部分开始依次比较并计算短片段的映射位置。然而，由于这种方法在映射一个短片段时消耗过多的时间并要求过高的计算能力，因此实际上要先把短片段的一部分所构成的片，即片段映射于参考序列而找出整个短片段的映射候选位置，然后将整个短片段映射于对应候选位置（Global Alignment）。

“种子序列（seed sequence）”指由短片段产生的片段中与参考序列相匹配的片段。即，在本发明的实施例中将由短片段产生的各片段分别与参考序列进行精确匹配（exact matching），并进行用于将其中不与所述参考序列精确匹配的片段排除在外的筛选过程，且将所述筛选过程中精确匹配的片段作为种子序列，而将这些种子序列的集合称为种子序列集合。此时，与所述参考序列相匹配的片段指与所述参考序列进行精确匹配（exact matching）时不一致的碱基数为预先设定的允许值以下的片段。此时，如果所述允许值为0，则种子序列集合中只包含与所述参考序列精确匹配（即，没有不一致的碱基）的片段。

图1为用于说明根据本发明一个实施例的碱基序列比对方法100的图。在本发明的实施例中，碱基序列比对方法100指通过将基因组测序仪（genome sequencer）中输出的一对短片段（双末端短片段或配对短片段）与参考序列进行比较而确定对应短片段在所述参考序列中的映射（或比对）位置的一系列过程。在以下实施例中，将构成所述一对短片段的两个短片段（5’短片段及3’短片段）分别称为第一短片段及第二短片段。

首先，当从基因组测序仪（genome sequencer）接收到第一短片段及第二短片段（步骤102）时，分别针对输入的两个短片段的正向序列及反向互补序列而计算最小误差估计值（MEB;Minimum Error Bound）（步骤104）。即，在本步骤中将分别计算包含第一短片段的正向序列、第一短片段的反向互补序列、第二短片段的正向序列、第二短片段的反向互补序列在内的四个序列的最小误差估计值。此时，所述最小误差估计值是指将所述各序列映射于参考序列时可能发生的误差的最小值。

图2为用于举例表示所述步骤104中的MEB计算过程的图。首先，如图2的（a）所示，将初始MEB设定为0，并从对象序列的第一个碱基开始向右侧逐个移动而尝试精确匹配。此时，如图2（b）所示，假定在对象序列的特定碱基（在图中以第二个T标记的地方）开始无法再实现精确匹配，则这种情况说明从序列的匹配起始位置到当前位置之间的区间中的某处出现了误差。因此，在这种情况下将MEB值增加1（MEB=1）之后在下一个位置上重新开始精确匹配（在图中标记为(c)）。如果在以后又遇到判断为无法精确匹配的情况，则是说明从重新开始精确匹配的位置到当前位置之间的区间某处又出现了误差，因此又将MEB值增加1（MEB=2）之后在下一个位置上重新开始精确匹配（在图中标记为(d)）。通过这种过程，到达序列末端时的MEB值将成为对应序列的MEB值。

通过如上所述的过程，将分别计算出包括第一短片段的正向序列、第一短片段的反向互补序列、第二短片段的正向序列、第二短片段的反向互补序列在内的共计4个序列各自的MEB值。

然后将计算出的4个MEB值与预先设定的最大误差允许值（maxError）进行比较（步骤106）。此时，如果计算出的4个MEB值均超过所述最大误差允许值，则判定针对对应短片段的比对失败。

与之相反，如果在所述步骤106中判断的结果至少有部分序列的MEB为所述最大误差允许值以下，则选择计算出的MEB为最大误差允许值以下的序列（步骤108），并构造出所选择序列各自的种子序列集合（步骤110）。然后将所述参考序列划分为多个区间，并按所述多个区间分别计算所述所选择序列的总映射值而生成映射直方图（步骤112），且利用所述映射直方图而将所述一对短片段比对到所述参考序列（步骤114）。

以下详细说明所述步骤110至步骤114的具体过程。

由选择的序列构成种子序列集合（步骤110）

该步骤为利用从所述步骤108中选择的短片段序列生成一个以上的种子序列的步骤。首先，考虑所选择序列的一部分或全部而生成多个片段。例如，可通过将所述序列的全部或特定区间分割为多个片或者将分割的片进行组合而生成片段。这种情况下，生成的片段可以连续地相连，然而并非一定要如此，也可以用序列内分离的片的组合构成片段。并且，生成的片段并非一定要具有相同的长度，在一个短片段内也可以生成多种长度的片段。总而言之，本发明中的由短片段序列生成片段的方法并不受到特别的局限，从短片段序列的一部分或全体中提取片段的各种算法均可不受限制地使用。

如果通过上述过程生成了与所选择的各序列分别对应的片段，接着便通过从生成的片段中除去与参考序列不匹配的片段的筛选过程而构成种子序列集合。即，尝试所生成的片段与所述参考序列之间的精确匹配（exact matching），然后用不一致的碱基数为预先设定的允许值以下的片段（种子序列）构成种子序列集合。此时，可通过适当考虑序列长度及从中提取的片段长度等而确定所述允许值。例如，序列长度较短的情况下（约为50bp以下），优选为只对与所述参考序列精确匹配的片段予以考虑，此时所述允许值可以为0。而随着序列长度变长，可将所述允许值提高为1或2等，从而防止映射的准确率过低。

生成映射直方图（步骤112）

如果通过上述过程构成了种子序列集合，接着便构成对应于各序列的映射直方图（histogram）。在本发明中，映射直方图为具有预定大小的整数阵列（integer array），整数阵列的值对应于将参考序列划分为具有相同大小的多个区间时的各区间。例如，将参考序列划分为具有65536(=2¹⁶)bp大小的区间时，参考序列的0～65535bp的区间对应于映射直方图h的第一个值h[0]，而65536～131071的区间对应于映射直方图h的第二个值h[1]。可通过这种方式使参考序列的各划分区间对应于映射直方图。

而且，映射直方图的各值h[i]中存储有在对应的参考序列区间当中按各短片段序列分别提取的种子序列的总映射值。此时，所述映射值可以是对应参考序列区间中的所述种子序列的总映射长度。例如，假设从特定短片段序列中提取的种子序列中53-67种子序列（从所述短片段序列的第53～第67个碱基提取的种子序列）及61-75种子序列被映射于映射直方图的第一个区间，则在这种情况下对应区间的直方图值将是23(=75-53+1)。

另外，所述映射值也可以是对应参考序列区间中的所述种子序列的总映射个数。在上面的例中，由于映射于映射直方图的第一个区间的种子序列个数为2，因此对应区间的直方图值将是2。并且，根据实施例的不同，也可以将多个区间各自的总映射长度及总映射个数作为所述映射值一并存储。

一对短片段的比对（步骤114）

如果通过上述过程生成了第一短片段及第二短片段的序列各自的映射直方图，则利用所生成的映射直方图而将所述一对短片段比对到所述参考序列。

图3为用于详细说明根据本发明一个实施例的比对步骤（步骤114）的顺序图。

首先，判断是否可以利用所述步骤106中选择的短片段序列构成序列对（Sequencepair）（步骤300）。

例如，在所述一对短片段为双末端短片段的情况下，判断MEB值为基准值（即最大误差允许值）以下的序列是否可以构成如下配对中的至少一个。

（第一短片段的正向序列-第二短片段的反向互补序列）

（第一短片段的反向互补序列-第二短片段的正向序列）

而如果所述一对短片段为配对短片段，则判断MEB值为基准值（即最大误差允许值）以下的序列是否可以构成如下配对中的至少一个。

（第一短片段的正向序列-第二短片段的正向序列）

（第一短片段的反向互补序列-第二短片段的反向互补序列）

如果在所述步骤300中判断的结果为至少可以在上述配对中实现一种构成，则对构成序列对的两个短片段序列的直方图值进行比较，从而判断是否存在两个序列的直方图值均为直方图切值（Histogram Cut）以上的参考序列区间（步骤302）。

如果在所述步骤302中判断的结果为存在两个序列的直方图值（映射值）均为直方图切值（Histogram Cut，H）以上的参考序列区间的情况下，将对应区间选择为映射对象区间（步骤304），并在所选区间内对构成所述序列对的两个短片段序列执行一次比对（步骤306、步骤308）。具体而言，在所述步骤306中，在所述映射对象区间内执行分别针对构成序列对的两个短片段序列的全局比对（global alignment），并将根据所述全局比对结果计算出的两个短片段序列的比对位置对当中满足预先设定的短片段间距范围（插入大小，insert size）的比对位置对（有效配对，valid pair）选择为所述第一短片段及所述第二短片段的比对位置。此时，所述有效配对要满足如下三个条件。

（1）两个序列的比对方向要与初始输入的一对短片段相同或对应。当输入的一对短片段为双末端短片段时各序列应具有反向互补关系。即，如果一个序列为正向序列，则另一个序列应该是反向互补序列。而且，当输入的一对短片段为配对短片段时，两个序列的比对方向应该相同。

（2）两个序列中的至少一个应具有最大误差允许值以下的误差。

（3）两个序列的比对位置间距应该在预先设定的可映射范围之内。此时，所述可映射范围可用如下数学式1确定。

[数学式1]

L₁-k·D≤L₂≤L₁+k·D

（L₁为构成序列对的第一个序列的映射位置，L₂为第二个序列的映射位置，作为加权值的k具有大于0而小于1.8的值，D为预先设定的序列间的距离差(插入大小)）

此时，之所以给所述插入大小赋予加权值k，是因为碱基序列的特性决定了一些碱基的插入或删除可能导致序列间的距离发生变化，故使用加权值k来反映。

如果举例说明搜寻所述有效配对的过程则如图4所示。假定在图示的映射对象区间内构成序列对的两个序列中的第一序列映射于位置A和B，而第二序列映射于位置C。这种情况下将生成如下的两个比对位置对。

(A,C)

(B,C)

假定所述A与C之间的插入大小d₁为1500bp，B与C之间的插入大小d₂为650bp，基于所述数学式1的可映射范围为-750bp～750bp。这种情况下，由于在两个比对位置对当中满足前述可映射范围的是(B，C)，因此所述第一短片段和第二短片段的比对位置将是B和C。

如上所述，将在所选区间内满足前述范围的比对位置对称为有效配对（validpair）。即，在上述例中有效配对是(B，C)，如果找到该有效配对，则对应双末端短片段的比对即成功。

而如果与之相反，在所述步骤304中在所选区间内进行一次比对的结果不存在有效配对，或者在所述步骤302中判断的结果为不存在两个序列的直方图值均为H以上的区间，则将构成序列对的两个序列中的某一个序列的直方图值为H以上的区间选择为映射对象区间（步骤310），并在所选映射对象区间当中执行二次比对（步骤312、步骤314）。

对所述二次比对过程进行如下更为详细的说明。首先，在两个序列中选择一个序列，并在所选序列的所述映射区间内计算比对位置。此时，所选的序列可以是两个序列当中在对应映射对象区间内的直方图值为H以上的序列。

然后，判断剩余序列在以计算出的所述比对位置为基准而设定的可映射范围内是否得到映射（局部比对，local alignment）。即，判断所述可映射范围内是否存在满足前述三个条件的有效配对。此时，所述可映射范围如同前述数学式1所示。即，在该二次比对过程中是将直方图值较大的序列作为一种锚点（anchor）来判断剩余序列在对应序列的周边是否得到映射。

如果进行所述映射的结果存在有效配对，则对应的一对短片段的比对成功完成。而如果与之不同，即执行所述步骤312、步骤314的结果发现不存在有效配对，则所述短片段的比对失败，在这种情况下将所述第一短片段及第二短片段分别在参考序列中进行全局比对，并输出进行所述全局比对的结果当中比对分数（alignment score）最高的比对位置（步骤322）。此时，由于与各短片段的全局比对及比对分数的计算有关的内容为本发明所属技术领域中的普通知识，因此省略其详细说明。

另外，如果在所述步骤300中判断的结果为无法构成两个序列的MEB均为最大误差允许值以下的序列对，接着便判断两个当中哪一序列的MEB是否为最大误差允许值以下（步骤316）。此时，当在所述步骤316中判断的结果某一序列的MEB为最大误差允许值以下的情况下，计算MEB为最大误差允许值以下的序列的相对所述参考序列的比对位置（步骤318，single end alignment）。然后，判断在以计算出的所述比对位置为基准设定的可映射范围内，是否存在剩余序列满足前述三个条件的有效配对（步骤320，local alignment）。此时，所述可映射范围如同前述数学式1所示。即，在该二次比对过程中是将MEB为最大误差允许值以下的序列作为一种锚点（anchor）来判断剩余序列在对应的序列周围是否得到映射。

如果进行所述映射的结果存在有效配对，则对应的一对短片段的比对成功完成。而如果与之不同，即执行所述步骤318、320的结果不存在有效配对，则所述一对短片段的比对失败，在这种情况下将所述第一短片段及第二短片段分别在参考序列中进行全局比对，并输出进行所述全局比对的结果当中比对分数（alignment score）最高的比对位置（步骤322）。而且，在所述步骤316中判断的结果所有序列的MEB值均超过最大误差允许值的情形也进行与此相同的处理。

计算直方图切值（Histogram Cut）

在上述实施例中，可通过如下方式计算直方图切值。

首先，以映射于对应区间的种子序列的个数定义所述直方图值（即，各区间内的映射值）的情况下，所述直方图切值至少应该是2。这是由于映射的基本单位为种子序列，因此只映射一个种子序列的区间发生短片段映射的可能性很低。即，以映射于各区间的种子序列的个数定义所述直方图值的情况下，可通过恰当考虑短片段长度、种子序列长度等而从2以上的整数值中确定所述直方图切值。

其次，以映射于对应区间的种子序列的长度定义所述直方图值的情况下，以如下方式计算直方图切值。在f表示片段大小、s表示为了生成片段而在短片段中移动的距离、L表示短片段长度、e表示短片段中允许的最大误差个数、H表示直方图切值的情况下，可用如下数学式求出短片段中不受误差影响的区域长度T。

T=L–f·e-s

此时，由于L和e为执行本发明时预先确定的值，因此由f和s的值决定T。即，算法的性能变化取决于f和s的值如何变化。

首先，在确定H的值时考虑以下两个条件。其中，必须条件为必须要满足的条件，而附加条件只在可能的情况下予以考虑。

必须条件：由于映射的基本单位为片段，因此无论直方图切值多小，至少要具有能够包含重叠（overlap）的两个以上片段的大小。例如图2所示，在f=15、s=4的情况下，由于重叠的两个片段的最小长度为15+4=19，因此H值应该至少为19。而且，由于要将所述H值设定为至少包含两个片段，因此H值至少要比f+s更大或相等。如后所述，f值应至少为15，因此将s值假定为其最小值1的情况下，H值至少为16(=15+1)。

附加条件：在理想情况下，通过设定H=T并搜寻映射了T以上的序列的直方图，便可以找到对应于给定误差的所有映射。然而，如前所述，在参考序列本身包含许多重复的情况下，根据情况可能会遇到需要增大片段长度的情形。因此，考虑到这一点，在确定H值时使用比T略小的T–s可能有利于映射率。如果假定H=T，则H=L-f·e-s，如果假定其中的e取最小值1（由于e=0的情况为与参考序列精确匹配的情形，因此此时将在前述步骤104中映射完毕），则有H=L-f-s。该值将是直方图值的最大值。如果假定L=75bp、f=15bp、s=1，H的最大值便成为75-15-1=59。

综上，所述H值应该满足如下范围。

f+s≤H≤L–(f+s)

然后，在满足以下两个条件的值当中选择较大值作为f值。必须条件仍然要必须满足，而附加条件只在可能的情况下考虑。

必须条件：f应该取15以上，这是由于如果片段长度为14以下，则参考序列中的映射位置的个数将急剧增加。

如下的表1表示根据片段长度的在人类基因组中的片段平均出现频率。

[表1]

片段长度	平均出现频率
		10	2726.1919
11	681.9731
		12	170.9185
13	42.7099
		14	10.6470
15	2.6617
		16	0.6654
17	0.1664

从以上的表中可知，片段长度为14以下的情况下各片段的频率为10以上，而在片段长度为15的情况下出现频率减少为3以下。即，相比于将片段长度设定为14以下的情况而言，将片段长度设定为15以上的情况能够大幅度减少片段重复。

附加条件：为了将T的长度确保为两个片段的大小以上，要满足f≤L/(e+2)。

例如，L=100、e=4的情况下，f要具有16以下的值。

综合以上条件，确定f、s、H的方法可整理如下。

——将s固定为4之后确定f和H。

——将15≤f≤L/(e+2)范围内的最大值确定为f（但必须满足f≥15）。

——H为通过如下数学式确定。

通过H=L–f·e–2s或H=f+s计算的值中的较大值（其中，H为基准值，L为短片段长度，f为片段长度，e为短片段的最大误差个数，s为各片段的移动间距）。

例1：当L=75、e=3时，

由于f=15～15，故f=15，

s=4，

H=75–3×15–2×4=22。

例2：当L=100、e=4时，

由于f=15～16，故f=16，

S=4，

H=100–4×16–2×4=36–8=28。

例3：当L=75、e=4时，

虽然f=15～12，然而由于f应大于等于15，故f=15，

s=4，

虽然H=75–4×15–2×4=15-8=7，然而由于f+s=19，故结果将是H=19。

图5为根据本发明一个实施例的碱基序列比对系统500的模块图。根据本发明一个实施例的碱基序列比对系统500为用于将具有同向关系或反向互补关系的第一序列及第二序列比对到参考序列的系统，包括种子序列生成单元502、映射值计算单元504、比对单元506。

种子序列生成单元502从所述第一序列及所述第二序列中分别生成一个以上的片段（fragment），并由此构成第一种子序列集合及第二种子序列集合。所述第一种子序列集合中只包括从所述第一序列提取的一个以上片段（fragment）中与所述参考序列匹配的片段，而所述第二种子序列集合只包括从所述第二序列提取的一个以上片段中与所述参考序列匹配的片段。并且，所述的与参考序列匹配的片段是指进行与所述参考序列之间的精确匹配（exact matching）的结果不一致的碱基数为设定数目以下的片段。

映射值计算单元504将所述参考序列划分为多个区间，并对多个所述区间分别计算第一映射值和第二映射值。此时，所述第一映射值可以是包含于所述第一种子序列集合中的种子序列在对应区间内的总映射长度，而所述第二映射值可以是包含于所述第二种子序列集合中的种子序列在对应区间内的总映射长度。并且，也可以将所述第一映射值定义为包含于所述第一种子序列集合中的种子序列在对应区间内的总映射个数，而所述第二映射值可以定义为包含于所述第二种子序列集合中的种子序列在对应区间内的总映射个数。

比对单元506选择计算出的所述第一映射值及所述第二映射值均为基准值以上的第一区间，并在所述第一区间内搜寻所述第一序列及所述第二序列的映射位置。具体而言，比对单元506在所述第一区间内执行针对所述第一序列及所述第二序列的全局比对（global alignment），并将进行所述全局比对的结果计算出的所述第一序列及所述第二序列的比对位置对当中满足预先设定的序列之间的距离范围的比对位置对选择为所述第一序列及所述第二序列的比对位置。

如果不存在所述第一映射值及所述第二映射值均为基准值以上的区间，比对单元506便在所述第一映射值与所述第二映射值中的某一个值为基准值以上的第二区间内搜寻所述第一序列及所述第二序列的映射位置。具体而言，比对单元506在所述第二区间内计算针对从所述第一序列与所述第二序列中选择的序列的比对位置，并在以计算出的所述比对位置为基准而设定的可映射范围内执行针对剩余序列的全局比对。

此时，选择的所述序列可以是所述第一序列和所述第二序列当中在所述第二区间内的值更大的序列。另外，所述可映射范围可以是以选择的所述序列的映射位置为基准向所述参考序列的前后端各有k×D（此时，k为加权值，D为预先设定的序列之间的距离）长度的区间，这种情况下所述加权值k可以为1.8以下。

图6为根据本发明另一实施例的碱基序列比对系统600的模块图。根据该实施例的碱基序列比对系统600为用于将具有同向关系或反向互补关系的第一序列及第二序列比对到参考序列的系统，包括误差估计单元602及比对单元604。

误差估计单元602分别计算所述第一序列及第二序列的最小误差估计值。具体而言，误差估计单元602由从所述第一序列与所述第二序列中选择的序列的第一个碱基开始每次移动一个碱基而将所述选择的序列与所述参考序列进行精确匹配，而在所选序列的特定位置上无法实现精确匹配的情况下，从对应位置的下一个碱基开始每次移动一个碱基并重新执行精确匹配，并将到达所述选择的序列的末尾碱基为止判断为无法实现精确匹配的位置的个数设定为所选序列的最小误差估计值。由于有关误差估计单元602中的最小误差估计值计算的内容已在图2及相关说明中充分说明，因此在此处省略重复说明。

比对单元604从所述第一序列与所述第二序列中选择计算出的所述最小误差估计值较小的序列，并计算该序列相对所述参考序列的比对位置，且在以计算出的所述比对位置为基准而设定的可映射范围内进行针对剩余序列的全局比对。此时，所述可映射范围可以是以所选序列的映射位置为基准向所述参考序列的前后端各有k×D（此时，k为加权值，D为预先设定的序列之间的距离）长度的区间，这种情况下所述加权值k可以为1.8以下。

另外，本发明的实施例可以包括记录有用于将本说明书中记载的方法在计算机上执行的程序的计算机可读记录介质。所述计算机可读记录介质可将程序命令、本地数据文件、本地数据结构等单独或组合而包括在内。所述介质既可以是为了本发明而特别设计并构成的，也可以是计算机软件领域中具有普通知识的人员所公知而能够使用的。计算机可读记录介质的实例中包括硬盘、软盘、磁带等磁介质；只读光盘（CD-ROM）、DVD等光记录介质；软盘等磁光介质；只读存储器、随机存储器、闪存等为了存储并执行程序命令而特意构成的硬件装置。程序命令的实例中不仅包括通过编译器（Compiler）制作的机器语言代码，而且还可以包括借助于解释器（Interpreter）等而能够在计算机上执行的高级语言代码。

以上通过代表性的实施例对本发明进行了详细说明，然而本发明所属技术领域中具有普通知识的人员即可明白在不脱离本发明范围的条件下对上述实施例能够进行多种多样的变形。

因此不能局限于上述实施例而确定本发明的权利范围，本发明的范围应当由权利要求书及其等价内容确定。

Claims

1.一种碱基序列比对系统，用于将包含第一序列及第二序列的一对碱基序列比对到参考序列，包括：

误差估计单元，针对所述第一序列及所述第二序列各自的正向序列及反向互补序列，从第一个碱基开始到末尾碱基为止以一个碱基为单位逐个移动而与所述参考序列进行精确匹配，并将判断为无法实现精确匹配的位置的个数设定为所述第一序列及所述第二序列各自的正向序列及反向互补序列的最小误差估计值；

种子序列生成单元，选择所述第一序列及所述第二序列各自的正向序列及反向互补序列中所述最小误差估计值为预先设定的最大误差允许值以下的序列，并从所述选择的各个序列中分别生成一个以上的片段，由此构成所述选择的各个序列的种子序列集合；

映射值计算单元，将所述参考序列划分为多个区间，并按所述多个区间分别计算包含于所述选择的各个序列的种子序列集合中的种子序列在对应区间内的映射值；

比对单元，从所述多个区间中选择包含于所述选择的各个序列的种子序列集合中的种子序列的映射值均为基准值以上的第一区间，并在所述第一区间内搜寻所述选择的各个序列的映射位置。

2.如权利要求1所述的碱基序列比对系统，其特征在于，所述种子序列集合只包括从所述选择的序列提取的一个以上的片段中与所述参考序列相匹配的片段。

3.如权利要求2所述的碱基序列比对系统，其特征在于，与所述参考序列相匹配的片段为进行与所述参考序列之间的精确匹配的结果，不一致的碱基数为设定个数以下的片段。

4.如权利要求1所述的碱基序列比对系统，其特征在于，所述映射值计算单元基于包含在所述种子序列集合中的种子序列在对应区间内的总映射长度而计算所述映射值。

5.如权利要求1所述的碱基序列比对系统，其特征在于，所述映射值计算单元基于包含在所述种子序列集合中的种子序列在对应区间内的总映射个数而计算所述映射值。

6.如权利要求1所述的碱基序列比对系统，其特征在于，所述比对单元在所述第一区间内执行针对所述选择的各个序列的全局比对，并将进行所述全局比对的结果计算出的比对位置对当中满足预先设定的序列之间的距离范围的比对位置对选择为所述选择的序列的比对位置。

7.如权利要求1所述的碱基序列比对系统，其特征在于，当无法选择所述第一区间时，所述比对单元便选择包含在所述选择的各个序列的种子序列集合中的种子序列的映射值中的某一映射值为基准值以上的第二区间，并在选择的所述第二区间内搜寻所述选择的各个序列的映射位置。

8.如权利要求7所述的碱基序列比对系统，其特征在于，所述比对单元在所述第二区间内计算相对所述选择的序列中的一个序列的比对位置，并在以计算出的所述比对位置为基准而设定的可映射范围内执行针对剩余序列的全局比对。

9.如权利要求8所述的碱基序列比对系统，其特征在于，所述选择的序列为所述选择的序列当中在所述第二区间内的映射值更大的序列。

10.如权利要求8所述的碱基序列比对系统，其特征在于，所述可映射范围是以所述计算的比对位置为基准向所述参考序列的前后端各延伸k×D长度的区间，其中，k为加权值，D为预先设定的序列之间的距离。

11.如权利要求10所述的碱基序列比对系统，其特征在于，所述加权值k为1.8以下。

12.一种碱基序列比对系统，用于将包含第一序列及第二序列的一对碱基序列比对到参考序列，包括：

误差估计单元，分别计算所述第一序列及所述第二序列的最小误差估计值；

比对单元，从所述第一序列与所述第二序列中选择计算出的所述最小误差估计值较小的序列，并计算该序列的相对所述参考序列的比对位置，且在以计算出的所述比对位置为基准而设定的可映射范围内执行针对剩余序列的全局比对，

其中，所述误差估计单元由从所述第一序列与所述第二序列中选择的序列的第一个碱基开始以一个碱基为单位逐个移动而将所述选择的序列与所述参考序列进行精确匹配，而在所述选择的序列的特定位置上无法实现精确匹配的情况下，从对应位置的下一个碱基开始以一个碱基为单位逐个移动的同时重新执行精确匹配，并在到达所述选择的序列的末尾碱基时，将判断为无法实现精确匹配的位置的个数设定为所述选择的序列的最小误差估计值。

13.一种碱基序列比对方法，用于在碱基序列比对系统中将包含第一序列及第二序列的一对碱基序列比对到参考序列，包括如下步骤：

在误差估计单元中，针对所述第一序列及所述第二序列各自的正向序列及反向互补序列，从第一个碱基开始到末尾碱基为止以一个碱基为单位逐个移动而与所述参考序列进行精确匹配，并将判断为无法实现精确匹配的位置的个数设定为所述第一序列及所述第二序列各自的正向序列及反向互补序列的最小误差估计值；

在种子序列生成单元中，选择所述第一序列及所述第二序列各自的正向序列及反向互补序列中所述最小误差估计值为预先设定的最大误差允许值以下的序列，并从所述选择的各个序列中分别生成一个以上的片段，由此构成所述选择的各个序列的种子序列集合；

在映射值计算单元中，将所述参考序列划分为多个区间，并按所述多个区间分别计算包含于所述选择的各个序列的种子序列集合中的种子序列在对应区间内的映射值；

在比对单元中，从所述多个区间中选择包含于所述选择的各个序列的种子序列集合中的种子序列的映射值均为基准值以上的第一区间，并在所述第一区间内搜寻所述选择的各个序列的映射位置。

14.如权利要求13所述的碱基序列比对方法，其特征在于，所述种子序列集合只包括从所述选择的序列提取的一个以上的片段中与所述参考序列相匹配的片段。

15.如权利要求14所述的碱基序列比对方法，其特征在于，与所述参考序列相匹配的片段为进行与所述参考序列之间的精确匹配的结果，不一致的碱基数为设定个数以下的片段。

16.如权利要求13所述的碱基序列比对方法，其特征在于，在进行所述计算的步骤中，基于包含在所述种子序列集合中的种子序列在对应区间内的总映射长度而计算所述映射值。

17.如权利要求13所述的碱基序列比对方法，其特征在于，在进行所述计算的步骤中，基于包含在所述种子序列集合中的种子序列在对应区间内的总映射个数而计算所述映射值。

18.如权利要求13所述的碱基序列比对方法，其特征在于，搜寻所述映射位置的过程包括如下步骤：

在所述第一区间内执行针对所述选择的各个序列的全局比对；

将进行所述全局比对的结果计算出的比对位置对当中满足预先设定的序列之间的距离范围的比对位置对选择为所述选择的序列的比对位置。

19.如权利要求13所述的碱基序列比对方法，其特征在于，搜寻所述映射位置的过程还包括如下步骤：

当无法选择所述第一区间时，选择包含在所述选择的各个序列的种子序列集合中的种子序列的映射值中的某一映射值为基准值以上的第二区间，并在选择的所述第二区间内搜寻所述选择的各个序列的映射位置。

20.如权利要求19所述的碱基序列比对方法，其特征在于，在搜寻所述映射位置的过程中，在所述第二区间内计算相对所述选择的序列中的一个序列的比对位置，并在以计算出的所述比对位置为基准而设定的可映射范围内执行针对剩余序列的全局比对，且所述选择的序列为所述选择的序列当中在所述第二区间内的映射值更大的序列。