CN115641911B - 一种针对序列间重叠检测的方法 - Google Patents

Abstract

一种针对序列间重叠检测的方法，具体涉及一种基于第三代测序技术ONT数据的序列间重叠检测的方法，本发明为了解决第三代测序技术在重叠检测时是全部序列进行比对，导致比对速度慢、内存消耗大，且第三代测序技术具有较高的错误率，产生的序列间重叠检测结果准确率不高，可靠性低的问题，它首先将获取的测序序列分为种子序列和映射序列；再构建种子序列的二级哈希索引表；计算种子序列的过滤阈值；最终利用种子序列的二级哈希索引表和过滤阈值将种子序列和映射序列进行比对，得到种子序列和映射序列的重叠信息，再以迭代的方式在每轮迭代开始时重新选取种子序列和映射序列，并重复上述操作，得到测序序列的重叠信息。属于序列重叠检测领域。

Description

一种针对序列间重叠检测的方法

技术领域

本发明涉及一种序列重叠检测方法，具体涉及一种基于第三代测序技术ONT数据的序列间重叠检测的方法，属于序列重叠检测领域。

背景技术

从头基因组拼接是指根据测序序列对物种基因组序列进行重建的过程。通过检测测序序列间的重叠关系，利用图模型发明，确定测序序列间的连接顺序，重构基因组上连续序列的过程，称为重叠群拼接。重叠群拼接是完成从头基因组拼接的第一步。测序序列间重叠检测发明是指给定一个测序序列的集合，通过测序序列间的相互比对找到序列间重叠位置的过程。该类发明为基因组序列拼接、基因组序列纠错等方法的上游工作，在生物信息学研究中具有重要地位。

与高通量测序技术(NGS)产生的测序序列(平均长度为100-250bp、测错率小于1％)相比，第三代测序技术(TGS)可以产生长度更长但错误率更高的测序序列(长度可达10kbp甚至100kbp以上、测错率约为12％-25％)。其中以PacBio公司的单分子实时测序技术(SMRT)和牛津纳米孔测序技术(ONT)为代表。第三代测序序列凭借其读长优势，能够覆盖基因组上多数大型结构变异和重复序列区域，因而可以大大提高检测结构变异的能力、比对的正确性和基因组拼接的连续性等。但其较高的错误率也对比对方法提出了更多的挑战，产生的序列间重叠检测结果准确率不高，可靠性低，且第三代测序技术在重叠检测时是全部序列进行比对，导致比对速度慢、内存消耗大。

发明内容

本发明为了解决第三代测序技术在重叠检测时是全部序列进行比对，导致比对速度慢、内存消耗大，且第三代测序技术具有较高的错误率，产生的序列间重叠检测结果准确率不高，可靠性低的问题，进而提出了一种针对序列间重叠检测的方法。

本发明采取的技术方案是：

它包括以下步骤：

S1、获取测序序列；

S2、将测序序列分为种子序列和映射序列，将种子序列和映射序列进行比对，得到测序序列的重叠信息，具体过程为：

S21、将测序序列作为映射序列，再按照序列长度由长到短进行排序，选取一条或多条排序后长度最长的测序序列作为种子序列；

S22、构建种子序列的二级哈希索引表，得到构建好的二级哈希索引表；

S23、计算种子序列的过滤阈值；

S24、利用种子序列的二级哈希索引表和过滤阈值对种子序列与映射序列进行比对，得到种子序列与映射序列的重叠集合；

S25、将种子序列与映射序列的重叠集合进行对称记录，得到每个重叠在种子序列和映射序列上对应的信息，根据得到的所述信息分别计算种子序列和映射序列上重叠区域的覆盖度；

S26、将所有覆盖度不足的测序序列作为映射序列，选取一定数量的覆盖度不足的测序序列作为种子序列，重复执行S21-S25，直至覆盖度不足的测序序列的减少数量低于阈值，得到测序序列的重叠信息。

进一步地，所述S22中构建种子序列的二级哈希索引表，得到构建好的二级哈希索引表，具体过程为：

S221、从测序序列中依次读取每条种子序列，得到每条种子序列的信息，将得到的每条种子序列的信息合并，得到种子序列的信息集合；

所有种子序列的碱基数为：

其中，M表示自定义的内存上限，单位为GB；

#reads表示测序序列总数；

top_n表示每条测序序列每轮最多保留的重叠信息数量；

t表示线程数；

B₁表示存储每条种子序列、映射序列标记信息和重叠信息所占字节数；

B₂表示存储每个具体的重叠信息所占的字节数；

B₃表示存储索引表中的一个元素所占的字节数；

S₁表示比对过程中多线程同时读入的情况下，S₁个碱基所占的字节数；

S₂表示索引过程中多线程同时读入的情况下，S₂个碱基所占的字节数；

S222、根据种子序列的信息集合构建种子序列的线性列表和指针列表，得到构建好的线性列表和指针列表；

S223、根据构建好的线性列表和指针列表构建种子序列的二级哈希索引表。

进一步地，所述S221中每条种子序列的信息均由四元组read(id,length,name,seq)表示，其中，id表示每条种子序列在测序序列中的顺位，且作为每条种子序列对应的标识符，length表示每条种子序列的长度，name表示每条种子序列在测序序列中的名称，seq表示每条种子序列具体的碱基序列。

进一步地，所述S222中根据种子序列的信息集合构建种子序列的线性列表和指针列表，得到构建好的线性列表和指针列表，具体过程为：

S2221、定义一个长度为w的窗口，所述窗口从每条种子序列的第一个碱基位置开始向右移动，直至最后一个碱基结束，利用所述窗口重复遍历每条种子序列，每个窗口内均有w个短字符串；

S2222、使用哈希函数分别处理每个窗口中的短字符串，得到每个窗口中每个短字符串的哈希值，选取每个窗口中哈希值最小的短字符串作为minimizer，令minimizer表示当前的窗口，当两个相邻的minimizer的哈希值一致时，两个minimizer位于同一条种子序列上的相同位置，则两个相邻窗口保留同一个minimizer，得到所有种子序列的minimizer，将所有minimizer根据得到的顺序进行排序，生成线性列表mi；

S2223、指针列表mi_num的每个元素均记录具有相同前l位的kmer的数量。

进一步地，每个minimizer均包括四元组minimizer(value,rid,pos,strand)，其中，value表示该minimizer对应的哈希值，rid表示种子序列标识符，pos表示该minimizer在种子序列上的位置，strand表示该minimizer对应的正负链。

进一步地，所述S223中根据构建好的线性列表和指针列表构建种子序列的二级哈希索引表，具体过程为：

S2231、根据指针列表mi_num计算得到指针列表mi_count；

S2232、根据每个minimizer对应的哈希值大小对线性列表mi进行排序，生成有序线性列表mi′；

S2233、将指针列表mi_count和有序线性列表mi′合并，组成种子序列的二级哈希索引表。

进一步地，所述S23中计算种子序列的过滤阈值，具体过程为：

计算种子序列的minimizer的过滤阈值：

