测序结果比对方法及其应用

摘要

本发明提出了一种测序结果比对方法。该方法包括：将参考序列按照第一预定长度进行k‑mer切分，构建索引库；基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所述种子序列集合中所述种子序列的数目大于比对所允许的错配数，并且所述种子序列的长度不超过所述第一预定长度；将所述种子序列集合的至少一部分与所述索引库进行匹配；将测序读段与匹配上的参考序列进行全局比对。

权利要求

1.一种测序结果比对方法，其特征在于，包括：将参考序列按照第一预定长度进行k-mer切分，构建索引库；

基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所述种子序列集合中所述种子序列的数目大于比对所允许的错配数，并且所述种子序列的长度不超过所述第一预定长度；

将所述种子序列集合的至少一部分与所述索引库进行匹配；

将测序读段与匹配上的参考序列进行全局比对。

2.根据权利要求1所述的方法，其特征在于，所述全局比对是通过如下方式进行的：获取待比对的来源于所述测序读段的第一核酸序列与来源于所述匹配上的参考序列的第二核酸序列各位置上的基本单元信息；

基于所述基本单元信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列的基本单元数目，n为所述第二核酸序列的基本单元数目，其中所述得分矩阵中的元素Mij表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j个基本单元的比对得分；

基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得经过所述第一核酸序列与所述第二核酸序列的比对结果，

其中，

所述元素Mij是基于下列公式确定的：其中，

Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；

Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；

Mi-1,j表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j个基本单元的比对得分；

g表示小于零的第一预定数值；

S(Ri,Sj)是基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相同时，S(Ri,Sj)为第三预定数值，所述第三预定数值小于所述第二预定数值。

3.根据权利要求1或2所述的方法，其特征在于，所述回溯处理是根据下列步骤确定的；

(a)确定所述矩阵Mmn中的最大值所对应的回溯起始位置；

(b)基于所述回溯起始位置上游相邻三个位置的数值，确定下一回溯位置，其中，所述上游相邻三个位置包括行相邻位置、对角线相邻位置和列相邻位置，其中，选择数值最大的位置作为所述下一回溯位置，并且优先选择所述对角线相邻位置；

(c)重复步骤(b)，直到步骤(b)中所确定的所述下一回溯位置的行号和列号的至少之一为0；

(d)基于步骤(a)-(c)中所确定的回溯路线，确定所述第一核酸序列与所述第二核酸序列的比对结果。

4.根据权利要求2所述的方法，其特征在于，所述第一预定数值为不小于-10的整数，优选-5。

5.根据权利要求2所述的方法，其特征在于，所述第二预定数值为1。

6.根据权利要求2所述的方法，其特征在于，所述第三预定数值为-2。

7.根据权利要求2所述的方法，其特征在于，所述基本单元为碱基。

8.根据权利要求1所述的方法，其特征在于，所述第一预定长度为不超过20的整数，优选不超过15，更优选不超过10，最优选9。

9.根据权利要求1所述的方法，其特征在于，所述种子序列集合是通过下列步骤确定的：

(1)将所述测序读段的5’末端和3’末端各去除第二预定长度的碱基；

(2)将经过步骤(1)处理的所述测序读段划分为多个长度相同的子片段，所述子片段的数目超过所述比对所允许的错配数；

(3)在每个子片段中，由5’末端起始基于所述第一预定长度，确定所述种子序列，其中，如果所述子片段的长度小于所述第一预定长度，则将所述子片段作为所述种子序列，如果所述子片段的长度大于所述第一预定长度，则由所述子片段作的5’末端起始延伸所述第一预定长度作为所述种子序列。

10.根据权利要求9所述的方法，其特征在于，所述第二预定长度为不超过2的整数。

11.根据权利要求9所述的方法，其特征在于，所述子片段的数目＝所述比对所允许的错配数+1。

12.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现如权利要求1-11中任一所述的测序结果比对方法。

13.一种电子设备，其特征在于，包括存储器、处理器；

其中，所述处理器通过读取所述存储器中存储的可执行程序代码来运行与所述可执行程序代码对应的程序，以用于实现如权利要求1-11中任一所述的测序结果比对方法。

