逆向加权的序列比对种子生成方法、装置、设备和存储器

技术领域

本申请涉及计算机技术领域，特别是涉及一种逆向加权的序列比对种子生成方法、装置、设备和存储器。

背景技术

序列比对将测序仪输出的读段序列与参考基因组比较，找出读段在基因组上最可能的起源位置。序列比对是测序数据分析中基础而重要的环节，比对的结果是其他步骤的前提。基于种子-扩展的序列比对算法可以有效缩小碱基级比对的范围，较少计算量，提高比对速度。然而基因组重复区域的大量重复片段会造成种子定位错误，产生假阳性命中位点，影响序列比对的准确度。

现有的种子生成方法中，最小子种子(minimizer)方法具有时间和空间上的双重优势，应用越来越普遍。该方法将基因组均匀分为若干长度相同的窗口，在每个窗口中，生成连续的k-mer(长度为k个碱基的片段)，通过事先定义的哈希函数计算每个k-mer的哈希值，并选择哈希值最小的k-mer作为每个窗口的最小子种子。然而真核生物等基因组中存在很多重复片段，如着丝点附近的卫星DNA等串联重复序列，这会导致基因组重复区域选取的最小子种子频率高于其他区域。当比对时，这些高频的最小子种子会引导映射算法定位到重复区域，带来很多假阳性位点，并在碱基比对阶段造成大量徒劳计算和空间浪费。

Heng Li在论文《Minimap2:pairwise alignment for nucleotide sequences》中公布了一种方法，统计所有最小子种子的频率，并丢弃频率最高的前0.02％的最小子种子。Chirag Jain等在论文《A Fast Approximate Algorithm for Mapping Long Reads toLarge Reference Databases》中公布了一种方法，丢弃频率最高的前0.001％的最小子种子。这两种方法虽然在一定程度上避免了假阳性位点，但是却破坏了最小子种子算法的完整性，而且也降低了基因组重复区域的比对准确性，会对下游数据分析造成影响，已有相关研究论文表示，上述方案造成了基因组重复区域的同源性检测率降低，高拷贝质粒的组装出错等问题。

发明内容

基于此，有必要针对上述技术问题，提供一种逆向加权的序列比对种子生成方法、装置、设备和存储器。

一种逆向加权的序列比对种子生成方法，所述方法包括：

读取参考序列，并将所述参考序列分割为N个窗口，其中前N-1个窗口均包含w个碱基，第N个窗口的碱基数量不大于w，N、w均为大于1的整数。

根据每个所述窗口，为对应窗口生成连续重叠的k-mer。

对所有窗口中的每条k-mer的出现次数进行统计，得到每条k-mer的出现次数。

根据每条k-mer的出现次数，计算每条k-mer的出现频率。

根据每条k-mer的出现频率，采用预设权重函数，计算每条k-mer的权重值；所述k-mer的权重值随着k-mer的出现频率的增大而减小。

根据每条k-mer的碱基序列信息和出现频率，采用预设哈希函数，计算得到每条k-mer的哈希值。

分别比较每个窗口中所有k-mer的哈希值，生成每个窗口的最小子种子。

在其中一个实施例中，所述根据每条k-mer的碱基序列信息和出现频率采用预设哈希函数，计算每条k-mer的哈希值，步骤中所述预设哈希函数为：