rep_n＝rep_kmer[#minimizer*r_n] (2)

其中，rep_n表示过滤阈值；

#minimizer表示不同哈希值的minimizer的数量；

r_n表示自定义的过滤重复次数最多的minimizer的百分比；

rep_kmer[]表示存储每个不同minimizer数量的递减有序数组。

进一步地，所述S24中利用种子序列的二级哈希索引表和过滤阈值对种子序列与映射序列进行比对，得到种子序列与映射序列的重叠集合，具体过程为：

S241、从测序序列中依次读取每条映射序列，得到每条映射序列的信息，将得到的每条映射序列的信息合并，得到映射序列的信息集合；

所述每条映射序列的信息均由四元组read′(id′,length′,name′,seq′)表示，其中，id′表示每条映射序列在文件中的顺位，且作为每条映射序列对应的标识符，length′表示每条映射序列的长度，name′表示每条映射序列在文件中的名称，seq′表示每条映射序列具体的碱基序列；

所有映射序列的碱基数为：

其中，M表示自定义的内存上限，单位为GB；

#reads表示测序序列总数；

top_n表示每条测序序列每轮最多保留的重叠信息数量；

t表示线程数；

B₁表示存储每条种子序列、映射序列标记信息和重叠信息所占字节数；

B₂表示存储每个具体的重叠信息所占的字节数；

B₃表示存储索引表中的一个元素所占的字节数；

S₁表示比对过程中多线程同时读入的情况下，S₁个碱基所占的字节数；

S₂表示索引过程中多线程同时读入的情况下，S₂个碱基所占的字节数；

S242、根据得到的每条映射序列的信息分别构建每条映射序列的线性列表二，得到每条映射序列的线性列表二，具体过程为：

S2421、定义一个长度为w的窗口，所述窗口从每条映射序列的第一个碱基位置开始向右移动，直至最后一个碱基结束；利用所述窗口重复遍历每条映射序列，每个窗口内均有w个短字符串；

S2422、使用哈希函数处理每条映射序列上每个窗口中的短字符串，得到每个窗口中每个短字符串的哈希值，选取每个窗口中哈希值最小的短字符串作为minimizer′，令minimizer′表示当前的窗口，当两个相邻窗口的minimizer′的哈希值一致时，两个minimizer′位于同一条映射序列上的相同位置，则两个相邻窗口保留同一个minimizer′，得到所述映射序列的所有minimizer′，将所有minimizer′合并生成线性列表二mi_q；

每个minimizer′均包括四元组minimizer′(value′,rid′,pos′,strand′)，其中，value′表示该minimizer′对应的哈希值，rid′表示该minimizer′对应的映射序列标识符，pos′表示该minimizer′在映射序列上的位置，strand′表示该minimizer′对应的正负链；

S243、利用种子序列的二级哈希索引表查找每条映射序列和目标种子序列的匹配块，得到每条映射序列和目标种子序列的匹配块集合，具体过程为：

S2431、根据二级哈希索引表中的指针列表mi_count定位minimizer′的lmer在有序线性列表mi′中的起始位置和终止位置；

S2432、利用二分查找找到与minimizer′哈希值相同的minimizer在有序线性列表mi′中的位置pos；

S2433、再利用二分查找找到位置pos的上下两个方向该哈希值的起始位置和终止位置，得到映射序列和目标种子序列对应的每个匹配块，如果查找得到匹配块且匹配块的数量大于过滤阈值rep_n，则忽略此匹配块；否则将此匹配块加入匹配块集合中，得到每条映射序列和目标种子序列的初步匹配块集合；

S2434、根据初步匹配块集合中匹配块的种子序列标识符和匹配块在映射序列上的起始位置对匹配块进行排序，若排序后相邻的两个匹配块符合共线性条件，则合并为一个匹配块，得到每条映射序列和目标种子序列的匹配块集合；

所述共线性条件：

ts₁和ts₂表示相邻的两个匹配块分别在目标种子序列上的起始位置；

qs₁和qs₂表示相邻的两个匹配块分别在映射序列上的起始位置；

T表示基因组上相邻的两个可以合并的kmer的距离的最大值；

tid₁和tid₂表示相邻的两个匹配块对应的目标种子序列的标识符；

所述匹配块集合中的每个匹配块均用一个六元组MB(qs,qe,tid,ts,te,cov′)表示，其中，qs表示该匹配块在映射序列上的起始位置，qe表示该匹配块在映射序列上的结束位置，tid表示该匹配块对应的目标种子序列的标识符，ts表示该匹配块在目标种子序列上的起始位置，te表示该匹配块在目标种子序列上的结束位置，cov′表示匹配的不重叠的碱基数量；

S244、根据所述匹配块集合利用稀疏动态规划算法得到每条映射序列和目标种子序列的比对骨架，具体过程为：

以每条映射序列和目标种子序列的匹配块为顶点构建有向无环图，则有向无环图的任意两顶点MB₃(qs₃,qe₃,tid₃,ts₃,te₃,cov₃′)和MB₄(qs₄,qe₄,tid₄,ts₄,te₄,cov′₄)间存在边的条件：

根据有向无环图的边和递归公式(6)，计算有向无环图中所有路径的单独得分：

S(MB_j)＝max(S(MB_i)+w(MB_i→MB_j)-p(MB_i→MB_j)) (6)

其中，S表示当前顶点最优路径的得分；

w(MB_i→MB_j)表示边MB_i→MB_j的权值，

p(MB_i→MB_j)表示边MB_i→MB_j的惩罚系数，

p(MB_i→MB_j)＝abs((qs_j-qe_i)-(ts_j-te_i))/w(MB_i→MB_j)；

路径是将很多个匹配块连接，同一个映射序列可能与多个目标种子序列存在匹配，针对每一个存在匹配的目标种子序列以得分最高的路径作为每条映射序列和对应目标种子序列的比对骨架；

S245、根据动态规划矩阵和所述比对骨架得到每条映射序列和种子序列的重叠信息，过滤不符合条件的重叠，得到符合条件的重叠集合，具体过程为：

S2451、根据动态规划矩阵和所述比对骨架进行回溯，得到每条映射序列和每个目标种子序列的重叠信息，将所有重叠信息进行合并，得到每条映射序列和所有目标种子序列的重叠信息集合；所述动态规划矩阵是有向无环中每个顶点在所有路径中得分最高的得分值的组合矩阵；

S2452、根据条件公式过滤重叠，得到符合条件的重叠；

所述条件公式：

其中，matching bases表示映射序列和目标种子序列匹配块的不重叠碱基数；

T₁表示自定义的最小匹配碱基数；

matching length表示映射序列和目标种子序列得到的重叠的长度；

length′表示映射序列的长度；

ratio表示自定义的重叠区域与映射序列长度的比率；

ove_min表示自定义的最小重叠长度；

overhang表示映射序列和目标种子序列两端不匹配的碱基数；

T₂表示自定义的两端不匹配的最大碱基数；

S2453、对于每条映射序列，每轮迭代自定义保留一定数量的重叠骨架得分最高的重叠，重复执行S242-S245得到映射序列与种子序列间的所有重叠，将所有重叠合并，得到重叠集合。

进一步地，所述S2451得到的每个重叠信息为ove(qid,ql,qs′,qe′,tid,tl,ts′,te′,rev,mbp,mln,score)，其中，qid表示映射序列的标识符；ql表示映射序列的长度；tl表示目标种子序列的长度；qs′表示映射序列上与目标种子序列重叠的起始位置；qe′表示映射序列上与目标种子序列重叠的结束位置；ts′表示目标种子序列上与映射序列重叠的起始位置；te′表示目标种子序列上与映射序列重叠的结束位置；rev表示映射序列和目标种子序列的正负情况，若在同一条链上则为“+”，反之为“-”；mbp表示映射序列和目标种子序列重叠的碱基数量；mln表示映射序列和目标种子序列重叠的长度；score表示映射序列和目标种子序列重叠骨架的得分。

进一步地，所述S25中将种子序列与映射序列的重叠集合进行对称记录，得到每个重叠在种子序列和映射序列上对应的信息，根据得到的所述信息分别计算种子序列和映射序列上重叠区域的覆盖度，具体过程为：

S251、将S2453得到的重叠集合进行对称记录，得到包含重叠信息的种子序列和映射序列；

S252、将每个映射序列和每个种子序列均分为多个1000bp的块，以块为单位分别计算S251得到的每个映射序列和每个种子序列的覆盖度，具体过程为：

对于种子序列B，通过重叠集合找到与种子序列B直接连接的种子序列P和映射序列Q，P和Q的数量均大于等于1；根据重叠信息计算种子序列P和映射序列Q重叠区域的起始位置和终止位置所在块的坐标范围，若所述坐标范围在种子序列B的范围内，则将种子序列B范围内的覆盖度加一；若所述坐标范围超出了种子序列B的范围，则无覆盖度，继续查找下一条重叠信息；

对于映射序列G，通过重叠集合找到与映射序列G直接连接的种子序列Y，Y的数量大于等于1，根据重叠信息计算块的坐标范围，令在映射序列G的坐标范围内的块的覆盖度等于相连的种子序列Y的对应块的覆盖度，若映射序列G与多个种子序列Y相连，则取多个种子序列Y的平均覆盖度作为映射序列G和种子序列Y间每块的覆盖度，不在映射序列G的坐标范围内的块的覆盖度不计算。

有益效果：

(1)本发明将测序序列分为种子序列和映射序列，利用迭代的方法检测种子序列和映射序列间的重叠，每轮迭代的开始均需重新选择种子序列和映射序列，避免了所有序列间的一一比对，减少计算资源的消耗，能够更加快速、准确地得到所需重叠信息。能够在每轮覆盖基因组上的部分区域通过将映射序列比对到每轮的种子序列上，即能完成该基因组区域的测序序列间重叠检测。

在第一轮迭代的开始按照测序序列长度由长到短进行排序，将测序序列作为映射序列，选取一条或多条排序后长度最长的测序序列作为种子序列；构建种子序列的二级哈希索引表；计算种子序列的过滤阈值；利用种子序列的二级哈希索引表和过滤阈值对种子序列与映射序列进行比对，得到种子序列与映射序列的重叠集合，此时，第一轮迭代结束。在第二轮及以上轮次迭代的开始将上一轮得到的种子序列与映射序列的重叠集合进行对称记录，得到每个重叠在种子序列和映射序列上对应的信息，根据得到的所述信息分别计算种子序列和映射序列上重叠区域的覆盖度；选取覆盖度不足的序列作为本轮需要检测的映射序列，随机选取部分覆盖度不足的序列作为种子序列，重复第一轮的操作，直至覆盖度不足的测序序列的减少数量低于阈值，得到测序序列的重叠信息。

(2)通过最小哈希策略在每轮迭代时构建种子序列的二级哈希索引表，提高种子检索效率，实现对序列间精确匹配块的快速查找。再通过稀疏动态规划算法进行种子序列和映射序列间的比对，能够有效地提高比对的速度，可以同时检测出映射序列与所有种子序列之间的重叠信息。

(3)设计了控制内存消耗的运行策略，可由用户在运行时设定当前计算机可用内存的上限，本发明根据所述上限计算每轮或每批能够处理的碱基总数。该策略使得本发明内存可控，可以在个人计算机和服务器等多种计算机上运行。实现了快速检测、内存可控的序列间重叠检测。

附图说明

图1是本发明的整体流程图；

具体实施方式

具体实施方式一：结合图1说明本实施方式，本实施方式所述一种针对序列间重叠检测的方法，它包括以下步骤：

S1、获取测序序列。

获取测序序列的集合，所述集合中包括多条测序序列。

S2、将测序序列分为种子序列和映射序列，将种子序列和映射序列利用迭代的方法进行比对，得到测序序列的重叠信息，具体过程为：

S21、将测序序列作为映射序列，再按照序列长度由长到短进行排序，选取一条或多条排序后长度最长的测序序列作为种子序列。

在每轮迭代的开始，均需要重新确定种子序列和映射序列，所有测序序列存储在fastq或fasta格式的文件中。可以理解为，测序序列集合中的每条序列都有两个标记，分别标记在当前迭代轮次中是否为种子序列和是否为映射序列。所述一条或多条排序后长度最长的测序序列为选取排序后的测序序列前几位的测序序列，选取的测序序列可能等长，也可能不等长。根据迭代的轮次不同，选取种子序列和映射序列的方法可以分为第一轮迭代和第二轮及以上迭代两种情况。

在第一轮迭代时，所有测序序列都未被检测，所以将所有测序序列作为映射序列，再对所有测序序列的长度进行由长到短的排序，选取一定数量的(所述数量由用户定义)最长的测序序列作为种子序列，如此只是为了在第一轮迭代时选取长度最长的部分序列，将它们标记为种子序列。

S22、构建种子序列的二级哈希索引表，得到构建好的二级哈希索引表，具体过程为：

确定种子序列后，对选定的种子序列进行二级哈希索引表的构建，构建过程可以分为以下几个步骤：

S221、从fastq或fasta格式的测序序列中依次读取每条种子序列，得到每条种子序列的信息，将得到的每条种子序列的信息合并，得到种子序列的信息集合，所有种子序列的碱基数约为batch_size，根据用户提供的内存上限M计算得到，M单位为GB。

其中，M表示自定义的内存上限，单位为GB；#reads表示测序序列总数；top_n表示每条测序序列每轮最多保留的重叠信息数量；t表示线程数；B₁表示存储每条种子序列、映射序列标记信息和重叠信息所占字节数；B₂表示存储每个具体的重叠信息所占的字节数；B₃表示存储索引表中的一个元素所占的字节数；S₁表示比对过程中多线程同时读入的情况下，S₁个碱基所占的字节数；S₂表示索引过程中多线程同时读入的情况下，S₂个碱基所占的字节数。

每条种子序列均用四元组read(id,length,name,seq)表示，其中，id表示每条种子序列在测序序列中的顺位，且作为每条种子序列对应的标识符，length表示每条种子序列的长度，name表示每条种子序列在测序序列中的名称，seq表示每条种子序列具体的碱基序列。

S222、根据种子序列的信息集合构建种子序列的线性列表和指针列表，得到构建好的线性列表和指针列表，具体过程为：

S2221、定义一个长度为w的窗口，所述窗口从每条种子序列的第一个碱基位置开始向右移动，直至最后一个碱基结束，利用所述窗口重复遍历每条种子序列，每个窗口内均有w个短字符串kmer。

窗口长度为w，表示其中有w个短字符串，分别作为w个短字符串kmer的开头第一个字符。窗口每向右移动一次，第一个短字符串kmer就会被移出，新的短字符串kmer被补入，所以窗口中一直有w个短字符串。

S2222、使用哈希函数分别处理每个窗口中的短字符串kmer，得到每个窗口中每个短字符串kmer的哈希值，选取每个窗口中哈希值最小的短字符串kmer作为minimizer，令minimizer表示当前的窗口，当两个相邻的minimizer的哈希值一致时，两个minimizer位于同一条种子序列上的相同位置，则两个相邻窗口保留同一个minimizer，得到所有种子序列的minimizer，将所有minimizer根据得到的顺序进行排序，生成线性列表mi，线性列表mi为长度不定的线性列表，用于存储所有种子序列的minimizer。

每个minimizer均用四元组minimizer(value,rid,pos,strand)表示，其中，value表示该minimizer对应的哈希值，rid表示该minimizer对应的种子序列标识符，pos表示该minimizer在种子序列上的位置，strand表示该minimizer对应的正负链。

S2223、指针列表mi_num的每个元素均记录具有相同前l位的kmer的数量。每得到一个表示当前窗口的minimizer，根据minimizer的哈希值，取所述哈希值的前l位lmer，lmer为长度为l的短字符串，对应指针列表mi_num的第lmer个元素加一。

S223、根据构建好的线性列表和指针列表构建种子序列的二级哈希索引表，具体过程为：

S2231、根据指针列表mi_num计算得到指针列表mi_count。

构建的二级哈希索引表表由列表mi和指针列表mi_count共同组成，指针列表mi_count用于记录所有前l位为相应键值的kmer在线性列表mi中的起始位置，指针列表mi_count是长度为4^l的哈希表，指针列表mi_count的键值为lmer(l＜k)，元素的数据类型为无符号整型。指针列表mi_count由指针列表mi_num直接计算得到。

S2232、根据每个minimizer对应的哈希值大小对线性列表mi进行排序，生成有序线性列表mi′。

S2233、将指针列表mi_count和有序线性列表mi′合并，组成种子序列的二级哈希索引表。

S23、计算种子序列的minimizer的过滤阈值，具体过程为：

rep_n＝rep_kmer[#minimizer*r_n] (2)

式中，rep_n表示过滤阈值；#minimizer表示不同哈希值的minimizer的数量；r_n表示过滤重复次数最多的minimizer的百分比(可由用户通过参数r定义)；rep_kmer[]表示存储每个不同minimizer数量的递减有序数组。

为了提高比对的准确性，本发明根据以上公式计算得到过滤阈值rep_n。根据有序线性列表mi′得到不同哈希值的minimizer的数量(将相同哈希值的minimizer划为一组，得到共有多少组不同的minimizer，每组minimizer的数量是多少)和每个minimizer出现的次数，如此是为了得到一个阈值，在下述过程过滤出现次数多的minimizer，即出现次数大于过滤阈值rep_n的minimizer在下述比对过程中不予考虑，出现次数过于频繁的minimizer可能来自基因组上的重复序列区域。

S24、利用种子序列的二级哈希索引表和过滤阈值对种子序列与映射序列进行比对，得到种子序列与映射序列的重叠集合，具体过程为：

S241、从fastq或fasta格式的测序序列中依次读取每条映射序列，得到每条映射序列的信息，将得到的每条映射序列的信息合并，得到映射序列的信息集合，所有映射序列的碱基数约为batch_size，根据用户提供的内存上限M计算得到，M单位为GB。

其中，M表示自定义的内存上限，单位为GB；#reads表示测序序列总数；top_n表示每条测序序列每轮最多保留的重叠信息数量；t表示线程数；B₁表示存储每条种子序列、映射序列标记信息和重叠信息所占字节数；B₂表示存储每个具体的重叠信息所占的字节数；B₃表示存储索引表中的一个元素所占的字节数；S₁表示比对过程中多线程同时读入的情况下，S₁个碱基所占的字节数；S₂表示索引过程中多线程同时读入的情况下，S₂个碱基所占的字节数。

每条映射序列均由四元组read′(id′,length′,name′,seq′)表示，其中，id′表示每条映射序列在文件中的顺位，且作为每条映射序列对应的标识符，length′表示每条映射序列的长度，name′表示每条映射序列在文件中的名称，seq′表示每条映射序列具体的碱基序列。

S242、根据得到的每条映射序列的信息分别构建每条映射序列的线性列表二，得到每条映射序列的线性列表二，具体过程为：

S2421、定义一个长度为w的窗口，所述窗口从每条映射序列的第一个碱基位置开始向右移动，直至最后一个碱基结束；利用所述窗口重复遍历每条映射序列，每个窗口内均有w个短字符串kmer。

窗口长度为w，表示其中有w个短字符串，分别作为w个短字符串kmer的开头第一个字符。窗口每向右移动一次，第一个短字符串kmer就会被移出，新的短字符串kmer被补入，所以窗口中一直有w个短字符串。

S2422、使用哈希函数处理每条映射序列上每个窗口中的短字符串kmer，得到每个窗口中每个短字符串kmer的哈希值，选取每个窗口中哈希值最小的短字符串kmer作为minimizer′，令minimizer′表示当前的窗口，当两个相邻窗口的minimizer′的哈希值一致时，两个minimizer′位于同一条映射序列上的相同位置，则两个相邻窗口保留同一个minimizer′，得到所述映射序列的所有minimizer′，将所有minimizer′合并生成线性列表二mi_q，线性列表二mi_q为长度不定的线性列表，用于存储每个映射序列的minimizer′。

每个minimizer′均包括四元组minimizer′(value′,rid′,pos′,strand′)，其中，value′表示该minimizer′对应的哈希值，rid′表示该minimizer′对应的映射序列标识符，pos′表示该minimizer′在映射序列上的位置，strand′表示该minimizer′对应的正负链。

S243、利用种子序列的二级哈希索引表查找每条映射序列和目标种子序列的匹配块，得到每条映射序列和目标种子序列的匹配块集合，具体过程为：

匹配块长度很小，和上文kmer长度是一致的。将每个线性列表二mi_q中每个minimizer′在索引表中进行二分查找。根据二级哈希索引表中指针列表mi_count定位minimizer′的前lmer在有序线性列表mi′中的起始位置和终止位置，再通过二分查找找到与minimizer′哈希值相同的minimizer在有序线性列表mi′中的位置pos，最后通过第二次二分查找找到位置pos的上下两个方向该哈希值的起始位置和终止位置(相同哈希值的minimizer可能有多个)，得到映射序列和目标种子序列对应的每个匹配块，如果查找得到匹配块且匹配块的数量大于过滤阈值rep_n，则忽略此匹配块；否则将此匹配块加入匹配块集合中，如此重复上述操作，最终得到每条映射序列和目标种子序列的初步匹配块集合。由于种子序列和映射序要一一对比，所以目标种子序列为当前和映射序列对比的种子序列。

将上述的初步匹配块集合中的每个匹配块(EMs，ExactMatches)均用一个四元组EM(tid,qs,ts,cov)表示，其中，tid表示目标种子序列的标识符，qs表示该匹配块在映射序列上的起始位置，ts表示该匹配块在目标种子序列上的起始位置，cov表示匹配块的长度，即匹配的碱基数量。

根据初步匹配块集合中匹配块的种子序列标识符和匹配块在映射序列上的起始位置对匹配块进行排序，若排序后相邻匹配块EM₁(tid₁,qs₁,ts₁,cov₁)和EM₂(tid₂,qs₂,ts₂,cov₂)符合下式的共线性条件，则合并为一个匹配块(MB，Match Block)，如此可得到每条映射序列和目标种子序列的最终匹配块集合，所述匹配块集合中的每个匹配块均用一个六元组MB(qs,qe,tid,ts,te,cov′)表示，其中，qe表示该匹配块在映射序列上的结束位置，te表示该匹配块在目标种子序列上的结束位置，cov′表示匹配的不重叠的碱基数量。

T表示基因组上两个相邻的可以合并的kmer的距离的最大值，可由现有方法计算得到；

S244、根据所述匹配块集合利用稀疏动态规划算法得到每条映射序列和目标种子序列的比对骨架，具体过程为：

根据得到的匹配块集合，使用稀疏动态规划算法得到每条映射序列和目标种子序列间的比对骨架，具体步骤如下：

以每条映射序列和目标种子序列的匹配块MB为顶点构建有向无环图DAG，则有向无环图DAG的任意两顶点MB₃(qs₃,qe₃,tid₃,ts₃,te₃,cov′₃)和MB₄(qs₄,qe₄,tid₄,ts₄,te₄,cov′₄)间存在边的条件如下式：

根据有向无环图的边和递归公式(6)，计算有向无环图DAG中所有路径的单独得分，路径是将很多个匹配块连接：

S(MB_j)＝max(S(MB_i)+w(MB_i→MB_j)-p(MB_i→MB_j)) (6)

其中，S表示当前顶点最优路径的得分，w(MB_i→MB_j)表示边MB_i→MB_j的权值，p(MB_i→MB_j)表示边MB_i→MB_j的惩罚系数。

p(MB_i→MB_j)＝abs((qs_j-qe_i)-(ts_j-te_i))/w(MB_i→MB_j) (8)

由于同一个映射序列可能与多个目标种子序列存在匹配，所以针对每一个存在匹配的目标种子序列，以得分最高的路径作为每条映射序列和所述目标种子序列的比对骨架。

S245、根据动态规划矩阵和所述比对骨架得到每条映射序列和种子序列的重叠信息，过滤不符合条件的重叠，得到符合条件的重叠集合，具体过程为：

动态规划矩阵是在计算所有路径得分时得到的，动态规划矩阵中的一个元素是该顶点在所有路径中得分最高的得分值，即动态规划矩阵是有向无环中每个顶点在所有路径中得分最高的得分值的组合矩阵。

S2451、根据动态规划矩阵和所述比对骨架(路径的最高得分值)进行回溯，得到每条映射序列和每个目标种子序列的重叠信息，将所有重叠信息进行合并，得到每条映射序列和所有目标种子序列的重叠信息集合；每个重叠信息表示为ove(qid,ql,qs′,qe′,tid,tl,ts′,te′,rev,mbp,mln,score)，其中，qid表示映射序列的标识符；ql表示映射序列的长度；tl表示目标种子序列的长度；qs′表示映射序列上与目标种子序列重叠的起始位置；qe′表示映射序列上与目标种子序列重叠的结束位置；ts′表示目标种子序列上与映射序列重叠的起始位置；te′表示目标种子序列上与映射序列重叠的结束位置；rev表示映射序列和目标种子序列的正负情况，若在同一条链上则为“+”，反之为“-”；mbp表示映射序列和目标种子序列重叠的碱基数量；mln表示映射序列和目标种子序列重叠的长度；score表示映射序列和目标种子序列重叠骨架的得分。

S2452、根据条件公式过滤重叠，得到符合条件的重叠，当重叠信息ove满足下列条件时保留，否则丢弃，条件公式为：

其中，matching bases表示映射序列和目标种子序列匹配块的不重叠碱基数；T₁表示可由用户定义的最小匹配碱基数；matching length表示映射序列和目标种子序列得到的重叠的长度；length′表示映射序列的长度；ratio表示由用户定义的重叠区域与映射序列长度的比率；ove_min为可由用户定义的最小重叠长度；overhang表示映射序列和目标种子序列两端不匹配的碱基数；T₂为可由用户定义的两端不匹配的最大碱基数。最终得到符合条件的重叠集合ove_collection。

S2453、对于每条映射序列，每轮迭代保留n个(n可由用户定义)重叠集合ove_collection中重叠骨架得分最高的重叠，获得对应的重叠信息ove。重复执行S242-S245得到映射序列与种子序列间的所有重叠，将所有重叠合并，得到重叠集合。

S25、将种子序列与映射序列的重叠集合进行对称记录，得到每个重叠在种子序列和映射序列上对应的信息，根据得到的所述信息分别计算种子序列和映射序列上重叠区域的覆盖度，具体过程为：

在第二轮及以上迭代时，第一轮迭代已经检测过部分测序序列间的重叠信息。本发明通过随机选取种子序列，能够在每轮覆盖基因组上的部分区域通过将映射序列比对到每轮的种子序列上，即能完成该基因组区域的测序序列间重叠检测。同时映射序列间的重叠信息可以通过种子序列传递计算得出。如果来自基因组某区域的序列没有被选取作为种子序列，那么同样来自该区域的种子序列和映射序列的覆盖度会低于该测序序列集合的平均覆盖度甚至为零。故本发明通过计算所有序列上重叠区域的覆盖度，选取覆盖度不足的序列进入第二轮及以上轮数比对。满足以下条件的序列r称为覆盖度不足的序列：序列r中m(m为程序参数，可由用户定义，默认值等于3)个以上连续块的平均覆盖度小于某个阈值T(T为程序参数，可由用户定义，为测序序列集合的覆盖度的百分比)。

S251、在第二轮及以上轮次迭代时估计覆盖度前，先将上一轮得到的重叠集合ove_collection进行对称记录，得到包含重叠信息的种子序列和映射序列，若标识符为qid的映射序列能够比对到标识符为tid的种子序列上，即重叠集合ove_collection中存在记录重叠信息ove(qid,ql,qs′,qe′,tid,tl,ts′,te′,rev,mbp,mln,score)表示映射序列和目标种子序列间的重叠信息，那么在重叠集合中ove_collection记录其对称的重叠信息ove(qid,ql,qs′,qe′,tid,tl,ts′,te′,rev,mbp,mln,score)，使得映射序列和种子序列可以双向查找到对方。

S252、为了避免对计算资源的浪费，本发明不对覆盖度进行碱基水平的估计，而是将每个映射序列和每个种子序列均分为多个1000bp的块，以块为单位分别计算S251得到的每个映射序列和每个种子序列的覆盖度。

对于种子序列B，通过以下方法计算其覆盖度：通过重叠集合ove_collection找到与种子序列B直接连接的种子序列P和映射序列Q，P和Q的数量均大于等于1；根据重叠信息计算种子序列P和映射序列Q重叠区域的起始位置和终止位置所在块的坐标范围，若所述坐标范围在种子序列B的范围内，则将种子序列B范围内的覆盖度加一；若所述坐标范围超出了种子序列B的范围，则无覆盖度，继续查找下一条重叠信息。

对于与种子序列B相连且在种子序列B范围内的种子序列P，使用堆栈(辅助结构)查找能找到通过种子序列B传递的与种子序列P不直接相连但存在重叠区域的映射序列。

对于映射序列G，通过以下方法计算其覆盖度：假设检测出的种子序列和映射序列间的重叠都是正确的，则任意两条序列间的重叠区域都为基因组上相同的区域，故重叠区域的覆盖度相等。通过重叠集合ove_collection找到与映射序列G直接连接的种子序列Y，Y的数量大于等于1，根据重叠信息计算块的坐标范围，令在映射序列G的坐标范围内的块的覆盖度等于相连的种子序列Y的对应块的覆盖度，若映射序列G与多个种子序列Y相连，则取多个种子序列Y的平均覆盖度作为映射序列G和种子序列Y间每块的覆盖度，不在映射序列G的坐标范围内的块的覆盖度不计算。

在实际情况中，由于测序错误和基因组上重复序列的存在，可能检测出错误的序列间重叠信息，此时来自基因组上不同区域的序列通过该错误重叠信息的传递连接在一起，导致计算得到的相应的序列的覆盖度高于实际覆盖度。

S26、将所有覆盖度不足的测序序列作为映射序列，选取一定数量(所述数量由用户定义)的覆盖度不足的测序序列作为种子序列，重复执行S21-S25，直至某轮迭代中减少的覆盖度不足的测序序列的数量低于某个阈值，得到测序序列的重叠信息。

所有迭代轮次结束后，序列间重叠信息的结果将以PAF格式输出至硬盘。

本发明在应用时通过双层并行架构实现多个线程。以三个主线程为例，该架构根据计算的需求，首先生成三个主线程，用于以下功能：

步骤一：从硬盘中读取数据至线程缓冲区。

步骤二：进行所需计算过程，将计算结果暂存至线程缓冲区。

步骤三：将线程缓冲区的计算结果写入硬盘。

三个主线程分别依次执行步骤一二三，且同一时刻不能执行相同的步骤。如果一个主线程某个步骤执行完毕，而其下一个要执行的步骤正在被其他线程占用，则会等待其他线程执行完毕再执行。每个线程执行完步骤三后，转至重新执行步骤一，重复执行步骤一至步骤三，直至将所有数据读取完毕。

该多线程的实现方法可以将读取数据、计算数据等步骤分离开同时进行，避免由于硬盘I/O速度过慢导致的整体速度太慢。

本发明的不同功能实现多线程的具体方法如下。

(1)构建索引的多线程方法

首先生成两个主线程，每个主线程都有两个步骤：

步骤一：从文件中读取种子序列。

步骤二：索引种子序列，构建二级哈希索引表表。

(2)比对的多线程方法

首先生成三个主线程，每个主线程都有三个步骤：

步骤一：从文件中读取映射序列；分配该主线程所需的内存空间。

步骤二：进行映射序列的对比。该步骤执行时，主线程生成t个子线程(该参数可由用户定义)，每个子线程单独处理一条映射序列，处理完毕后分配下一条。

步骤三：释放该主线程使用的内存空间。

(3)计算覆盖度的多线程方法

计算覆盖度时只使用一个主线程，主线程生成t个子线程(该参数可由用户定义)，每个子线程单独处理一条映射序列，处理完毕后分配下一条。

实施例

为了测试本发明在不同基因组序列上检测序列间重叠的性能，使用PBSIM软件模拟了8种不同物种(大肠杆菌、酵母、秀丽隐杆线虫、拟南芥、果蝇、玉米、小鼠、人类)的测序深度为50x、平均长度为15000、测序错误率为13％(***错误、删除错误和替换错误的比例为23：31：46)的R103 ONT数据集，并下载了4种不同物种(大肠杆菌、秀丽隐杆线虫、果蝇、人类)的不同测序深度的真实ONT数据集。模拟物种的基因组信息及生成的数据集信息如下表1(a)，真实数据集信息如下表1(b)。

表1(a):

表1(b)

/>

本次实验从序列水平的比对准确性、召回率两个方面将本发明和minimap2工具进行比较。根据程序运行过程中所需时间和内存峰值来评估本发明对计算资源的耗费程度。各评价标准定义如下。

(1)准确性

序列间重叠位置正确需要满足以下条件：即根据两序列在参考基因组上的真实位置计算得到重叠在参考基因组上的起始位置和结束位置，若两起始位置和结束位置的差值在一定范围内，则认为该重叠是正确的。

其中，

分别为两序列间的重叠在参考基因组上的起始位置和结束位置，seq_length为序列长度，error_rate为测序错误率。

假定满足以上条件的重叠信息数目为num_correct，输出的重叠信息数目为num_reported，则序列间重叠的总体准确性定义为：

(2)召回率

由于本发明对于每条序列只保留得分最高的重叠，每条序列与其他序列间的全部重叠信息可以通过传递的方法进行计算。假定经过上述传递操作计算得到的重叠数目为num_transitive，序列间真实重叠数目为num_true，序列间重叠的总体召回率定义为：

实验结果

实验利用本发明和minimap2对8个模拟数据集进行重叠检测，实验结果如下表2，统计了程序运行所需时间和内存峰值，以及程序得到的结果的准确性和召回率。

表2

程序得到的结果的准确性和召回率如下表3，其中，由于minimap2致力于得到所有序列间的重叠信息，故不根据上节所述计算传递重叠信息得到召回率的方法计算，而是使用其直接得到的重叠信息。

表3

实验利用本发明和minimap2对4个真实数据集进行重叠检测，实验结果如下表4，统计了程序运行所需时间和内存峰值，以及程序得到的结果的准确性和召回率

表4

/>

程序得到的结果的准确性和召回率如下表5。

表5

由上述表格2、表格4可以看出，本发明提出的序列间重叠检测方法的检测的比对效率在所有物种的测序数据上都超过minimap2，为该工具的3-6倍；内存消耗为minimap2的10％-20％。

由上述表格3、表格5可以看出，在重复序列较少的小规模基因组上，本发明的检测准确率可以达到在99％以上；在模拟数据集上的召回率达到98％以上，在真实数据集上的召回率在87％以上，比检测工具minimap2低1％-12％；而在重复序列较多的较大规模基因组上，检测准确率可以维持在90％以上，远超其他的检测工具，而召回率在94％以上。

Claims (9)

1.一种针对序列间重叠检测的方法，其特征在于：它包括以下步骤：

S1、获取测序序列；

S2、将测序序列分为种子序列和映射序列，将种子序列和映射序列进行比对，得到测序序列的重叠信息，具体过程为：

S21、将测序序列作为映射序列，再按照序列长度由长到短进行排序，选取一条或多条排序后长度最长的测序序列作为种子序列；

S22、构建种子序列的二级哈希索引表，得到构建好的二级哈希索引表；

S23、计算种子序列的过滤阈值；

S24、利用种子序列的二级哈希索引表和过滤阈值对种子序列与映射序列进行比对，得到种子序列与映射序列的重叠集合，具体过程为：

S241、从测序序列中依次读取每条映射序列，得到每条映射序列的信息，将得到的每条映射序列的信息合并，得到映射序列的信息集合；

所述每条映射序列的信息均由四元组read′(id′,length′,name′,seq′)表示，其中，id′表示每条映射序列在文件中的顺位，且作为每条映射序列对应的标识符，length′表示每条映射序列的长度，name′表示每条映射序列在文件中的名称，seq′表示每条映射序列具体的碱基序列；

所有映射序列的碱基数为：

其中，M表示自定义的内存上限，单位为GB；

#reads表示测序序列总数；

top_n表示每条测序序列每轮最多保留的重叠信息数量；

t表示线程数；

B₁表示存储每条种子序列、映射序列标记信息和重叠信息所占字节数；

B₂表示存储每个具体的重叠信息所占的字节数；

B₃表示存储索引表中的一个元素所占的字节数；

S₁表示比对过程中多线程同时读入的情况下，S₁个碱基所占的字节数；

S₂表示索引过程中多线程同时读入的情况下，S₂个碱基所占的字节数；

S242、根据得到的每条映射序列的信息分别构建每条映射序列的线性列表二，得到每条映射序列的线性列表二，具体过程为：

S2421、定义一个长度为w的窗口，所述窗口从每条映射序列的第一个碱基位置开始向右移动，直至最后一个碱基结束；利用所述窗口重复遍历每条映射序列，每个窗口内均有w个短字符串；

S2422、使用哈希函数处理每条映射序列上每个窗口中的短字符串，得到每个窗口中每个短字符串的哈希值，选取每个窗口中哈希值最小的短字符串作为minimizer′，令minimizer′表示当前的窗口，当两个相邻窗口的minimizer′的哈希值一致时，两个minimizer′位于同一条映射序列上的相同位置，则两个相邻窗口保留同一个minimizer′，得到所述映射序列的所有minimizer′，将所有minimizer′合并生成线性列表二mi_q；

每个minimizer′均包括四元组minimizer′(value′,rid′,pos′,strand′)，其中，value′表示该minimizer′对应的哈希值，rid′表示该minimizer′对应的映射序列标识符，pos′表示该minimizer′在映射序列上的位置，strand′表示该minimizer′对应的正负链；

S243、利用种子序列的二级哈希索引表查找每条映射序列和目标种子序列的匹配块，得到每条映射序列和目标种子序列的匹配块集合，具体过程为：

S2431、根据二级哈希索引表中的指针列表mi_count定位minimizer′的lmer在有序线性列表mi′中的起始位置和终止位置；

S2432、利用二分查找找到与minimizer′哈希值相同的minimizer在有序线性列表mi′中的位置pos；

S2433、再利用二分查找找到位置pos的上下两个方向该哈希值的起始位置和终止位置，得到映射序列和目标种子序列对应的每个匹配块，如果查找得到匹配块且匹配块的数量大于过滤阈值rep_n，则忽略此匹配块；否则将此匹配块加入匹配块集合中，得到每条映射序列和目标种子序列的初步匹配块集合；

S2434、根据初步匹配块集合中匹配块的种子序列标识符和匹配块在映射序列上的起始位置对匹配块进行排序，若排序后相邻的两个匹配块符合共线性条件，则合并为一个匹配块，得到每条映射序列和目标种子序列的匹配块集合；

所述共线性条件：

ts₁和ts₂表示相邻的两个匹配块分别在目标种子序列上的起始位置；

qs₁和qs₂表示相邻的两个匹配块分别在映射序列上的起始位置；

T表示基因组上相邻的两个可以合并的kmer的距离的最大值；

tid₁和tid₂表示相邻的两个匹配块对应的目标种子序列的标识符；

所述匹配块集合中的每个匹配块均用一个六元组MB(qs,qe,tid,ts,te,cov′)表示，其中，qs表示该匹配块在映射序列上的起始位置，qe表示该匹配块在映射序列上的结束位置，tid表示该匹配块对应的目标种子序列的标识符，ts表示该匹配块在目标种子序列上的起始位置，te表示该匹配块在目标种子序列上的结束位置，cov′表示匹配的不重叠的碱基数量；

S244、根据所述匹配块集合利用稀疏动态规划算法得到每条映射序列和目标种子序列的比对骨架，具体过程为：

以每条映射序列和目标种子序列的匹配块为顶点构建有向无环图，则有向无环图的任意两顶点MB₃(qs₃,qe₃,tid₃,ts₃,te₃,cov′₃)和MB₄(qs₄,qe₄,tid₄,ts₄,te₄,cov′₄)间存在边的条件：

根据有向无环图的边和递归公式(6)，计算有向无环图中所有路径的单独得分：

S(MB_j)＝max(S(MB_i)+w(MB_i→MB_j)-p(MB_i→MB_j)) (6)

其中，S表示当前顶点最优路径的得分；

w(MB_i→MB_j)表示边MB_i→MB_j的权值，

p(MB_i→MB_j)表示边MB_i→MB_j的惩罚系数，p(MB_i→MB_j)＝abs((qs_j-qe_i)-(ts_j-te_i))/w(MB_i→MB_j)；

路径是将很多个匹配块连接，同一个映射序列可能与多个目标种子序列存在匹配，针对每一个存在匹配的目标种子序列以得分最高的路径作为每条映射序列和对应目标种子序列的比对骨架；

S245、根据动态规划矩阵和所述比对骨架得到每条映射序列和种子序列的重叠信息，过滤不符合条件的重叠，得到符合条件的重叠集合，具体过程为：

S2451、根据动态规划矩阵和所述比对骨架进行回溯，得到每条映射序列和每个目标种子序列的重叠信息，将所有重叠信息进行合并，得到每条映射序列和所有目标种子序列的重叠信息集合；所述动态规划矩阵是有向无环中每个顶点在所有路径中得分最高的得分值的组合矩阵；

S2452、根据条件公式过滤重叠，得到符合条件的重叠；

所述条件公式：

其中，matchingbases表示映射序列和目标种子序列匹配块的不重叠碱基数；

T₁表示自定义的最小匹配碱基数；

matchinglength表示映射序列和目标种子序列得到的重叠的长度；

length′表示映射序列的长度；

ratio表示自定义的重叠区域与映射序列长度的比率；

ove_min表示自定义的最小重叠长度；

overhang表示映射序列和目标种子序列两端不匹配的碱基数；

T₂表示自定义的两端不匹配的最大碱基数；

S2453、对于每条映射序列，每轮迭代自定义保留一定数量的重叠骨架得分最高的重叠，重复执行S242-S245得到映射序列与种子序列间的所有重叠，将所有重叠合并，得到重叠集合；

S25、将种子序列与映射序列的重叠集合进行对称记录，得到每个重叠在种子序列和映射序列上对应的信息，根据得到的所述信息分别计算种子序列和映射序列上重叠区域的覆盖度；

S26、将所有覆盖度不足的测序序列作为映射序列，选取一定数量的覆盖度不足的测序序列作为种子序列，重复执行S21-S25，直至覆盖度不足的测序序列的减少数量低于阈值，得到测序序列的重叠信息。

2.根据权利要求1中所述的一种针对序列间重叠检测的方法，其特征在于：所述S22中构建种子序列的二级哈希索引表，得到构建好的二级哈希索引表，具体过程为：

S221、从测序序列中依次读取每条种子序列，得到每条种子序列的信息，将得到的每条种子序列的信息合并，得到种子序列的信息集合；

所有种子序列的碱基数为：

其中，M表示自定义的内存上限，单位为GB；

#reads表示测序序列总数；

top_n表示每条测序序列每轮最多保留的重叠信息数量；

t表示线程数；

B₁表示存储每条种子序列、映射序列标记信息和重叠信息所占字节数；

B₂表示存储每个具体的重叠信息所占的字节数；

B₃表示存储索引表中的一个元素所占的字节数；

S₁表示比对过程中多线程同时读入的情况下，S₁个碱基所占的字节数；

S₂表示索引过程中多线程同时读入的情况下，S₂个碱基所占的字节数；

S222、根据种子序列的信息集合构建种子序列的线性列表和指针列表，得到构建好的线性列表和指针列表；

S223、根据构建好的线性列表和指针列表构建种子序列的二级哈希索引表。

3.根据权利要求1中所述的一种针对序列间重叠检测的方法，其特征在于：所述S221中每条种子序列的信息均由四元组read(id,length,name,seq)表示，其中，id表示每条种子序列在测序序列中的顺位，且作为每条种子序列对应的标识符，length表示每条种子序列的长度，name表示每条种子序列在测序序列中的名称，seq表示每条种子序列具体的碱基序列。

4.根据权利要求3中所述的一种针对序列间重叠检测的方法，其特征在于：所述S222中根据种子序列的信息集合构建种子序列的线性列表和指针列表，得到构建好的线性列表和指针列表，具体过程为：

S2221、定义一个长度为w的窗口，所述窗口从每条种子序列的第一个碱基位置开始向右移动，直至最后一个碱基结束，利用所述窗口重复遍历每条种子序列，每个窗口内均有w个短字符串；

S2222、使用哈希函数分别处理每个窗口中的短字符串，得到每个窗口中每个短字符串的哈希值，选取每个窗口中哈希值最小的短字符串作为minimizer，令minimizer表示当前的窗口，当两个相邻的minimizer的哈希值一致时，两个minimizer位于同一条种子序列上的相同位置，则两个相邻窗口保留同一个minimizer，得到所有种子序列的minimizer，将所有minimizer根据得到的顺序进行排序，生成线性列表mi；

S2223、指针列表mi_num的每个元素均记录具有相同前l位的kmer的数量。

5.根据权利要求4中所述的一种针对序列间重叠检测的方法，其特征在于：每个minimizer均包括四元组minimizer(value,rid,pos,strand)，其中，value表示该minimizer对应的哈希值，rid表示种子序列标识符，pos表示该minimizer在种子序列上的位置，strand表示该minimizer对应的正负链。

6.根据权利要求5中所述的一种针对序列间重叠检测的方法，其特征在于：所述S223中根据构建好的线性列表和指针列表构建种子序列的二级哈希索引表，具体过程为：

S2231、根据指针列表mi_num计算得到指针列表mi_count；

S2232、根据每个minimizer对应的哈希值大小对线性列表mi进行排序，生成有序线性列表mi′；

S2233、将指针列表mi_count和有序线性列表mi′合并，组成种子序列的二级哈希索引表。

7.根据权利要求6中所述的一种针对序列间重叠检测的方法，其特征在于：所述S23中计算种子序列的过滤阈值，具体过程为：

计算种子序列的minimizer的过滤阈值：

rep_n＝rep_kmer[#minimizer*r_n] (2)

其中，rep_n表示过滤阈值；

#minimizer表示不同哈希值的minimizer的数量；

r_n表示自定义的过滤重复次数最多的minimizer的百分比；

rep_kmer[]表示存储每个不同minimizer数量的递减有序数组。

8.根据权利要求7中所述的一种针对序列间重叠检测的方法，其特征在于：所述S2451得到的每个重叠信息为ove(qid,ql,qs′,qe′,tid,tl,ts′,te′,rev,mbp,mln,score)，其中，qid表示映射序列的标识符；ql表示映射序列的长度；tl表示目标种子序列的长度；qs′表示映射序列上与目标种子序列重叠的起始位置；qe′表示映射序列上与目标种子序列重叠的结束位置；ts′表示目标种子序列上与映射序列重叠的起始位置；te′表示目标种子序列上与映射序列重叠的结束位置；rev表示映射序列和目标种子序列的正负情况，若在同一条链上则为“+”，反之为“-”；mbp表示映射序列和目标种子序列重叠的碱基数量；mln表示映射序列和目标种子序列重叠的长度；score表示映射序列和目标种子序列重叠骨架的得分。

9.根据权利要求8中所述的一种针对序列间重叠检测的方法，其特征在于：所述S25中将种子序列与映射序列的重叠集合进行对称记录，得到每个重叠在种子序列和映射序列上对应的信息，根据得到的所述信息分别计算种子序列和映射序列上重叠区域的覆盖度，具体过程为：

S251、将S2453得到的重叠集合进行对称记录，得到包含重叠信息的种子序列和映射序列；

S252、将每个映射序列和每个种子序列均分为多个1000bp的块，以块为单位分别计算S251得到的每个映射序列和每个种子序列的覆盖度，具体过程为：

对于种子序列B，通过重叠集合找到与种子序列B直接连接的种子序列P和映射序列Q，P和Q的数量均大于等于1；根据重叠信息计算种子序列P和映射序列Q重叠区域的起始位置和终止位置所在块的坐标范围，若所述坐标范围在种子序列B的范围内，则将种子序列B范围内的覆盖度加一；若所述坐标范围超出了种子序列B的范围，则无覆盖度，继续查找下一条重叠信息；

对于映射序列G，通过重叠集合找到与映射序列G直接连接的种子序列Y，Y的数量大于等于1，根据重叠信息计算块的坐标范围，令在映射序列G的坐标范围内的块的覆盖度等于相连的种子序列Y的对应块的覆盖度，若映射序列G与多个种子序列Y相连，则取多个种子序列Y的平均覆盖度作为映射序列G和种子序列Y间每块的覆盖度，不在映射序列G的坐标范围内的块的覆盖度不计算。

Images