14.一种测序结果比对系统，其特征在于，包括：构建索引库装置，所述构建索引库装置用于将参考序列按照第一预定长度进行k-mer切分，构建索引库；

种子序列集合确定装置，所述种子序列集合确定装置用于基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所述种子序列集合中所述种子序列的数目大于比对所允许的错配数，并且所述种子序列的长度不超过所述第一预定长度；

匹配装置，所述匹配装置用于将所述种子序列集合的至少一部分与所述索引库进行匹配；

比对装置，所述比对装置用于将测序读段与匹配上的参考序列进行全局比对。

15.根据权利要求14所述的比对系统，其特征在于，所述比对装置包括：获取基本单元信息单元，所述获取基本单元信息单元用于获取待比对的来源于所述测序读段的第一核酸序列与来源于所述匹配上的参考序列的第二核酸序列各位置上的基本单元信息；

构建得分矩阵单元，所述构建得分矩阵单元与所述获取基本单元信息单元相连，用于基于所述基本单元信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列的基本单元数目，n为所述第二核酸序列的基本单元数目，其中所述得分矩阵中的元素Mij表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j个基本单元的比对得分；

回溯单元，所述回溯单元与所述构建得分矩阵单元相连，用于基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得经过所述第一核酸序列与所述第二核酸序列的比对结果，其中，

所述元素Mij是基于下列公式确定的：其中，

Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；

Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；

Mi-1,j表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j个基本单元的比对得分；

g表示小于零的第一预定数值；

S(Ri,Sj)是基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相同时，S(Ri,Sj)为第三预定数值，所述第三预定数值小于所述第二预定数值。

16.根据权利要求14所述的系统，其特征在于，所述种子序列集合确定装置进一步包括：

末端去除单元，所述末端去除单元用于将所述测序读段的5’末端和3’末端各去除第二预定长度的碱基；

划分子片段单元，所述划分子片段单元用于将经过末端去除单元处理的所述测序读段划分为多个长度相同的子片段，所述子片段的数目超过所述测序的平均错配数；

种子序列确定单元，所述种子序列确定单元用于在每个子片段中，由5’末端起始基于所述第一预定长度，确定所述种子序列，其中，如果所述子片段的长度小于所述第一预定长度，则将所述子片段作为所述种子序列，如果所述子片段的长度大于所述第一预定长度，则由所述子片段作的5’末端起始延伸所述第一预定长度作为所述种子序列。

17.根据权利要求15所述的系统，其特征在于，所述回溯单元进一步包括；

确定回溯起始位置模块，所述确定回溯起始位置模块用于确定所述矩阵Mmn中的最大值所对应的回溯起始位置；

确定下一回溯位置模块，所述确定下一回溯位置模块用于基于所述回溯起始位置上游相邻三个位置的数值，确定下一回溯位置，其中，所述上游相邻三个位置包括行相邻位置、对角线相邻位置和列相邻位置，其中，选择数值最大的位置作为所述下一回溯位置，并且优先选择所述对角线相邻位置；重复确定下一回溯位置，直到所确定的所述下一回溯位置的行号和列号的至少之一为0；

比对结果输出模块，所述比对结果输出模块用于基于所确定的回溯路线，确定所述第一核酸序列与所述第二核酸序列的比对结果。

说明书

测序结果比对方法及其应用

技术领域

[0001] 本发明涉及核酸信息领域，具体地，本发明涉及测序结果比对方法、计算机可读存储介质、电子设备以及测序结果比对系统。

背景技术

[0002] 目前主流的短序列比对工具对于片段比较短(譬如18个碱基)的情况，比对的准确性较低，尤其是这类短片段中存在小的插入、缺失时，基本很难比对上。而对于小RNA分析，
由于部分物种参考序列不一定完全准确，或者测序过程中引入的误差(譬如illumina机器
在高GC区时测序不准确)等原因，导致测序得到的片段存在一些小的突变，却难以准确定位
到相应的参考序列，使得最终的定量结果存在偏差。