其中，h(k

在其中一个实施例中，分别比较每个窗口中所有k-mer的哈希值，生成每个窗口的最小子种子，包括：

将第一个窗口作为当前窗口；将当前窗口中的第一条k-mer作为当前k-mer；设置当前窗口中k-mer的最小哈希值为1；将当前窗口中当前k-mer的哈希值与最小哈希值进行比较，得到比较结果；当所述比较结果为当前k-mer的哈希值小于最小哈希值时，则将最小哈希值更新为当前k-mer的哈希值，并记录对应的k-mer为当前窗口的备选最小子种子；当所述比较结果为当前k-mer的哈希值大于最小哈希值时，如果当前k-mer不是当前窗口中最后一条k-mer，则处理下一条k-mer；如果当前k-mer是当前窗口中最后一条k-mer，则将当前窗口的备选最小子种子作为当前窗口的最小子种子，将当前窗口更新为下一窗口，继续进行下一轮迭代，直到遍历所有窗口为止，完成所有窗口的最小子种子生成过程。

在其中一个实施例中，对所有窗口中的每条k-mer的出现次数进行统计，得到每条k-mer的出现次数，包括：

设置频率记录表为空，所述频率记录表中的信息包括：k-mer、k-mer的出现次数、k-mer的出现频率；读取一条k-mer，记为k

一种逆向加权的序列比对种子生成装置，所述装置包括：

参考序列读取模块，用于读取参考序列，并将所述参考序列分割为N个窗口，其中前N-1个窗口均包含w个碱基，第N个窗口的碱基数量不大于w，N、w均为大于1的整数。

k-mer出现频率统计模块，用于根据每个所述窗口，为对应窗口生成连续重叠的k-mer；对所有窗口中的每条k-mer的出现次数进行统计，得到每条k-mer的出现次数；根据每条k-mer的出现次数，计算每条k-mer的出现频率。

k-mer权重计算模块，根据每条k-mer的出现频率，采用预设权重函数，计算每条k-mer的权重值；所述k-mer的权重值随着k-mer的出现频率的增大而减小。

k-mer哈希值计算模块，用于根据每条k-mer的碱基序列信息和出现频率，采用预设哈希函数，计算每条k-mer的哈希值。

最小子种子生成模块，用于分别比较每个窗口中所有k-mer的哈希值，生成每个窗口的最小子种子。

在其中一个实施例中，k-mer哈希值计算模块中所述预设哈希函数为：

其中，h(k

在其中一个实施例中，最小子种子生成模块，还用于将第一个窗口作为当前窗口；将当前窗口中的第一条k-mer作为当前k-mer；设置当前窗口中k-mer的最小哈希值为1；将当前窗口中当前k-mer的哈希值与最小哈希值进行比较，得到比较结果；当所述比较结果为当前k-mer的哈希值小于最小哈希值时，则将最小哈希值更新为当前k-mer的哈希值，并记录对应的k-mer为当前窗口的备选最小子种子；当所述比较结果为当前k-mer的哈希值大于最小哈希值时，如果当前k-mer不是当前窗口中最后一条k-mer，则处理下一条k-mer；如果当前k-mer是当前窗口中最后一条k-mer，则将当前窗口的备选最小子种子作为当前窗口的最小子种子，将当前窗口更新为下一窗口，继续进行下一轮迭代，直到遍历所有窗口为止，完成所有窗口的最小子种子生成过程。

在其中一个实施例中，k-mer出现频率统计模块，还用于设置频率记录表为空，所述频率记录表中的信息包括：k-mer、k-mer的出现次数、k-mer的出现频率；读取一条k-mer，记为k

一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现上述方法的步骤。

一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述方法的步骤。

上述逆向加权的序列比对种子生成方法、装置、设备和存储器，所述方法包括：读取参考序列，并将参考序列分割为N个窗口，其中每一个窗口包含w个碱基；根据每个窗口，为对应窗口生成连续重叠的k-mer；对所有窗口中的每条k-mer的出现次数进行统计，得到每条k-mer的出现次数；根据每条k-mer的出现次数，计算每条k-mer的出现频率；根据每条k-mer的出现频率，采用预设权重函数，计算每条k-mer的权重值，k-mer的权重值随着k-mer的出现频率的增大而减小；根据每条k-mer的碱基序列信息和出现频率，采用预设哈希函数，计算每条k-mer的哈希值；分别比较每个窗口中所有k-mer的哈希值，生成每个窗口的最小子种子。本方法避免了通过删除高频种子来降低假阳性位点，并且保证了种子生成算法的完整性，实现了对基因组重复区域逆向加权种子生成，支持了精准比对算法的设计和实现。

附图说明

图1为一个实施例中逆向加权的序列比对种子生成方法的流程示意图；

图2为一个实施例中逆向加权的序列比对种子生成装置的结构框图；

图3为一个实施例中计算机设备的内部结构图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

基因组重复区域存在假阳性定位的原因在于，重复区域的最小子种子频率过高，如果能够避免选择频率过高的k-mer作为最小子种子，则可以解决此问题。因此，本发明提出采用带权重的种子选取方法，给每个k-mer赋予一个权重，权重越高，这个k-mer被选为最小子种子的概率就越大，权重越低，被选为最小子种子的概率就越小。对于频率较高的k-mer，赋予较低的权重，实现逆向加权，这样高频k-mer被选为最小子种子的概率会降低，也就不会带来过多假阳性位点。

在一个实施例中，如图1所示，提供了一种逆向加权的序列比对种子生成方法，该方法包括以下步骤：

步骤100：读取参考序列，并将参考序列分割为N个窗口，其中前N-1个窗口均包含w个碱基，第N个窗口的碱基数量不大于w，N、w均为大于1的整数。

具体的，在基因比对领域，所述参考序列是有多年积累建立的基因碱基序列样版，又称标准基因库，代表着目前已知的基因与基因效果的对应关系，通过将待比对序列与参考序列的对比可以预测待比对序列的基因效果。

步骤102：根据每个窗口，为对应窗口生成连续重叠的k-mer。

具体的，读取一个窗口，为该窗口生成连续重叠的k-mer，若不是第N个窗口，则为w-k+1条k-mer，记为k

步骤104：对所有窗口中的每条k-mer的出现次数进行统计，得到每条k-mer的出现次数。

步骤106：根据每条k-mer的出现次数，计算每条k-mer的出现频率。

具体的，读取频率记录表中的一条k-mer及其出现次数n

步骤108：根据每条k-mer的出现频率，采用预设权重函数，计算每条k-mer的权重值；k-mer的权重值随着k-mer的出现频率的增大而减小。

具体的，读取频率记录表中的一条k-mer及其出现频率f(k

其中，预设权重函数采用如下方法设计：预设权重函数φ的输入为k-mer的出现频率f(k

步骤110：根据每条k-mer的碱基序列信息和出现频率，采用预设哈希函数，计算每条k-mer的哈希值。

具体的，读取频率记录表中一条k-mer k

步骤112：分别比较每个窗口中所有k-mer的哈希值，生成每个窗口的最小子种子。

上述逆向加权的序列比对种子生成方法中，所述方法包括：读取参考序列，并将参考序列分割为N个窗口；根据每个窗口，为对应窗口生成连续重叠的k-mer；对所有窗口中的每条k-mer的出现次数进行统计，得到每条k-mer的出现次数；根据每条k-mer的出现次数，计算每条k-mer的出现频率；根据每条k-mer的出现频率，采用预设权重函数，计算每条k-mer的权重值，k-mer的权重值随着k-mer的出现频率的增大而减小；根据每条k-mer的碱基序列信息和出现频率，采用预设哈希函数，计算每条k-mer的哈希值；分别比较每个窗口中所有k-mer的哈希值，生成每个窗口的最小子种子。本方法避免了通过删除高频种子来降低假阳性位点，并且保证了种子生成算法的完整性，实现了对基因组重复区域逆向加权种子生成，支持了精准比对算法的设计和实现。

在其中一个实施例中，步骤110中预设哈希函数为：

其中，h(k

在其中一个实施例中，步骤112包括：将第一个窗口作为当前窗口；将当前窗口中的第一条k-mer作为当前k-mer；设置当前窗口中k-mer的最小哈希值为1；将当前窗口中当前k-mer的哈希值与最小哈希值进行比较，得到比较结果；当比较结果为当前k-mer的哈希值小于最小哈希值时，则将最小哈希值更新为当前k-mer的哈希值，并记录对应的k-mer为当前窗口的备选最小子种子；当比较结果为当前k-mer的哈希值大于最小哈希值时，如果当前k-mer不是当前窗口中最后一条k-mer，则处理下一条k-mer；如果当前k-mer是当前窗口中最后一条k-mer，则将当前窗口的备选最小子种子作为当前窗口的最小子种子，将当前窗口更新为下一窗口，继续进行下一轮迭代，直到遍历所有窗口为止，完成所有窗口的最小子种子生成过程。

在其中一个实施例中，步骤104包括：设置频率记录表为空，频率记录表中的信息包括：k-mer、k-mer的出现次数、k-mer的出现频率；读取一条k-mer，记为k

表2频率记录表

应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。

在一个实施例中，如图2所示，提供了一种逆向加权的序列比对种子生成装置，包括：参考序列读取模块、k-mer出现频率统计模块、k-mer权重计算模块、k-mer哈希值计算模块和最小子种子生成模块，其中：

参考序列读取模块，用于读取参考序列，并将参考序列分割为N个窗口，其中前N-1个窗口均包含w个碱基，第N个窗口的碱基数量不大于w，N、w均为大于1的整数。

k-mer出现频率统计模块，用于根据每个窗口，为对应窗口生成连续重叠的k-mer；对所有窗口中的每条k-mer的出现次数进行统计，得到每条k-mer的出现次数；根据每条k-mer的出现次数，计算每条k-mer的出现频率。

k-mer权重计算模块，根据每条k-mer的出现频率，采用预设权重函数，计算每条k-mer的权重值；k-mer的权重值随着k-mer的出现频率的增大而减小。

k-mer哈希值计算模块，用于根据每条k-mer的碱基序列信息和出现频率，采用预设哈希函数，计算每条k-mer的哈希值。

最小子种子生成模块，用于分别比较每个窗口中所有k-mer的哈希值，生成每个窗口的最小子种子。

在其中一个实施例中，k-mer哈希值计算模块中预设哈希函数为：

其中，h(k

在其中一个实施例中，最小子种子生成模块，还用于将第一个窗口作为当前窗口；将当前窗口中的第一条k-mer作为当前k-mer；设置当前窗口中k-mer的最小哈希值为1；将当前窗口中当前k-mer的哈希值与最小哈希值进行比较，得到比较结果；当比较结果为当前k-mer的哈希值小于最小哈希值时，则将最小哈希值更新为当前k-mer的哈希值，并记录对应的k-mer为当前窗口的备选最小子种子；当比较结果为当前k-mer的哈希值大于最小哈希值时，如果当前k-mer不是当前窗口中最后一条k-mer，则处理下一条k-mer；如果当前k-mer是当前窗口中最后一条k-mer，则将当前窗口的备选最小子种子作为当前窗口的最小子种子，将当前窗口更新为下一窗口，继续进行下一轮迭代，直到遍历所有窗口为止，完成所有窗口的最小子种子生成过程。

在其中一个实施例中，k-mer出现频率统计模块，还用于设置频率记录表为空，频率记录表中的信息包括：k-mer、k-mer的出现次数、k-mer的出现频率；读取一条k-mer，记为k

关于逆向加权的序列比对种子生成装置的具体限定可以参见上文中对于逆向加权的序列比对种子生成方法的限定，在此不再赘述。上述逆向加权的序列比对种子生成装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。

在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图3所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种逆向加权的序列比对种子生成方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。

本领域技术人员可以理解，图3中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。

在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。

在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。