[0003] 目前主要应用于小RNA的比对工具，如Bowtie2、BWA等，都是基于种子序列延伸的策略进行比对；对于片段较短且可能存在的突变，导致这些比对工具在进行种子搜索时，容
易检索不到完全匹配的片段，从而导致序列不能被比对上，进而影响了下游的小RNA定量。

[0004] 种子序列定位及延伸算法，通过扫描参考基因组序列，对参考基因组序列建立哈希表，将序列分成一定长度的小片段，这种小片段也被称之为种子。然后，在目标序列中查
找和种子序列相同的片段并标记，以这些标记点为锚点向左右按一定规律延伸比对，将不
合条件的舍弃，符合条件的结果将输出保存。如果种子中存在错配，整条read的比对就不会
被进行，所以就会被认定为非比对上，导致丢失部分比对信息。

[0005] 因此，对于小片段核酸序列的比对方法仍需要科研工作者的进一步开发和改进。

发明内容

[0006] 本发明旨在至少在一定程度上解决相关技术中的技术问题之一。

[0007] 在本发明的第一方面，本发明提出了一种测序结果比对方法。根据本发明的实施例，所述方法包括：将参考序列按照第一预定长度进行k-mer(序列中包含k个碱基的字符
串)切分，构建索引库；基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由
多个种子序列构成，所述种子序列集合中所述种子序列的数目大于比对所允许的错配数，
并且所述种子序列的长度不超过所述第一预定长度；将所述种子序列集合的至少一部分与
所述索引库进行匹配，将测序读段与匹配上的参考序列进行全局比对。根据本发明实施例
的上述测序结果比对方法通过采取比测序的平均错配数多的多个种子序列检索模式，保证
了每条测序读段至少有一条种子序列的匹配，避免了因为种子序列错配而导致匹配不上的
情况；同时，通过种子序列与索引库先匹配，将参考读段与匹配上的参考序列进行全局比
对，而非每条测序读段与所有参考序列进行全局比对，减少了比对的运算量，极大提高了比
对速度；并且通过全局比对，极大地满足了各种类型的小核酸序列比对，如目的片段比参考
序列长、目的片段中存在突变等；提高了小核酸数据的利用率，同时提升了小核酸定量的准
确性。

[0008] 根据本发明的实施例，上述方法还可以进一步包括如下附加技术特征至少之一：

[0009] 根据本发明的实施例，所述全局比对是通过如下方式进行的：获取待比对的来源于所述测序读段的第一核酸序列与来源于所述匹配上的参考序列的第二核酸序列各位置
上的基本单元信息；基于所述基本单元信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列
的基本单元数目，n为所述第二核酸序列的基本单元数目，其中所述得分矩阵中的元素Mij表
示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j个基本单元的比对得分；
基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得经过所述第一核酸序列与所述第二
核酸序列的比对结果，其中，所述元素Mij是基于下列公式确定的：

[0010]

[0011] 其中，Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核
酸序列中第j-1个基本单元的比对得分；Mi-1,j表示所述第一核酸序列中第i-1个基本单元与
所述第二核酸序列中第j个基本单元的比对得分；g表示小于零的第一预定数值；S(Ri,Sj)是
基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的
数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相同时，S(Ri,Sj)为第三
预定数值，所述第三预定数值小于所述第二预定数值。

[0012] 需要说明的是，本申请所述的“比对所允许的错配数”是指在具体比对时，所允许的容错减基数。例如，在小RNA测序中，小RNA比较短，比对所允许的错配数一般为1，如果条
件更严格一些，就是不允许错配；如果所允许的错配数大于1，比对结果可能不太可靠，导致
后续定量等分析结果相对不那么准确。所以在小RNA测序中，比对所允许的错配数一般是1
个错配。

[0013] 根据本发明的实施例，所述第一预定数值为不小于-10的整数，优选-5。

[0014] 根据本发明的实施例，所述第二预定数值为1。

[0015] 根据本发明的实施例，所述第三预定数值为-2。

[0016] 模拟数据测试时发明人发现，所述第一预定数值为-5、所述第二预定数值为1、所述第三预定数值为-2的预定数值组合下，数据的比对效果最佳。

[0017] 根据本发明的实施例，所述基本单元为碱基。

[0018] 根据本发明的实施例，所述第一预定长度为不超过20的整数，优选不超过15，更优选不超过10，最优选9。发明人发现，种子序列的长度过长，可能会因为碱基错配导致整条比
对不上，从而丢失部分比对结果；种子长度过短会增加运算时间，比对速度较慢，但结果更
准确。本方法通过比较不同长度后优选设置长度为9，平衡了比对速度和结果准确性。

[0019] 根据本发明的实施例，所述种子序列集合是通过下列步骤确定的：(1)将所述测序读段的5’末端和3’末端各去除第二预定长度的碱基；(2)将经过步骤(1)处理的所述测序读
段划分为多个长度相同的子片段，所述子片段的数目超过所述测序的平均错配数；(3)在每
个子片段中，由5’末端起始基于所述第一预定长度，确定所述种子序列，其中，如果所述子
片段的长度小于所述第一预定长度，则将所述子片段作为所述种子序列，如果所述子片段
的长度大于所述第一预定长度，则由所述子片段作的5’末端起始延伸所述第一预定长度作
为所述种子序列。

[0020] 根据本发明的实施例，通过数据测试，发明人发现种子序列的长度为9比较合适。发明人发现，种子序列的长度决定了比对的速度和准确性，长度越短比对越准但更慢，长度
越长则快但可能会丢失部分比对信息；

[0021] 根据本发明的实施例，所述第二预定长度为不超过2的整数。

[0022] 根据本发明的实施例，所述子片段的数目＝所述比对所允许的错配数+1。进而保证了一条测序读段至少有一个种子序列能够比对到参考序列。

[0023] 在本发明的第二方面，本发明提出了一种计算机可读存储介质，其上存储有计算机程序。根据本发明的实施例，该程序被处理器执行时实现前面所述的测序结果比对方法。

[0024] 在本发明的第三方面，本发明提出了一种电子设备。根据本发明的实施例，所述电子设备包括存储器、处理器；其中，所述处理器通过读取所述存储器中存储的可执行程序代
码来运行与所述可执行程序代码对应的程序，以用于实现前面所述的测序结果比对方法。

[0025] 在本发明的第四方面，本发明提出了一种测序结果比对系统。根据本发明的实施例，所述系统包括：构建索引库装置，所述构建索引库装置用于将参考序列按照第一预定长
度进行k-mer切分，构建索引库；种子序列集合确定装置，所述种子序列集合确定装置用于
基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所
述种子序列集合中所述种子序列的数目大于比对所允许的错配数，并且所述种子序列的长
度不超过所述第一预定长度；匹配装置，所述匹配装置用于将所述种子序列集合的至少一
部分与所述索引库进行匹配；比对装置，所述比对装置用于将测序读段与匹配上的参考序
列进行全局比对。

[0026] 根据本发明的实施例，所述比对装置包括：

[0027] 获取基本单元信息单元，所述获取基本单元信息单元用于获取待比对的来源于所述测序读段的第一核酸序列与来源于所述匹配上的参考序列的第二核酸序列各位置上的
基本单元信息；

[0028] 构建得分矩阵单元，所述构建得分矩阵单元与所述获取基本单元信息单元相连，用于基于所述基本单元信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列的基本单元数
目，n为所述第二核酸序列的基本单元数目，其中所述得分矩阵中的元素Mij表示所述第一核
酸序列中第i个基本单元与所述第二核酸序列中第j个基本单元的比对得分；

[0029] 回溯单元，所述回溯单元与所述构建得分矩阵单元相连，用于基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得经过所述第一核酸序列与所述第二核酸序列的比对结
果，其中，

[0030] 所述元素Mij是基于下列公式确定的：

[0031]

[0032] 其中，

[0033] Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；

[0034] Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；

[0035] Mi-1,j表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j个基本单元的比对得分；

[0036] g表示小于零的第一预定数值；

[0037] S(Ri,Sj)是基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相
同时，S(Ri,Sj)为第三预定数值，所述第三预定数值小于所述第二预定数值。

[0038] 根据本发明的实施例，所述种子序列集合确定装置进一步包括：末端去除单元，所述末端去除单元用于将所述测序读段的5’末端和3’末端各去除第二预定长度的碱基；划分
子片段单元，所述划分子片段单元用于将经过末端去除单元处理的所述测序读段划分为多
个长度相同的子片段，所述子片段的数目超过所述测序的平均错配数；种子序列确定单元，
所述种子序列确定单元用于在每个子片段中，由5’末端起始基于所述第一预定长度，确定
所述种子序列，其中，如果所述子片段的长度小于所述第一预定长度，则将所述子片段作为
所述种子序列，如果所述子片段的长度大于所述第一预定长度，则由所述子片段作的5’末
端起始延伸所述第一预定长度作为所述种子序列。

[0039] 根据本发明的实施例，所述回溯单元进一步包括；确定回溯起始位置模块，所述确定回溯起始位置模块用于确定所述矩阵Mmn中的最大值所对应的回溯起始位置；确定下一回
溯位置模块，所述确定下一回溯位置模块用于基于所述回溯起始位置上游相邻三个位置的
数值，确定下一回溯位置，其中，所述上游相邻三个位置包括行相邻位置、对角线相邻位置
和列相邻位置，其中，选择数值最大的位置作为所述下一回溯位置，并且优先选择所述对角
线相邻位置；重复确定下一回溯位置，直到所确定的所述下一回溯位置的行号和列号的至
少之一为0；比对结果输出模块，所述比对结果输出模块用于基于所确定的回溯路线，确定
所述第一核酸序列与所述第二核酸序列的比对结果。

[0040] 根据本发明实施例的上述测序结果比对系统适于执行上述测序结果比对方法，其附加技术特征以及附加技术特征所带来的优势效果与前面类似，在此不再赘述。

[0041] 根据本发明实施例的测序结果比对方法和系统以及计算机可读存储介质和电子设备具有以下有益效果：

[0042] 1)通过采取比平均错配数至少多一个的多个种子序列检索模式，保证了种子序列的匹配，避免了因为种子序列错配而导致匹配不上的情况；

[0043] 2)通过全局比对，极大地满足了各种类型的小核酸序列的比对，如目的片段比参考序列长、目的片段中存在突变等；

[0044] 3)提高了小核酸序列数据的利用率，同时提升了其定量的准确性。

附图说明

[0045] 图1是根据本发明实施例的比对装置的结构示意图；

[0046] 图2是根据本发明实施例的回溯单元的结构示意图；

[0047] 图3是根据本发明实施例的测序结果比对系统的结构示意图；

[0048] 图4是根据本发明实施例的种子序列集合确定装置的结构示意图；以及

[0049] 图5是根据本发明实施例的比对结果图。

具体实施方式

[0050] 下面详细描述本发明的实施例，所述实施例的示例在附图中示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

[0051] 一方面，本发明提出了测序结果的比对方法。根据本发明的实施例，所述方法可描述为：

[0052] 1)参考序列按设定种子长度进行序列切分(连续位置截取相应长度的片段)，构建索引库；

[0053] 2)对需要比对的测序读段，记录其5‘端第三个位置到3’端第三个位置的片段(首末端分别允许至多两个滑动)为S，比对时允许的错配数为m；分别从S的m+1等分的起点，截
取种子序列长度的片段，如果到S末端的长度小于种子序列长度时，则截取S末端种子序列
长度的片段，确定种子序列集合；

[0054] 3)将种子序列集合的至少一部分与索引库进行匹配；

[0055] 4)将测序读段与匹配上的参考序列进行Needleman-Wunsh全局比对，具体如下所述：

[0056] a)假定参考序列R的长度为n，测序片段S的长度为m，计算得分矩阵Mmn；

[0057]

[0058] 其中g为indel罚分-5；

[0059]

[0060] 在计算得分矩阵时，记录最大得分值Ss，和相应的位置Sr和Sc；

[0061] Mij是记录参考序列和测序片段每个碱基比对时对应的得分矩阵，其中i和j分别是参考序列和测序片段中第i和j位；Ri指参考序列中第i个碱基，Sj指测序片段中第j个碱基，S
(Ri，Sj)是根据Ri和Sj的碱基的匹配情况给不同的分值。

[0062] 打分值与常规的全局比对不一样(常规一般是匹配1，错配0，indel 0)，因为在核酸体内indel(小的插入、缺失)发生的概率是很低的，错配可以理解为发生了突变，而突变
频率也是很低的，所以在进行测序读段比对时，遇到错配或者indel应该是小概率事件，所
以这样进行罚分。

[0063] b)回溯得分矩阵Mmn，从最大得分值Sr和Sc往前回溯，位置Prc的计算规则如下：

[0064]

[0065] Prc是回溯得分矩阵时，下一个回溯位置的的行和列，其中r和c分别代表对应的行和列；第一个判断条件是指如果对角线的值最大，则往对角线回溯；其他两种情况是分别对
参考序列和测序片段开gap。

[0066] 在回溯过程中，记录累计错配数，如果超过设定阈值，则退出比对；

[0067] 选择最大得分值回溯，可以减少矩阵回溯运算量，同时确保保留最佳比对信息；

[0068] 回溯过程统计累积错配数，可以减少非最佳参考的比对时间，从而提升软件的运行效率。

[0069] c)计算比对得分值，假定匹配数为i，错配数为j，空位数为k，参考序列长度为n，比对得分为：Score＝i-2*j-5*k+i/n；

[0070] 4)按照比对得分进行排序，输出最高得分的比对结果。

[0071] 另一方面，本发明提出了一种核酸序列比对装置。根据本发明的实施例，参考图1，所述比对装置包括：获取基本单元信息单元410，所述获取基本单元信息单元410用于获取
待比对的第一核酸序列与第二核酸序列各位置上的基本单元信息；构建得分矩阵单元420，
所述构建得分矩阵单元420与所述获取基本单元信息单元410相连，用于基于所述基本单元
信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列的基本单元数目，n为所述第二核酸序
列的基本单元数目，其中所述得分矩阵中的元素Mij表示所述第一核酸序列中第i个基本单
元与所述第二核酸序列中第j个基本单元的比对得分；回溯单元430，所述回溯单元430与所
述构建得分矩阵单元420相连，用于基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得
经过所述第一核酸序列与所述第二核酸序列的比对结果，其中，

[0072] 所述元素Mij是基于下列公式确定的：

[0073]

[0074] 其中，

[0075] Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；

[0076] Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；

[0077] Mi-1,j表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j个基本单元的比对得分；

[0078] g表示小于零的第一预定数值；

[0079] S(Ri,Sj)是基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相
同时，S(Ri,Sj)为第三预定数值，所述第三预定数值小于所述第二预定数值。

[0080] 根据本发明的实施例，参考图2，所述回溯单元430进一步包括；确定回溯起始位置模块431，所述确定回溯起始位置模块431用于确定所述矩阵Mmn中的最大值所对应的回溯起
始位置；确定下一回溯位置模块432，所述确定下一回溯位置模块432用于基于所述回溯起
始位置上游相邻三个位置的数值，确定下一回溯位置，其中，所述上游相邻三个位置包括行
相邻位置、对角线相邻位置和列相邻位置，其中，选择数值最大的位置作为所述下一回溯位
置，并且优先选择所述对角线相邻位置；重复确定下一回溯位置，直到所确定的所述下一回
溯位置的行号和列号的至少之一为0；比对结果输出模块433，所述比对结果输出模块433用
于基于所确定的回溯路线，确定所述第一核酸序列与所述第二核酸序列的比对结果。

[0081] 再一方面，本发明提出了一种测序结果比对系统。根据本发明的实施例，参考图3，所述系统包括：构建索引库装置100，所述构建索引库装置100用于将参考序列按照第一预
定长度进行k-mer切分，构建索引库；；种子序列集合确定装置200，所述种子序列集合确定
装置200用于基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子
序列构成，所述种子序列集合中所述种子序列的数目大于测序的平均错配数，并且所述种
子序列的长度不超过所述第一预定长度；匹配装置300，所述匹配装置300用于将所述种子
序列集合的至少一部分与所述索引库进行匹配；比对装置400，所述比对装置400用于将所
述种子序列集合的至少一部分与所述索引库进行比对，其中所述核酸序列比对装置如前面
所限定的。

[0082] 根据本发明的实施例，参考图4，所述种子序列集合确定装置200进一步包括：末端去除单元210，所述末端去除单元210用于将所述测序读段的5’末端和3’末端各去除第二预
定长度的碱基；划分子片段单元220，所述划分子片段单元220用于将经过末端去除单元处
理的所述测序读段划分为多个长度相同的子片段，所述子片段的数目超过所述测序的平均
错配数；种子序列确定单元230，所述种子序列确定单元230用于在每个子片段中，由5’末端
起始基于所述第一预定长度，确定所述种子序列，其中，如果所述子片段的长度小于所述第
一预定长度，则将所述子片段作为所述种子序列，如果所述子片段的长度大于所述第一预
定长度，则由所述子片段作的5’末端起始延伸所述第一预定长度作为所述种子序列。

[0083] 下面详细描述本发明的实施例，所述实施例的示例在附图中示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

[0084] 实施例

[0085] 1)模拟数据

[0086] a)对Human的mature、hairpin和piRNA共28615条进行模拟(其中，Mature、hairpin从miRBase数据库下载，piRNA从piRBase中下载的)，具体规则如下：

[0087] b)从参考序列上随机位置截取长度为18～50的片段，如果片段小于50，则最大为片段长度；输出完全匹配片段；

[0088] c)截取片段随机一个位置错配，输出错配的片段；截取片段的首末端进行随机(1～2)个碱基的Clip片段，输出clip片段；

[0089] d)b)～c)的组合，产生错配和clip的片段；

[0090] e)每一条参考序列，重复1～4步骤随机深度次(同一条参考序列，随机生产模拟数据的次数；例如深度2，第一次执行1-4步骤生产了位置为8～20几个片段，第二次执行1-4步
骤生成了1-19的几个片段，深度的意思就是对同一个参考序列，进行多少次1-4步骤。10以
内)。

[0091] 序列名称为该片段来源以及相应的位置和匹配模式，如hsa_piR_011099_piRNA_1_27M代表来源于hsa_piR_011099_piRNA的第一个位置，匹配模式为27M，其中M代表匹配，X
代表错配，S代表首末端滑动。

[0092] 2)比对结果

[0093] 对目前小RNA比对中常用的比对软件，按照不同的参数对1)中模拟数据进行了比对，具体的软件运行参数见表2：

[0094] 错配数1，种子长度9，线程数1(本方法支持多线程运行)；

[0095] 默认的种子长度是9，发明人之前测试了7～13，发现9的性价比最高(比对率高和运行时间短)，所以默认选的9。

[0096] 表2：软件运行参数

[0097]

[0098]

[0099] 其中，sa为本发明对应的软件(https://github.com/zhuqianhua/sa)。

[0100] 根据各个软件的比对结果，对比其性能，具体结果见图5。

[0101] 其中，perfectly mapped是要求比对的序列和位置都与真实的一致，而part perfectly mapped是相应软件能比对上就算；不管是总体比对率，还是精确位置的比对，本
发明的性能都是优于现有的比对方法的。

[0102] 需要说明的是，这里的perfectly mapped要求比对起始位置和序列名称正确；在真实数据比较时，Star比对上而本方法比对不上的序列，基本上都是首末端有2个以上碱基
滑动,(这个在本申请的方法里面是不认为比对上的，本申请的方法至多允许首末端2个滑
动)。这是因为在真实的小RNA分析时，也不太可能说首末端多出3个甚至更多的碱基，还认
为是同一个小RNA的片段。

[0103] 此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者
隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三
个等，除非另有明确具体的限定。

[0104] 在本发明中，除非另有明确的规定和限定，术语“相连”、“连接”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电
连接或彼此可通讯；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部
的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而
言，可以根据具体情况理解上述术语在本发明中的具体含义。

[0105] 在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特
点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不
必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任
一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技
术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结
合和组合。

[0106] 尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述
实施例进行变化、修改、替换和变型。