CN112397148B - 序列比对方法、序列校正方法及其装置 - Google Patents

Abstract

本发明提供了一种序列比对方法、序列校正方法及其装置，涉及生物信息技术领域。该序列比对方法主要包括区间比对、碱基比对，配合序列压缩、比对结果优化存储手段，实现了高效、高准确度的比对。该序列校正方法是基于比对结果，通过正序计分、反序回溯，进一步的对低质量区进行校正，实现高效、高准确度的比对。

Description

序列比对方法、序列校正方法及其装置

技术领域

本发明涉及生物信息技术领域，尤其是涉及一种序列比对方法、序列校正方法及装置。

背景技术

随着高通量测序技术的飞速发展，其已成为研究物种进化、诊断疾病成因等生物学问题 的重要手段。

其中Illumina测序平台为主的第二代测序技术以其数据产出量高、成本低廉、碱基准确 率高的特点，在基因组学研究领域得到广泛的应用。基因组组装是指通过测序序列结果重建 完整的基因组序列，组装结果准确性和完整性对于后续基因注释分析、SNP和结构变异检测 等方面产生重要影响。由于存在重复序列和杂合度等复杂情况，采用二代测序技术产生的短 片段序列，难以组装得到完整的基因组序列。以人类基因组为例，现有通过二代测序组装的 人类参考基因组序列中仍存在5％～10％未明确的复杂区域。

第三代测序技术平台，以Pacific Biosciences和Oxford NanoporeTechnologies为代表，具 有长读长的显著特点，能够产生大于10k的长片段序列，并且对GC含量较高的区域测序结 果无偏好性，为基因组研究领域带来了重大的变革。三代测序技术在提升序列长度的同时， 与二代测序相比(错误率0.1％)却带来了碱基序列错误率较高的缺陷(错误率15％～20％)， 测序序列中的错误信息对后续进行数据分析造成了重大的影响。对于三代测序数据，通常采 用Overlap-Layout-Consensus(OLC)算法的策略进行基因组组装，但由于序列错误率高的缺 陷，极大增加了组装的复杂性和计算序列间相互关系的资源消耗。现有针对三代测序数据的 基因组组装软件，或采用先校正测序数据后组装基因组的方法，或采用先组装基因组后对组 装结果进行校正的方法。上述两种方法均涉及序列间的相互比对运算，并基于比对结果校正 测序数据中错误信息。现有比对方法中存在大量的无效比对，该无效比对将连锁影响校正的 准确性，同时比对结果需要占据巨大的存储空间，因此比对及校正步骤在基因组组装中占据 了最长的时间消耗，并对于内存需求量巨大。以组装哺乳动物基因组三代测序数据为例，现 有软件需要消耗数万CPU核时，并且不能处理超大型基因组的组装，严重阻碍了三代测序技 术的广泛应用。因此亟待开发针对性处理长读长测序数据的比对及校正方法，有效较少现有 方法的计算资源消耗，提升运行的效率，充分解锁长读长测序技术在基因组学研究领域的潜 力。所述长读长测序包括，但不限于，单分子实时测序(SMRT)技术和纳米孔测序技术。

发明内容

本发明的目的在于提供一种序列比对方法、序列校正方法及装置，可以较好的减少无效 比对数量，并提高序列校正的准确性。

本发明提供的一种序列对比方法，所述方法包括：S101、获取测序序列，根据第一预设 数据长度在所述测序序列中筛选第一序列，根据第二预设长度在所述测序序列中筛选第二序 列；S102、将所述第一序列与第二序列进行比对，得到第一参考序列、第一比对序列；S103、 计算所述第一比对序列对所述第一参考序列的覆盖度，基于所述覆盖度筛选所述第一参考序 列、第一比对序列，获得第二参考序列、第二比对序列；S104、计算所述第二参考序列与第 二比对序列之间的比对路径，基于编辑距离筛选所述第二参考序列、第二比对序列，获得第 三参考序列、第三比对序列。

进一步，所述方法包括基于二进制数码对序列进行转换和存储，所述序列包括所述测序 序列、所述第一序列、所述第二序列、所述第一参考序列、所述第一比对序列、所述第二参 考序列、所述第二比对序列、所述第三参考序列和所述第三比对序列中的一种或多种；

优选的，所述基于二进制数码对序列进行转换和存储的步骤，包括：按照预设分组将所 述序列划分为多个碱基组合；其中，所述预设分组包括将相邻的四个碱基确定为一个碱基组 合；根据碱基与二进制数码之间预设的转换关系，将所述序列的碱基转换为二进制数码；对 所述序列中少于四个碱基的组合，采用指定的二进制数码对所述少于四个碱基的组合进行补 位扩充，得到满足所述预设分组的二进制序列。

进一步，所述方法包括比对结果的存储优化：纪录比对序列编号、比对方向、比对序列 比对区间起始、比对序列比对区间终止、参考序列编号、参考序列比对区间起始、参考序列 比对区间终止存储比对结果；优选的，纪录比对序列编号差值、参考序列编号差值、参考序 列比对区间长度与比对序列比对区间长度差值存储比对结果；其中所述比对序列包括第一比 对序列、第二比对序列、第三比对序列中的一种或多种，所述参考序列包括第一参考序列、 第二参考序列、第三参考序列中的一种或多种。

进一步，所述存储比对结果的方法包括：按照4字节进行存储；每字节使用7比特位纪 录数值，不足7比特位的使用特定二进制数码进行填充；剩余1比特位使用特定二进制数码 标识比对结果是否终止。

进一步，步骤S103中：所述覆盖度包括窗口覆盖度、整体覆盖度；其中，所述窗口覆盖 度为，按照预设第一碱基数将所述第一参考序列划分为多个窗口，计算所述第一比对序列对 每个窗口的覆盖度；所述整体覆盖度为，所述第一参考序列的平均覆盖度；所述基于所述筛 选所述第一参考序列、第一比对序列的方法为，将满足预设第一覆盖条件的所述第一比对序 列确定为所述第二比对序列，至所述第二比对序列对所述第一参考序列的所述整体覆盖度达 到预设第一覆盖度；其中，所述第一覆盖条件为，所述窗口覆盖度小于预设第二覆盖度，且 小于所述整体覆盖度的第一预设倍数；优选地，根据所述第一比对序列与所述第一参考序列 重叠长度由长至短的顺序，进行所述覆盖度计算及所述筛选。

进一步，步骤S103中还包括对嵌合体序列的过滤，即剔除确定为嵌合体序列的第一参考 序列及其比对信息，其中所述嵌合体序列为，所述第一参考序列中满足预设第二覆盖条件的 序列；其中，所述第二覆盖条件为，至少存在一个窗口，所述窗口的窗口覆盖度小于预设第 三覆盖度且与所述窗口毗邻的指定个数的窗口的窗口覆盖度均大于预设第四覆盖度；或者， 至少存在一个窗口，所述窗口的窗口覆盖度小于预设第五覆盖度；其中，所述第五覆盖度阈 值为与当前窗口毗邻的指定个数的窗口的窗口覆盖度平均值的第二预设倍数。

进一步，步骤S103中还包括对包含序列的过滤，即剔除确定为包含序列的第一参考序列 及其比对信息，其中所述包含序列为，满足下述条件的所述第一参考序列：至少存在一个第 一比对序列，所述第一比对序列与所述第一参考序列比对区间的起始距所述第一参考序列的 起始的碱基数、所述第一比对序列与所述第一参考序列比对区间的终止距所述第一参考序列 的终止的碱基数均在预设第二碱基数范围内。

进一步，步骤S104具体为，在所述第二参考序列与所述第二比对序列比对区间内，执行 以下步骤：S141、使用贪婪算法计算所述第二参考序列与第二比对序列之间的比对路径，选 择最优比对路径并确定其编辑距离；S142、从所述第二比对序列中，选择不大于预设编辑距 离条件的第二比对序列为第三比对序列，其对应的参考序列为第三参考序列；其中，所述预 设编辑距离条件为，所述比对区间内第一参考序列与第二比对序列长度之和的第三预设倍数； S143、循环步骤S141至S142，直至第三比对序列对第三参考序列的覆盖度达到预设第六覆 盖度。

优选的，步骤S141中，贪婪算法计算时约束线范围为，在所述比对路径的编辑距离的约 束下，最大延伸距离减去指定距离值时最小约束线与最大约束线之间。

优选的，步骤S143中，计算所述覆盖度的区间为，根据所述最优比对路径，回溯确定第 三比对序列与第三参考序列的比对区间，查找所述比对区间内的第一个预设第三碱基数完全 匹配的起始位置、最后一个预设第三碱基数完全匹配的终止位置，将所述起始位置、所述终 止位置确定为所述区间的起始位置、终止位置。

本发明提供的一种序列校正方法，所述方法包括：S201、对多个测序片段之间进行比对， 得到参考序列、比对序列；S202、以k-mer为单位，对所述参考序列进行划分，对所述比对 序列进行对应的划分，得到多个k-mer片段；S203、从正序，以k-mer片段为单位，计算每 个位置上所有碱基的得分；S204、从反序，确定最末碱基为最末位置上最高得分碱基，根据 所述最末碱基对应的所述k-mer片段，确定倒数第二位碱基；根据所述倒数第二位碱基对应 的所述k-mer片段，确定倒数第三位碱基；依次类推，获得k-mer得分链；S205、使用k-mer 得分链替换所述参考序列，获得第一校正结果。

优选的，所述序列校正方法中k值为3。

进一步，步骤S203中，根据公式score(p,b)＝max{score(p-1,b)+count_k-mer}-C×0.3计算 每个位置上所有碱基的得分；其中，score(p,b)为所述参考序列p位置处碱基的得分，b为碱 基A、T、G、C或缺失，count_k-mer为所述k-mer得分链在所述比对序列中的出现次数；C为 所述参考序列p位置处所述比对序列对所述参考序列的覆盖度。

进一步，所述校正方法还包括：S206、对所述第一校正结果中的低质量序列，使用步骤 S201比对结果中覆盖所述低质量序列的所述比对序列的比对区间序列作为种子序列，基于所 述种子序列使用有向图算法进行校正；使用校正后的序列替换所述低质量序列，获得第二校 正结果；其中，所述低质量序列为满足以下所有条件的序列，以连续碱基数目不小于预设第 一数值的低质量碱基为起始、区间平均质量值小于预设第二数值的最长序列、区间长度超过 预设第三数值、区间内连续高质量碱基数目不大于预设第一数值。

优选的，所述有向图算法校正中所述种子序列数目不超过预设第五数值，优选为长度排 序位于中间预设第五数值的所述比对序列包含的所述种子序列。

进一步，步骤S206中，所述有向图算法进行校正的具体步骤为：S261、任选一条所述种 子序列作为基序绘制有向图；S262、根据得分矩阵计算所述种子序列中碱基得分，更新有向 图；S263、对有向图进行拓扑排序；S264、根据拓扑排序结果，并确定每个节点的碱基为该 节点上出现概率最高的碱基，获得校正后的序列。

进一步，步骤S262中，对所有的所述种子序列执行以下步骤：a.根据公式score(x)＝max{score(x-1)}-GP对节点得分进行初始化，其中score(x)为节点得分，x为节点索引，GP为空位罚分；设定无前任节点的节点得分为0；b.遍历所述种子序列的每一个碱基，根据公式更新得分矩阵，

其中x为节点索引，y为所述种子序列中碱基索引，GP为空位罚分，M为错配或匹配得 分；c.以出度为0的最高得分节点为起点，选择最高得分的前任节点，进行回溯，以节点索 引为0或碱基索引为0为终点；d.若所述终点的碱基索引大于0，将所述种子序列中位于所 述终点的碱基索引之前的碱基加入有向图；若所述起点的碱基索引小于所述种子序列长度， 将所述种子序列中位于所述起点的碱基索引之后的碱基加入有向图。

进一步，在进行所述有向图算法校正前，按照以下任一条件对所述第一校正结果进行过 滤，即对满足以下任一条件的第一校正结果不进行有向图算法校正：a.所述低质量序列累计 长度超过预设第四数值的所述第一校正结果；b.所述种子序列数目小于预设第一数值的所述 第一校正结果；c.长度超过预设第四数值的所述比对序列占比超过第四预设倍数的所述第一 校正结果。

进一步，所述校正方法还包括：S207、以所述种子序列为所述比对序列，以所述第二校 正结果为所述参考序列，使用前述步骤S201、S202、S203、S204所述方法再次进行校正， 得到第三校正结果。

进一步的，步骤S201中，使用如前述步骤S101、S102、S103、S104所述的比对方法进行比对；其中，所述参考序列包括所述第一参考序列、所述第二参考序列和所述第三参考序列中的一种或多种，所述比对序列包括所述第一比对序列、所述第二比对序列和所述第三比 对序列中的一种或多种。

本发明提供的一种序列对比装置，所述装置包括：序列获取模块，用于获取测序序列， 根据第一预设数据长度在所述测序序列中筛选第一序列，根据第二预设长度在所述测序序列 中筛选第二序列；第一比对模块，用于将所述第一序列与第二序列进行比对，得到第一参考 序列、第一比对序列；第二比对模块，用于计算所述第一比对序列对所述第一参考序列的覆 盖度，基于所述覆盖度筛选所述第一参考序列、第一比对序列，获得第二参考序列、第二比 对序列；第三比对模块，用于计算所述第二参考序列与第二比对序列之间的比对路径，基于 编辑距离筛选所述第二参考序列、第二比对序列，获得第三参考序列、第三比对序列。

本发明提供的一种序列校正装置，所述装置包括：序列比对模块，用于对多个测序片段 之间进行比对，得到参考序列、比对序列；序列划分模块，用于以k-mer为单位，对所述参 考序列、所述比对序列进行划分，得到k-mer片段；得分计算模块，用于从正序，以k-mer片段为单位，计算每个位置上所有碱基的得分；碱基确定模块，用于从反序，确定最末碱基为最末位置上最高得分碱基，根据所述最末碱基对应的所述k-mer片段，确定倒数第二位碱基；根据所述倒数第二位碱基对应的所述k-mer片段，确定倒数第三位碱基；依次类推，获得k-mer得分链；校正模块，用于使用k-mer得分链替换所述参考序列，获得第一校正结果。

本发明提供的序列比对方法及装置，该方法包括：S101、获取测序序列，根据第一预设 数据长度在测序序列中筛选第一序列，根据第二预设长度在测序序列中筛选第二序列；S102、 将第一序列与第二序列进行比对，得到第一参考序列、第一比对序列；S103、计算第一比对 序列对第一参考序列的覆盖度，基于覆盖度筛选第一参考序列、第一比对序列，获得第二参 考序列、第二比对序列；S104、计算第二参考序列与第二比对序列之间的比对路径，基于编 辑距离筛选第二参考序列、第二比对序列，获得第三参考序列、第三比对序列。本发明通过 上述序列比对方式，先区间比对后碱基比对，并对参考序列和比对序列进行多次筛选，可以 减少无效比对数量、提高比对准确性、减少数据存储和运算空间。此外，本发明通过上述方 式提供测序序列之间的比对结果，可以与各种基于比对结果进行校正的校正方法进行组合， 对序列校正方法不进行限制。

本发明提供的序列校正方法及装置，该方法包括：S201、对多个测序片段之间进行比对， 得到参考序列、比对序列；S202、以k-mer为单位，对所述参考序列进行划分，对所述比对 序列进行对应的划分，得到多个k-mer片段；S203、从正序，以k-mer片段为单位，计算每 个位置上所有碱基的得分；S204、从反序，确定最末碱基为最末位置上最高得分碱基，根据 所述最末碱基对应的所述k-mer片段，确定倒数第二位碱基；根据所述倒数第二位碱基对应 的所述k-mer片段，确定倒数第三位碱基；依次类推，获得k-mer得分链；S205、使用k-mer 得分链替换所述参考序列，获得第一校正结果。本发明通过上述序列校正方式先对测序片段 进行比对，然后再基于k-mer片段获得校正结果，可以避免因只统计单碱基比对所导致的校 正结果存在偏好性，从而提高序列校正的准确性，还能够提高序列校正的处理效率。此外， 本发明提供的上述方式是基于序列比对结果进行的校正，不对获得比对结果的比对方法进行 限制，可以与各种序列比对方法进行组合。

可以理解的是，本发明所提供的上述序列校比对方法与上述序列校正方法可以进行组合， 以提高序列校正的准确性。

需要说明的，采用第一、第二、第三等来标记序列，采用S101、S102、S201、S202等对各步骤进行编号，采用第一、第二、第三等对各预设条件或预设值进行标记，均仅为了更清楚的描述。其中，每一种序列以集合的形式描述，意指对该集合中的每一条序列进行操作。 例如，步骤S102中，将所述第一序列与第二序列进行比对，是将第一序列中的所有序列依次 与第二序列中所有序列进行比对。又如，步骤S202中，对所述参考序列进行划分，对所述比 对序列进行对应的划分，是对每一条参考序列均进行k-mer划分，并对与参考序列具有比对 关系的所有比对序列进行对应的划分。

另需要说明的是，本发明提供的一种序列比对方法、序列校正方法及其装置、设备，适 用于对长读长测序结果的处理。对长读长测序结果的来源不做限制，其来源不限于第三代测 序技术。第三代测序技术平台不限于Pacific Biosciences的Sequel、Sequel II平台，Oxford Nanopore Technologies的MinION、PromethION、GridION平台。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方 式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本 发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可 以根据这些附图获得其他的附图。

图1为本发明实施例一提供的序列比对方法流程图；

图2为本发明实施例二提供的二进制形式转换后的序列的示意图；

图3为本发明实施例三提供的二进制形式比对结果优化存储的示意图；

图4为本发明实施例四提供的序列校正方法流程图；

图5为本发明实施例四提供的查找低质量序列的示意图；

图6为本发明实施例七提供的序列比对装置的结构框图；

图7为本发明实施例八提供的序列校正装置的结构框图。

具体实施方式

下面将结合实施例对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例 是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人 员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的，实施例中提及的具体数值，如第一碱基数64、第一覆盖度40、第一预设倍 数1.3、第一数值4、空位罚分2仅为举例说明，不应理解为限制。

目前，超大基因组组装所需的数据量往往达到Tb级别，现有技术中的软件(如Falcon、 Canu)根本无法使用；还考虑到现有针对大型及超大型的基因组的序列校正方案，存在无效 比对数量大且准确性低的问题，本发明实施例提供的一种序列比对、序列校正方法及装置， 可以减少无效比对数量，提高序列校正的准确性。该技术可应用于超大基因组的序列比对和/ 或序列校正的场景中，其中，超大基因组一般是指其基因组大小超过30GB。

实施例一

本实施例提一种序列比对方法。参照如图1所示的序列比对方法流程图，该方法包括如 下步骤：

S101、获取测序序列，根据第一预设数据长度在测序序列中筛选第一序列，根据第二预 设长度在测序序列中筛选第二序列；

S102、将第一序列与第二序列进行比对，得到第一参考序列、第一比对序列；

S103、计算第一比对序列对第一参考序列的覆盖度，基于覆盖度筛选第一参考序列、第 一比对序列，获得第二参考序列、第二比对序列；

S104、计算第二参考序列与第二比对序列之间的比对路径，基于编辑距离筛选第二参考 序列、第二比对序列，获得第三参考序列、第三比对序列。上述第三参考序列和上述第三比 对序列为所述测序序列的比对结果。参考序列与比对序列因序列相似性存在比对关系，具有 比对关系的区间即为比对区间，即体现比对结果；其中，所述参考序列包括第一参考序列、 第二参考序列和第三参考序列中的一种或多种，所述比对序列包括第一比对序列、第二比对 序列和第三比对序列中的一种或多种。

本发明实施例提供的序列比对方法包括：获取测序序列，根据第一预设数据长度在测序 序列中筛选第一序列，根据第二预设长度在测序序列中筛选第二序列；将第一序列与第二序 列进行比对，得到第一参考序列、第一比对序列；计算第一比对序列对第一参考序列的覆盖 度，基于覆盖度筛选第一参考序列、第一比对序列，获得第二参考序列、第二比对序列；计 算第二参考序列与第二比对序列之间的比对路径，基于编辑距离筛选第二参考序列、第二比 对序列，获得第三参考序列、第三比对序列。本实施例通过上述序列比对方式，经两步比对， 并对参考序列和比对序列进行多次筛选，可以减少无效比对数量、提高比对准确性、减少数 据存储和运算空间。此外，本实施例通过上述方式提供测序序列之间的比对结果，可以与各 种基于比对结果进行校正的校正方法进行组合，对序列校正方法不进行限制。

在本实施例步骤S101中，具体实现可参考如下内容。在测序序列reads中选取数据长度 大于第一预设长度(如最长的40×序列长度)的测序序列作为第一序列，并将第一序列输出 至第一序列文件；在测序序列中选取第二预设长度(如高于1k)的测序序列作为第二序列， 并将第二序列输出至第二序列文件。

为便于后续的多任务并行运算，可以将第一序列文件和第二序列文件进行切割，本实施 例对所切割的份数不做限定。

在本实施例步骤S102中，基于上述第一序列和第二序列，提供了一种区间比对方法，将 第一序列中所有reads与第二序列中所有reads进行两两比对，得到第一参考序列、第一比对 序列。此步骤可以使用现有比对软件实现，例如minimap2。现有技术中通常仅过滤掉短序列 (如小于1K)，对过滤后的序列之间全部进行比对，即第二序列与第二序列进行比对，这样 就包括了第二序列中非第一序列之间的比对，此类序列较短(如最长的40×序列长度、1K 之间)且数量多，其比对结果相对冗余。本实施例仅将第一序列与第二序列进行比对，可以 避免冗余比对、提高比对效率，并避免冗余比对对后续校正准确性的负影响。

在本实施例步骤S103中，为了有效减少用于校正的序列的数据集，同时降低因基因组中 重复序列、比对错误导致的多重比对以及比对噪音，基于覆盖度筛选第一参考序列、第一比 对序列，获得比对结果更准确的第二参考序列、第二比对序列。同时，经覆盖度筛选后，获 得的各第二参考序列的覆盖度水平相对一致，有利于后续的校正。

覆盖度包括窗口覆盖度、整体覆盖度；其中，窗口覆盖度为，按照预设第一碱基数(如 64)将第一参考序列划分为多个窗口，计算第一比对序列对每个窗口的覆盖度；整体覆盖度 为，第一参考序列的平均覆盖度。基于覆盖度筛选第一参考序列、第一比对序列的规则为： 将满足预设第一覆盖条件的第一比对序列确定为第二比对序列，且第二比对序列对第一参考 序列的整体覆盖度达到预设第一覆盖度(如40)；其中，第一覆盖条件为，窗口覆盖度小于 预设第二覆盖度(如50)，且小于整体覆盖度的第一预设倍数(如1.3)。为了更高效的进 行覆盖度计算及第一参考序列、第一比对序列筛选，按第一比对序列与第一参考序列比对区 间长度由长至短的顺序操作。具体实现可参考如下内容。

S131：计算第一比对序列与第一参考序列比对区间长度，将第一比对序列按比对区间由 长至短排列；对于比对区间长度相等的第一比对序列，按照其编号由大到小进行排序。

S132：对第一参考序列，按64个碱基长度划分为多个窗口，并将其窗口覆盖度、整体覆 盖度均初始化为0。

S133：对当前第一参考序列的每个窗口，按照S131排列顺序依次查找第一比对序列。如 果该窗口的窗口覆盖度小于预设第二覆盖度50，且低于当前第一参考序列整体覆盖度的第一 预设倍数1.3，则认为该第一比对序列有效，输出该第一比对序列，将与该第一比对序列对应 的窗口的窗口覆盖度加1；反之，表示该第一比对序列冗余，舍弃。直至，当前第一参考序 列整体覆盖度达到预设第一覆盖度40，即表示该第一参考序列有足够的第一比对序列，舍弃 后续所有的第一比对序列。

在步骤S103前，进一步的对第一参考序列中的嵌合体序列和包含序列分别进行过滤。通 过过滤，进一步去除了第一参考序列、第一比对序列中可信度低、冗余的比对结果。

该嵌合体序列为，第一参考序列中满足预设第二覆盖条件的序列；其中，第二覆盖条件 为，至少存在一个窗口，窗口的窗口覆盖度小于预设第三覆盖度(如3)且与窗口毗邻的指 定个数的窗口的窗口覆盖度均大于预设第四覆盖度(如20)；或者，至少存在一个窗口，窗 口的窗口覆盖度小于预设第五覆盖度；其中，第五覆盖度阈值为与当前窗口毗邻的指定个数 的窗口的窗口覆盖度平均值的第二预设倍数(如1/5)。

具体实现可参考如下内容：如果对于任一第一参考序列，其某窗口的窗口覆盖度小于第 三覆盖度3，且附近毗邻的多个窗口的窗口覆盖度大于第四覆盖度20；或者某个窗口的窗口 覆盖度小于附近毗邻的多个窗口的窗口覆盖度平均值的第二预设倍数1/5；则将该条第一参考 序列定义为嵌合体序列，并舍弃该第一参考序列及其比对信息。

该包含序列为，满足下述条件的第一参考序列：至少存在一个第一比对序列，第一比对 序列与第一参考序列比对区间的起始距第一参考序列的起始的碱基数、第一比对序列与第一 参考序列比对区间的终止距第一参考序列的终止的碱基数均在预设第二碱基数(如300)范 围内。

具体实现可参考如下内容：如果对于任一条第一参考序列，其比对区间的起点距第一参 考序列的起点在第二碱基数300bp以内，而比对区间的终点距第一参考序列的终点也在第二 碱基数300bp以内，则定义该第一参考序列为被包含序列，并舍弃该第一参考序列及其比对 信息。

在本实施例步骤S104中，基于步骤S102的区间比对、S103的覆盖度筛选，进一步的在 第二参考序列、第二比对序列的比对区间内进行更为精确的碱基比对。在区间比对结果上进 行碱基比对，而非直接进行碱基比对，其原因在于，可基于更少但更准的数据进行更高效的 比对。通过碱基比对，能够进一步保留可信的、剔除不可信的、去掉冗余的比对信息，进一 步缩减了比对结果、提高了比对精度。

本实施例给出了一种第二参考序列与第二比对序列之间碱基比对的具体实施方式，在第 二参考序列与第二比对序列比对区间内，执行以下步骤：

S141、使用贪婪算法计算第二参考序列与第二比对序列之间的比对路径，选择最优比对 路径并确定其编辑距离；

S142、从第二比对序列中，选择不大于预设编辑距离条件的序列为第三比对序列，其对 应的参考序列为第三参考序列；其中，预设编辑距离条件为，比对区间内第一参考序列与第 二比对序列长度之和的第三预设倍数(如0.4)；

S143、循环步骤S141至S142，直至第三比对序列对第三参考序列的覆盖度达到预设第 六覆盖度(如90)。

进一步的，为了加速步骤S141中贪婪算法的计算速度，将贪婪算法计算时约束线范围限 制为，在所述比对路径的编辑距离的约束下，最大延伸距离减去指定距离值(如150)时最 小约束线与最大约束线之间。其原因在于，针对每一个编辑距离d，都要计算所有的约束线k， 而大部分k线对于查找当前最优的碱基比对结果没有意义。

进一步的，为了进一步准确的计算覆盖度，重新确定计算覆盖度的比对区间为：根据最 优比对路径，回溯确定第三比对序列与第三参考序列的比对区间，查找比对区间内的第一个 预设第三碱基数完全匹配的起始位置、最后一个预设第三碱基数(如8)完全匹配的终止位 置，将起始位置、终止位置确定为区间的起始位置、终止位置。经此步骤可有效解决比对区 间边界(也即端点)错误率较高的问题。

具体实现可参考如下步骤(1)至(4)：

(1)利用Diff贪婪算法，在k∈[k_min，k_max]的约束范围内，计算各路径的编辑距离，选择最优比对路径并确定其编辑距离。其中k_min约束线为在当前编辑距离d的约束下，找到最大延伸距离的最大值减去指定距离值150时最小的k线。k_max为约束线为在当前编辑距离d的约束下，找到最大延伸距离的最大值减去指定距离值150时最大的k线。最优编辑距离为编辑操作次数最少时所对应的编辑距离。

(2)判断最优编辑距离是否大于预设编辑距离条件。预设编辑距离条件为(L1+L2)× 0.4，其中，L1是第二参考序列的数据长度，L2是第二比对序列的数据长度。

若是，表示第二参考序列与第二比对序列为错误比对，舍弃当前最优编辑距离对应的第 二比对序列；若否，表示第二参考序列与第二比对序列为正确比对，保留当前最优编辑距离 对应的第二比对序列，即为第三比对序列，其对应的为第三参考序列。

(3)根据最优编辑距离回溯至比对区间，在比对区间内，从前往后搜索，将第一个第三 参考序列与第三比对序列长度为预设第三碱基数8完全匹配的起点重新设置为比对区间起 点；从后往前搜索，将最后一个第三参考序列与第三比对序列长度为预设第三碱基数8完全 匹配的终点设置为比对区间重点。

(4)计算比对区间的覆盖度，直至第三比对序列对第三参考序列的覆盖度达到预设第六 覆盖度90。

综上，本实施例提供的序列比对方法，通过两次比对、多次筛选，有效的减少无效比对、 剔除错误比对，实现快速、准确、高效的比对。

实施例二

在实际应用中，测序序列的数据量非常庞大，特别是对于大型及超大型的基因组，需占 用较大的硬盘储存空间。将具体的序列转换为二进制形式，可以有效的进行压缩，节约存储 空间。基于此，本实施例所提供的序列对比方法还包括：基于二进制数码对序列进行转换和 存储；序列包括实施例1中的测序序列、第一序列、第二序列、第一参考序列、第一比对序 列、第二参考序列、第二比对序列、第三参考序列和第三比对序列中的一种或多种。

上述步骤包括：(1)按照预设分组将序列划分为多个碱基组合；其中，预设分组策略包 括将相邻的四个碱基确定为一个碱基组合；(2)根据碱基与二进制数码之间预设的转换关系， 将序列的碱基转换为二进制数码；(3)对序列中少于四个碱基的组合，采用指定的二进制数 码对少于四个碱基的组合进行补位扩充，得到满足分组的二进制序列；(4)按照划分的碱基 组合对二进制序列进行存储。

具体的，测序序列包含腺嘌呤(A)、胞嘧啶(C)、鸟嘌呤(G)和胸腺嘧啶(T)四 种形式，每一个碱基都需要1个字节(8个比特位)来存储数据。在本实施例中可以首先将 上述碱基采用二进制形式进行转换，诸如：A＝00，C＝01，G＝10，T＝11，即将每种碱基转换 为占用两个比特位的二进制表现形式。然后将相邻的四个碱基作为一组构成一个字节存储； 因此采用二进制形式压缩后的测序序列约为原始数据大小的四分之一，有效的减少了所需占用的内存空间。

对于测序序列中序列长度是四的倍数的序列，每四个碱基组合均可表示为完整的一个字 节形式。对于序列长度不是四的倍数的序列，最后一组的碱基组合不能占用八个比特位，可 以采用00对最后一组的碱基组合进行补位扩充，使其满足八个比特位即一个字节的形式。

为便于理解，以下给出一种对测序序列进行二进制形式转换的示例：

诸如原始测序序列为：

>6 36535

TTTGCTATTAGT

>252 33213

AAGTATTA AT

采用上述方式将该原始测序序列进行二进制形式转换，编号636535的序列转换后为 111111100111001111001011，编号25233213的序列转换后为000010110011110000110000。具 体的如图2所示，图2中以不同背景填充区别不同存储单元。

同时，可以采用位运算来压缩和解压数据，位运算对CPU消耗较小，可以提高写入和读 取速度。

实施例三

比对结果可以通过记录比对区间的特征进行存储，为了尽可能的缩小存储空间。本实施 例所提供的序列比对方法还包括：纪录比对序列编号、比对方向、比对序列比对区间起始、 比对序列比对区间终止、参考序列编号、参考序列比对区间起始、参考序列比对区间终止存 储比对结果；优选的，纪录比对序列编号差值、参考序列编号差值、参考序列比对区间长度 与比对序列比对区间长度差值存储比对结果；其中比对序列包括实施例一中第一比对序列、 第二比对序列、第三比对序列中的一种或多种，参考序列包括第一参考序列、第二参考序列、 第三参考序列中的一种或多种。此外，按照4字节进行存储；并且，每字节使用7比特位纪 录数值，不足7比特位的使用特定二进制数码进行填充；剩余1比特位使用特定二进制数码 标识比对结果是否终止。

为便于理解，可参照如下示例性描述。上述第一参考序列与第一比对序列的比对结果中 可以包括如下比对信息：第一比对序列编号(第一列)、比对方向(第二列；其表示比对上 的方向，按一个字节8比特位存储，前7位为0，第8位以0或1存储，1代表正向，0表示 反向；)、第一比对序列比对区间起始(第三列)、第一比对序列比对区间终止(第四列)、 第一参考序列编号(第五列)、第一参考序列比对区间起始(第六列)、第一参考序列比对 区间终止(第七列)。

在实际应用中，还可以对第一参考序列与第一比对序列的比对结果做进一步处理：第一 列记录当前第一比对序列编号与上一个比对序列编号的差值，如果该值小于零，第二列数值 加二，第五列记录当前第一参考序列编号与上一个第一参考序列编号的差值，如果该值小于 零，第二列数值加四，第三列和第六列数值不变，第四列记录第一比对序列比对区间的长度， 第七列记录第一比对序列和第一参考序列比对区间长度的差值，如果该值小于零，第二列数 值加八。

由于当前储存数值是按照4个字节(32个比特位)来储存的，而输出结果中的值一般都 很小，只需要2个字节或者更小，高位的比特位都为零，为无效位。因此为了进一步减小存 储空间，可以对第一比对序列中各列的数值做如下压缩处理：按每七个比特位作为一组进行 分割，丢弃前端分组中值为0的组合，对保留的分组除最后一组最高位置为0其他组最高位 均置为1，最后一组不足七比特位的使用0填充存储数值压缩后的字节组合。通过采用上述 方式压缩处理后的第一比对结果的大小通常为原始大小的1/3，明显占用存储空间更小。

为了便于理解，可参照如下具体的列数、数值、二进制的示例。

压缩举例，以第一比对结果中的两条记录为例进行说明。

比对结果：

简化结果为：

二进制存储可参照图3，具体为：

上述二进制存储的序列中双下划线突出数码为填充、标识数码。

实施例四

本实施例提供一种序列校正方法。参照如图4所示的序列校正方法流程图，该方法包括 如下步骤：

S201、对多个测序片段之间进行比对，得到参考序列、比对序列；

S202、以k-mer为单位，对所述参考序列进行划分，对所述比对序列进行对应的划分， 得到多个k-mer片段；

S203、从正序，以k-mer片段为单位，计算每个位置上所有碱基的得分；

S204、从反序，确定最末碱基为最末位置上最高得分碱基，根据最末碱基对应的k-mer 片段，确定倒数第二位碱基；根据倒数第二位碱基对应的k-mer片段，确定倒数第三位碱基； 依次类推，获得k-mer得分链；

S205、使用k-mer得分链替换参考序列，获得第一校正结果。

其中，k-mer，monomeric unit(mer)，相当于nt或者bp，100mer DNA相当于每一条链 有100nt，那么整条链就是100bp。一般长度为m的reads可以分成m-k+1个k-mers。以k-mer 为3举例说明，以k-mer为3对上述待校正序列进行划分，得到多个3-mer片段，比如位点1、 2、3为第一个3-mer，位点2、3、4为第二个3-mer，依次类推。

本发明实施例提供的序列校正方法包括：对多个测序片段之间进行比对，得到参考序列、 比对序列；以k-mer为单位，对所述参考序列进行划分，对所述比对序列进行对应的划分， 得到多个k-mer片段；从正序，以k-mer片段为单位，计算每个位置上所有碱基的得分；从 反序，确定最末碱基为最末位置上最高得分碱基，根据所述最末碱基对应的所述k-mer片段， 确定倒数第二位碱基；根据所述倒数第二位碱基对应的所述k-mer片段，确定倒数第三位碱 基；依次类推，获得k-mer得分链；使用k-mer得分链替换所述参考序列，获得第一校正结 果。本实施例通过上述序列校正方式先对测序片段进行比对，然后再基于k-mer片段对序列 获得校正结果，可以避免因只统计单碱基比对所导致的校正结果存在偏好性，从而提高序列 校正的准确性，还能够提高序列校正的处理效率。此外，本实施例提供的上述方式是基于序 列比对结果进行的校正，不对获得比对结果的比对方法进行限制，可以与各种序列比对方法 进行组合。

上述步骤S203中，计算每个位置上所有碱基的得分是通过如下计算公式：

score(p,b)＝max{score(p-1,b)+count_k-mer}-C×0.3

其中，score(p,b)为参考序列p位置处碱基的得分，b为碱基A、T、G、C或缺失，count_k-mer为k-mer片段在比对序列中的出现次数；C为参考序列p位置处比对序列对参考序列的覆盖 度。

在一种具体的实现方式中，k值可以为3，即k-mer为3-mer。

进一步的，考虑到长读长测序技术存在准确度分布不均匀的情况。以三代测序结果为例， 某些区间准确度极低(<70％)，该区间的比对结果可靠性差，因此基于比对结果进行校正获 得的校正序列准确度低。通过查找低准确度的区间、并针对该区间进行再次校正，可进一步 提高准确性。因此，序列校正方法还可以包括：

步骤S206、对第一校正结果中的低质量序列，使用步骤S201比对结果中覆盖低质量序 列的比对序列的比对区间序列作为种子序列，基于种子序列使用有向图算法进行校正；使用 校正后的序列替换低质量序列，获得第二校正结果；

其中，低质量序列为，以连续碱基数目不小于预设第一数值(如4)的低质量碱基为起 始、区间平均质量值小于预设第二数值(如40)的最长序列、区间长度超过预设第三数值(如 16)、区间内连续高质量碱基数目不大于预设第一数值(如4)。

在一种具体的实现方式中，定义碱基质量值Q为Q＝count_3-mer÷C×100；将碱基质量值大于预设的第一数值4的碱基确定为低质量碱基，将碱基质量值不大于第一数值4 的碱基确定为高质量碱基；区间平均质量值为该区间内所有碱基质量值的平均值。

可参照如下步骤(a)至(c)，基于低质量碱基和高质量碱基查找第一校正结果中的低 质量序列：

(a)对第一校正序列按反序进行低质量碱基查找，记录第一个低质量碱基位置lq_s。

(b)以lq_s位置为起始位置，按反序继续低质量碱基查找，记录查找的连续低质量碱 基数目lq_c，并判断lq_c与第一数值4的关系。

若lq_c不小于第一数值4，继续按反序查找，若查找到高质量碱基则记录连续的高质量 碱基数目hq_c。若hq_c大于第一数值4，初始化lq_c、hq_c为0，按反序重新查找低质量碱 基，重复步骤a、b。若hq_c不大于第一数值4，初始化hq_c为0，继续按反序查找，连续记录lq_c，执行步骤(c)；

若lq_c小于第一数值4，初始化lq_c为0，按反序重新查找低质量碱基，重复步骤a、b。

(c)计算区间平均质量值，当区间平均质量值为小于等于第二数值40的最大值时，判 断区间长度lq_l是否不小于第三数值16。

若大于等于第三数值，记录反序查找的最后一位碱基位置lq_e为终止位置，标记lq_s位 置与lq_e位置之间的序列为低质量序列。同时，初始化低质量区间的长度lq_l、低质量碱基 数目lq_c及高质量碱基数目hq_c为0，按反序重新查找低质量碱基，重复步骤a、b、c，直 至完成所有参考序列的查找。

若小于第三数值，初始化低质量区间的长度lq_l、低质量碱基数目lq_c及高质量碱基数 目hq_c为0，按反序重新查找低质量碱基，重复步骤a、b、c，直至完成所有参考序列的查 找。

为便于理解，可参照如下示例，对第一校正结果查找低质量序列，具体可参照图5：

其中小写部分为低质量碱基，大写部分为高质量碱基，划线部分为找到的低质量区间序 列。

在步骤S206中，有向图算法进行校正的具体步骤为：

S261、任选一条种子序列作为基序绘制有向图。

在本实施例中，为了减少无效数据的比对，可以在进行有向图算法校正前，按照以下任 一条件对第一校正结果进行过滤，即对满足以下任一条件的第一校正结果不进行有向图算法 校正：a.低质量序列累计长度超过预设第四数值(如10000)的第一校正结果；b.种子序列 数目小于预设第一数值(如4)的第一校正结果；c.长度超过预设第四数值(如10000)的 比对序列占比超过第四预设倍数(如1/3)的第一校正结果。过滤掉的数据以低可信度、超长 为特征，可以解决有向图算法中因超长序列造成计算速度极慢的瓶颈问题。

具体的，遍历每一个低质量序列，利用该低质量序列的起点位置和终点位置，查找出所 有覆盖该低质量序列的比对序列。如果查找到的比对序列大于60条，则停止查找低质量区间 序列的比对序列。

进一步的，将查找到的比对序列按序列长度从长到短排序，选取种子序列。在一种较佳 实施例中，种子序列数目不超过预设第五数值(如30)，优选为长度排序位于中间预设第五 数值(如30)的比对序列包含的种子序列。若查找到的比对序列不足30条，则将所有的比 对序列均作为种子序列S。

在种子序列S中任选一条作为基序绘制有向图G；其中，有向图的节点Node为序列中 的碱基，有向图的有向边Edge为任意两个连续碱基之间的连线，有向图的起始节点为序列的 起始碱基。

S262、根据得分矩阵计算种子序列中碱基得分，更新有向图，即将全部种子序列均加入 有向图。

在一种具体的实现方式中，按以下步骤(1)至(4)执行：

(1)根据公式score(x)＝max{score(x-1)}-GP对节点得分进行初始化，其中score(x) 为节点得分，x为节点索引，GP为空位罚分(如2)；设定无前任节点的节点得分为0。

(2)本实施例可以遍历有向图中的节点，基于每一个节点，遍历该节点的所有前任节点。 对于任意一个前任节点，遍历S中的每一个碱基，并按照如下计算公式更新矩阵得分：

其中，x为节点索引，y为种子序列中碱基索引，GP为空位罚分(如2)，M为错配得 分(如-2)或匹配得分(如1)。

(3)以出度为0的最高得分节点为起点，选择最高得分的前任节点，进行回溯，以节点 索引为0或碱基索引为0为终点。

具体的，遍历有向图中的节点，如果节点的出度为0，则定义该节点为终节点，记录该 终节点的得分和索引。找到所有的终结点后，查找出得分最大的终节点，将其确定为回溯起 点，并记录该起始的节点索引和对应的碱基索引。从起点开始循环往回追溯，查找任意节点 x的得分最高的前任节点和前任碱基，直到节点索引为0或者序列索引为0终止，将其确定 为回溯终点，记录终点的节点索引，以及对应的碱基索引。

(4)若终点的碱基索引大于0，将种子序列中位于终点的碱基索引之前的碱基加入有向 图；若起点的碱基索引小于种子序列长度，将种子序列中位于起点的碱基索引之后的碱基加 入有向图。

具体加入有向图的操作为：如果当前碱基在当前节点下为匹配，则跳过该碱基；如果为 错配，则将当前碱基加入到当前节点包含的碱基列表里；如果为插入，则以新节点和新边加 入到有向图中。

S263、对有向图进行拓扑排序；对于更新后的有向图还可以进行拓扑排序，使得更新后 的有向图中任意一对顶点u和v，若满足边(u,v)∈E(G)，则u在排序后的有向图中出现在v 之前。

S264、根据拓扑排序结果，并确定每个节点的碱基为该节点上出现概率最高的碱基，获 得校正后的序列。

在一种具体的实现方式中，按照拓扑排序结果，在每个节点上均选取出现次数最多的碱 基，由此确定校正后的低质量序列，替换第一校正结果的对应序列，获得第二校正结果。

为了进一步提高准确性，序列校正方法还可以包括：

S207、以种子序列为比对序列，以第二校正结果为参考序列，使用前述步骤S201至S205 的序列校正方法，再次进行校正，得到第三校正结果。具体细节不再赘述。

综上，上述实施例中所提供的序列校正方法，可以实现序列高效、高准确度的校正。

实施例五

在实施例4步骤S202中，进行k-mer划分时，k值越大，k-mer片段的特异性越好，校正越准确。但是，当k值越大，完全一致k-mer的概率越小。当测序准确率为p时，出现完 全一致k-mer的概率是p^k。

通过前述score(p,b)的计算公式可知，得分与完全一致k-mer数目相关，若数目过低，即 计入校正的序列过少，会严重影响校对效果。以现阶段三代测序平均准确率为85％为例，2-mer 完全一致的概率是72.25％，3-mer完全一致的概率是61.41％，4-mer完全一致的概率是52.20％。 显而易见的，4-mer及大于4-mer时，概率过低。为了进一步说明k-mer划分的技术效果，通 过示例进行演算，具体如下。

表1为示例数据，比对序列片段包括有Read1～Read7，设定位置1为起始位点，位置5 为结束位点，进行处理。

表1数据信息

示例1：

设定k-mer为3，按照本实施例方法进行校正，通过回溯确定的得分链从位置1至5依次 为GGAGT。

采用k-mer为2进行划分，参照本实施例方法进行校正，通过回溯确定的校正链从位置1 至5依次为GGCGT。

两校正结果相比，位置3的校正结果不同。如前，k值越大，k-mer的特异性越好，校正 结果越准确。因此，3-mer是最佳选择。

示例2：

简单的按单碱基的比对次数进行校正，每个位置选择比对次数最高的作为校正后序列， 各位置单碱基的比对次数统计如表2。

表2单碱基比对次数

示例2的结果为GGCGT，该结果与示例1中k＝2结果一致。可见只统计单碱基比对次数无法进行准确的校正。

在实施例4步骤S204中，从反序，确定最末碱基为最末位置上最高得分碱基，根据所述 最末碱基对应的所述k-mer片段，确定倒数第二位碱基；根据所述倒数第二位碱基对应的所 述k-mer片段，确定倒数第三位碱基；依次类推，获得k-mer得分链。

为了更清楚的说明，通过示例3进行演算，具体如下。

表3为示例数据，比对序列片段包括有Read1～Read5，设定位置1为起始位置，位置7 为结束位置，使用本实施例4的方法进行校正。其中，k-mer设定为3-mer，各位置有效覆盖度均是5。

表3数据信息

使用本实施例的方法，获得的k-mer得分链从位置1至7依次为ACAGACC。

若不进行反向回溯，仅考虑单个碱基，而没有考虑相邻碱基之间的连接关系，即通过最 高得分进行选择。位置3，最高得分碱基为A；位置4，最高得分碱基为G；位置5，最高得分碱基为T；位置6，最高得分碱基为G；位置7，最高得分碱基为C。最后确定的得分链位 置3至位置7为AGTGC，与回溯后的结果不一致，体现了回溯的效果。

由上述示例可见，采用本实施例4方法进行校正时，考虑前向序列与当前碱基的组合情 况，即考虑最可能存在的连续小片段作为校正结果，有效地避免了只统计单碱基比对次数对 于校正结果的偏好性影响。

实施例六

本实施例提供一种序列校正方法，使用实施例一所述的方法对测序序列进行比对；基于 第三参考序列、第三比对序列，使用实施例四所述的方法进行序列的校正。为描述的方便和 简洁，本实施例所描述的序列比对方法和序列校正方法的具体工作过程，可以参考前述方法 实施例一至五中的对应过程，在此不再赘述。

本发明实施例还提供了一种实际应用场景的示例，参照如下内容：

使用PacBio Sequal平台、NanoPore PromethION平台标准流程对一个人类的样本进行测 序，调取20号染色体的所有测序数据作为测试数据，使用该染色体公开的二代数据作为标准 数据，用以评价本发明比对方法、校正方法的技术效果，如表4所示。

表4

其中，校正后得率是指经校正后读长超过5kb的数据占原始数据的比例，本发明、Canu 及Falcon的“reads平均长度”、“最长读长”均针对校正后数据，“≥99％覆盖占比”是指 超过99％的长度能与标准数据相匹配的reads占比，“≥97％相似性占比”是指与标准数据相 匹配的相似性超过97％的reads占比，“≥99％覆盖的平均相似性”是指超过99％的长度能与 标准数据相匹配的reads与标准数据的平均相似性。

从前述20号染色体Nanopore测序结果中随机抽取34条第一序列，共计1588683bp。在 比对时，分别使用如实施例1所述的方法(优化方法)，实施例1中步骤S141中不限制约束线范围的方法(非优化方法)；在校正时，均使用实施例4所述的方法。结果显示，优化方 法运行时间为38.2s，非优化方法运行时间为813.6s。

实施例七

基于上述实施例提供的序列比对方法，本实施例提供一种序列比对装置，参照如图6所 示的序列比对装置的结构框图，该装置包括：

序列获取模块301，用于获取测序序列，根据第一预设数据长度在测序序列中筛选第一 序列，根据第二预设长度在测序序列中筛选第二序列；

第一比对模块302，用于将第一序列与第二序列进行比对，得到第一参考序列、第一比 对序列；

第二比对模块303，用于计算第一比对序列对第一参考序列的覆盖度，基于覆盖度筛选 第一参考序列、第一比对序列，获得第二参考序列、第二比对序列；

第三比对模块304，用于计算第二参考序列与第二比对序列之间的比对路径，基于编辑 距离筛选第二参考序列、第二比对序列，获得第三参考序列、第三比对序列。

在一种实施例中，上述序列比对装置还包括序列存储模块，该序列存储模块用于：基于 二进制数码对序列进行转换和存储；序列包括测序序列、第一序列、第二序列、第一参考序 列、第一比对序列、第二参考序列、第二比对序列、第三参考序列和第三比对序列中的一种 或多种；优选的，基于二进制数码对序列进行转换和存储的步骤，包括：按照预设分组将序 列划分为多个碱基组合；其中，预设分组包括将相邻的四个碱基确定为一个碱基组合；根据 碱基与二进制数码之间预设的转换关系，将序列的碱基转换为二进制数码；对序列中少于四 个碱基的组合，采用指定的二进制数码对少于四个碱基的组合进行补位扩充，得到满足预设 分组的二进制序列；按照划分的碱基组合对二进制序列进行存储。

在一种实施例中，上述序列比对装置还包括比对存储优化模块，该比对存储优化模块用 于：纪录比对序列编号、比对方向、比对序列比对区间起始、比对序列比对区间终止、参考 序列编号、参考序列比对区间起始、参考序列比对区间终止存储比对结果；优选的，纪录比 对序列编号差值、参考序列编号差值、参考序列比对区间长度与比对序列比对区间长度差值 存储比对结果；其中比对序列包括第一比对序列、第二比对序列、第三比对序列中的一种或 多种，参考序列包括第一参考序列、第二参考序列、第三参考序列中的一种或多种。

在一种实施例中，上述比对存储优化模块还用于：按照4字节进行存储；并且，每字节 使用7比特位纪录数值，不足7比特位的使用特定二进制数码进行填充；剩余1比特位使用 特定二进制数码标识比对结果是否终止。

在一种实施例中，上述第二比对模块303中，覆盖度包括窗口覆盖度、整体覆盖度；其 中，窗口覆盖度为，按照预设第一碱基数将第一参考序列划分为多个窗口，计算第一比对序 列对每个窗口的覆盖度；整体覆盖度为，第一参考序列的平均覆盖度；筛选的方法为，将满 足预设第一覆盖条件的第一比对序列确定为第二比对序列，且第二比对序列对第一参考序列 的整体覆盖度达到预设第一覆盖度；其中，第一覆盖条件为，窗口覆盖度小于预设第二覆盖 度，且小于整体覆盖度的第一预设倍数；优选地，按照第一比对序列与第一参考序列重叠长 度由长至短的顺序，进行覆盖度计算及筛选。

在一种实施例中，上述序列比装置还包括对嵌合体序列的过滤模块，其中嵌合体序列为， 第一参考序列中满足预设第二覆盖条件的序列；其中，第二覆盖条件为，至少存在一个窗口， 窗口的窗口覆盖度小于预设第三覆盖度且与窗口毗邻的指定个数的窗口的窗口覆盖度均大于 预设第四覆盖度；或者，至少存在一个窗口，窗口的窗口覆盖度小于预设第五覆盖度；其中， 第五覆盖度阈值为与当前窗口毗邻的指定个数的窗口的窗口覆盖度平均值的第二预设倍数。

在一种实施例中，上述序列比装置还包括对包含序列的过滤模块，其中包含序列为，满 足下述条件的第一参考序列：至少存在一个第一比对序列，第一比对序列与第一参考序列比 对区间的起始距第一参考序列的起始的碱基数、第一比对序列与第一参考序列比对区间的终 止距第一参考序列的终止的碱基数均在预设第二碱基数范围内。

在一种实施例中，上述第三比对模块304中，在第二参考序列与第二比对序列比对区间 内，执行以下步骤：S141、使用贪婪算法计算第二参考序列与第二比对序列之间的比对路径， 选择最优比对路径并确定其编辑距离；S142、从第二比对序列中，选择不大于预设编辑距离 条件的第二比对序列为第三比对序列，其对应的参考序列为第三参考序列；其中，预设编辑 距离条件为，比对区间内第一参考序列与第二比对序列长度之和的第三预设倍数；S143、循 环步骤S141至S142，直至第三比对序列对第三参考序列的覆盖度达到预设第六覆盖度；

优选的，步骤S141中，贪婪算法计算时约束线范围为，在比对路径的编辑距离的约束下， 最大延伸距离减去指定距离值时最小约束线与最大约束线之间；

优选的，步骤S143中，计算覆盖度的区间为，根据最优比对路径，回溯确定第三比对序 列与第三参考序列的比对区间，查找比对区间内的第一个预设第三碱基数完全匹配的起始位 置、最后一个预设第三碱基数完全匹配的终止位置，将起始位置、终止位置确定为区间的起 始位置、终止位置。

实施例八

基于上述实施例提供的序列校正方法，本实施例提供一种序列序列校正装置，参照如图 7所示的序列校正装置的结构框图，该装置包括：

序列比对模块401，用于对多个测序片段之间进行比对，得到参考序列、比对序列；

序列划分模块402，用于以k-mer为单位，对参考序列进行划分，对比对序列进行对应的 划分，得到多个k-mer片段；

得分计算模块403，用于从正序，以k-mer片段为单位，计算每个位置上所有碱基的得分；

碱基确定模块404，用于从反序，确定最末碱基为最末位置上最高得分碱基，根据最末 碱基对应的k-mer片段，确定倒数第二位碱基；根据倒数第二位碱基对应的k-mer片段，确 定倒数第三位碱基；依次类推，获得k-mer得分链；

校正模块405，用于使用k-mer得分链替换参考序列，获得第一校正结果。

在一种实施例中，上述得分计算模块403中，根据公式 score(p,b)＝max{score(p-1,b)+count_k-mer}-C×0.3计算每个位置上所有碱基的得分；

其中，score(p,b)为参考序列p位置处碱基的得分，b为碱基A、T、G、C或缺失，count_k-mer为k-mer片段在比对序列中的出现次数；C为参考序列p位置处比对序列对参考序列的覆盖 度；优选的，k值为3。

在一种实施例中，上述序列校正装置还包括：低质量序列校正模块，用于对第一校正结 果中的低质量序列，使用序列比对模块401比对结果中覆盖低质量序列的比对序列的比对区 间序列作为种子序列，基于种子序列使用有向图算法进行校正；使用校正后的序列替换低质 量序列，获得第二校正结果；其中，低质量序列为满足以下所有条件的序列，以连续碱基数 目不小于预设第一数值的低质量碱基为起始、区间平均质量值小于预设第二数值的最长序列、 区间长度超过预设第三数值、区间内连续高质量碱基数目不大于预设第一数值；优选的，有 向图算法校正中种子序列数目不超过预设第五数值，优选为长度排序位于中间预设第五数值 的比对序列包含的种子序列。

在一种实施例中，上述低质量序列校正模块中还用于：S261、任选一条种子序列作为基 序绘制有向图；S262、根据得分矩阵计算种子序列中碱基得分，更新有向图；S263、对有向 图进行拓扑排序；S264、根据拓扑排序结果，并确定每个节点的碱基为该节点上出现概率最 高的碱基，获得校正后的序列。

在一种实施例中，上述步骤S262中，对所有的种子序列执行以下步骤：a.根据公式score (x)＝max{score(x-1)}-GP对节点得分进行初始化，其中score(x)为节点得分，x为节点索 引，GP为空位罚分；设定无前任节点的节点得分为0；

b.遍历种子序列的每一个碱基，根据公式更新得分矩阵，

其中x为节点索引，y为种子序列中碱基索引，GP为空位罚分，M为错配或匹配得分；

c.以出度为0的最高得分节点为起点，选择最高得分的前任节点，进行回溯，以节点索 引为0或碱基索引为0为终点；d.若终点的碱基索引大于0，将种子序列中位于终点的碱基 索引之前的碱基加入有向图；若起点的碱基索引小于种子序列长度，将种子序列中位于起点 的碱基索引之后的碱基加入有向图。

在一种实施例中，在进行有向图算法校正前，按照以下任一条件对第一校正结果进行过 滤：a.低质量序列累计长度超过预设第四数值的第一校正结果；b.种子序列数目小于预设第 一数值的第一校正结果；c.长度超过预设第四数值的比对序列占比超过第四预设倍数的第一 校正结果。

在一种实施例中，上述序列校正装置还包括：再次校正模块，用于以种子序列为比对序 列，以第二校正结果为参考序列，使用如上任意一项的校正方法，再次进行校正，得到第三 校正结果。

在一种实施例中，上述序列比对模块401中，使用如上任意一项的比对方法进行比对； 其中，参考序列包括第一参考序列、第二参考序列和第三参考序列中的一种或多种，比对序 列包括第一比对序列、第二比对序列和第三比对序列中的一种或多种。

除非另外具体说明，否则在这些实施例中阐述的模块和步骤的相对步骤并不限制本发明 的范围。在这些实施例示出和描述的所有示例中，任何具体值应被解释为仅仅是示例性的， 而不是作为限制，因此，示例性实施例的其他示例可以具有不同的值。

附图中的流程图和框图显示了根据本发明的多个实施例的方法、装置和计算机程序产品 的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个 模块、程序段或代码的一部分，模块、程序段或代码的一部分包含一个或多个用于实现规定 的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也 可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行， 它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程 图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专 用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参 照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以 对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替 换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (15)

1.一种序列对比方法，其特征在于，所述方法包括：

S101、获取测序序列，根据第一预设数据长度在所述测序序列中筛选第一序列，在所述测序序列中选取数据长度大于第一预设长度的测序序列作为第一序列；根据第二预设长度在所述测序序列中筛选第二序列；

S102、将所述第一序列与第二序列进行比对，得到第一参考序列、第一比对序列；

S103、计算所述第一比对序列对所述第一参考序列的覆盖度，基于所述覆盖度筛选所述第一参考序列、第一比对序列，获得第二参考序列、第二比对序列；

S104、计算所述第二参考序列与第二比对序列之间的比对路径，基于编辑距离筛选所述第二参考序列、第二比对序列，获得第三参考序列、第三比对序列；

步骤S103中，所述覆盖度包括窗口覆盖度、整体覆盖度；其中，所述窗口覆盖度为，按照预设第一碱基数将所述第一参考序列划分为多个窗口，计算所述第一比对序列对每个窗口的覆盖度；所述整体覆盖度为，所述第一参考序列的平均覆盖度；

所述筛选的方法为，将满足预设第一覆盖条件的所述第一比对序列确定为所述第二比对序列，且所述第二比对序列对所述第一参考序列的所述整体覆盖度达到预设第一覆盖度；其中，所述第一覆盖条件为，所述窗口覆盖度小于预设第二覆盖度，且小于所述整体覆盖度的第一预设倍数；

按照所述第一比对序列与所述第一参考序列重叠长度由长至短的顺序，进行所述覆盖度计算及所述筛选；

步骤S104中，在所述第二参考序列与所述第二比对序列比对区间内，执行以下步骤：

S141、使用贪婪算法计算所述第二参考序列与第二比对序列之间的比对路径，选择最优比对路径并确定其编辑距离；其中，最优比对路径的最优编辑距离为：编辑操作次数最少时所对应的编辑距离；

S142、从所述第二比对序列中，选择不大于预设编辑距离条件的第二比对序列为第三比对序列，其对应的参考序列为第三参考序列；其中，所述预设编辑距离条件为，所述比对区间内第一参考序列与第二比对序列长度之和的第三预设倍数；

S143、循环步骤S141至S142，直至第三比对序列对第三参考序列的覆盖度达到预设第六覆盖度；

步骤S141中，贪婪算法计算时约束线范围为，在所述比对路径的编辑距离的约束下，最大延伸距离减去指定距离值时最小约束线与最大约束线之间；

步骤S143中，计算所述覆盖度的区间为，

根据所述最优比对路径，回溯确定第三比对序列与第三参考序列的比对区间，查找所述比对区间内的第一个预设第三碱基数完全匹配的起始位置、最后一个预设第三碱基数完全匹配的终止位置，将所述起始位置、所述终止位置确定为所述区间的起始位置、终止位置。

2.根据权利要求1所述的方法，其特征在于，所述方法包括：

基于二进制数码对序列进行转换和存储；所述序列包括所述测序序列、所述第一序列、所述第二序列、所述第一参考序列、所述第一比对序列、所述第二参考序列、所述第二比对序列、所述第三参考序列和所述第三比对序列中的一种或多种；

所述基于二进制数码对序列进行转换和存储的步骤，包括：

按照预设分组将所述序列划分为多个碱基组合；其中，所述预设分组包括将相邻的四个碱基确定为一个碱基组合；

根据碱基与二进制数码之间预设的转换关系，将所述序列的碱基转换为二进制数码；

对所述序列中少于四个碱基的组合，采用指定的二进制数码对所述少于四个碱基的组合进行补位扩充，得到满足所述预设分组的二进制序列。

3.根据权利要求1所述的方法，其特征在于，所述方法包括比对结果的存储优化：

纪录比对序列编号、比对方向、比对序列比对区间起始、比对序列比对区间终止、参考序列编号、参考序列比对区间起始、参考序列比对区间终止存储比对结果；

纪录比对序列编号差值、参考序列编号差值、参考序列比对区间长度与比对序列比对区间长度差值存储比对结果；

其中所述比对序列包括第一比对序列、第二比对序列、第三比对序列中的一种或多种，所述参考序列包括第一参考序列、第二参考序列、第三参考序列中的一种或多种。

4.根据权利要求3所述的方法，其特征在于，所述存储比对结果的方法包括：

按照4字节进行存储；

并且，每字节使用7比特位纪录数值，不足7比特位的使用特定二进制数码进行填充；剩余1比特位使用特定二进制数码标识比对结果是否终止。

5.根据权利要求1所述的方法，其特征在于，在步骤S103之前，所述方法还包括对嵌合体序列的过滤，其中所述嵌合体序列为，

所述第一参考序列中满足预设第二覆盖条件的序列；

其中，所述第二覆盖条件为，至少存在一个窗口，所述窗口的窗口覆盖度小于预设第三覆盖度且与所述窗口毗邻的指定个数的窗口的窗口覆盖度均大于预设第四覆盖度；或者，至少存在一个窗口，所述窗口的窗口覆盖度小于预设第五覆盖度；其中，所述第五覆盖度阈值为与当前窗口毗邻的指定个数的窗口的窗口覆盖度平均值的第二预设倍数。

6.根据权利要求1所述的方法，其特征在于，在步骤S103之前，所述方法还包括对包含序列的过滤，其中所述包含序列为，

满足下述条件的所述第一参考序列：至少存在一个第一比对序列，所述第一比对序列与所述第一参考序列比对区间的起始距所述第一参考序列的起始的碱基数、所述第一比对序列与所述第一参考序列比对区间的终止距所述第一参考序列的终止的碱基数均在预设第二碱基数范围内。

7.一种序列校正方法，其特征在于，所述方法包括：

S201、对多个测序片段之间进行比对，得到参考序列、比对序列；

S202、以k-mer为单位，对所述参考序列进行划分，对所述比对序列进行对应的划分，得到多个k-mer片段；

S203、从正序，以k-mer片段为单位，计算每个位置上所有碱基的得分；

S204、从反序，确定最末碱基为最末位置上最高得分碱基，根据所述最末碱基对应的所述k-mer片段，确定倒数第二位碱基；根据所述倒数第二位碱基对应的所述k-mer片段，确定倒数第三位碱基；依次类推，获得k-mer得分链；

S205、使用k-mer得分链替换所述参考序列，获得第一校正结果；

所述序列校正方法中k值为3；

步骤S201中，使用如权利要求1至6任意一项所述的序列比对方法进行比对。

8.根据权利要求7所述的方法，其特征在于，步骤S203中，

根据公式计算每个位置上所有碱基的得分；

其中，为所述参考序列p位置处碱基的得分，b为碱基A、T、G、C或缺失，count_k-mer为所述k-mer片段在所述比对序列中的出现次数；C为所述参考序列p位置处所述比对序列对所述参考序列的覆盖度。

9.根据权利要求7所述的方法，其特征在于，所述方法还包括：

S206、对所述第一校正结果中的低质量序列，使用步骤S201比对结果中覆盖所述低质量序列的所述比对序列的比对区间序列作为种子序列，基于所述种子序列使用有向图算法进行校正；使用校正后的序列替换所述低质量序列，获得第二校正结果；

其中，所述低质量序列为满足以下所有条件的序列，以连续碱基数目不小于预设第一数值的低质量碱基为起始、区间平均质量值小于预设第二数值的最长序列、区间长度超过预设第三数值、区间内连续高质量碱基数目不大于预设第一数值；

有向图算法校正中所述种子序列数目为长度排序位于中间预设第五数值的所述比对序列包含的所述种子序列。

10.根据权利要求9所述的方法，其特征在于，步骤S206中，所述有向图算法进行校正的具体步骤为：

S261、任选一条所述种子序列作为基序绘制有向图；

S262、根据得分矩阵计算所述种子序列中碱基得分，更新有向图；

S263、对有向图进行拓扑排序；

S264、根据拓扑排序结果，并确定每个节点的碱基为该节点上出现概率最高的碱基，获得校正后的序列。

11.根据权利要求10所述的方法，其特征在于，步骤S262中，对所有的所述种子序列执行以下步骤：

a. 根据公式score（x）=max{score(x-1)}-GP对节点得分进行初始化，其中score（x）为节点得分，x为节点索引，GP为空位罚分；设定无前任节点的节点得分为0；

b. 遍历所述种子序列的每一个碱基，根据公式更新如下得分矩阵：

其中x为节点索引，y为所述种子序列中碱基索引，GP为空位罚分，M为错配或匹配得分；

c. 以出度为0的最高得分节点为起点，选择最高得分的前任节点，进行回溯，以节点索引为0或碱基索引为0为终点；

d. 若所述终点的碱基索引大于0，将所述种子序列中位于所述终点的碱基索引之前的碱基加入有向图；若所述起点的碱基索引小于所述种子序列长度，将所述种子序列中位于所述起点的碱基索引之后的碱基加入有向图。

12.根据权利要求9所述的方法，其特征在于，在进行有向图算法校正前，按照以下任一条件对所述第一校正结果进行过滤：

a. 所述低质量序列累计长度超过预设第四数值的所述第一校正结果；

b. 所述种子序列数目小于预设第一数值的所述第一校正结果；

c. 长度超过预设第四数值的所述比对序列占比超过第四预设倍数的所述第一校正结果。

13.根据权利要求9所述的方法，所述方法还包括：

S207、以所述种子序列为所述比对序列，以所述第二校正结果为所述参考序列，使用如权利要求7至8任意一项所述的校正方法，再次进行校正，得到第三校正结果。

14.一种序列对比装置，其特征在于，所述装置包括：

序列获取模块，用于获取测序序列，根据第一预设数据长度在所述测序序列中筛选第一序列，在所述测序序列中选取数据长度大于第一预设长度的测序序列作为第一序列；根据第二预设长度在所述测序序列中筛选第二序列；

第一比对模块，用于将所述第一序列与第二序列进行比对，得到第一参考序列、第一比对序列；

第二比对模块，用于计算所述第一比对序列对所述第一参考序列的覆盖度，基于所述覆盖度筛选所述第一参考序列、第一比对序列，获得第二参考序列、第二比对序列；

第三比对模块，用于计算所述第二参考序列与第二比对序列之间的比对路径，基于编辑距离筛选所述第二参考序列、第二比对序列，获得第三参考序列、第三比对序列；

第二比对模块中，所述覆盖度包括窗口覆盖度、整体覆盖度；其中，所述窗口覆盖度为，按照预设第一碱基数将所述第一参考序列划分为多个窗口，计算所述第一比对序列对每个窗口的覆盖度；所述整体覆盖度为，所述第一参考序列的平均覆盖度；

第二比对模块还用于，将满足预设第一覆盖条件的所述第一比对序列确定为所述第二比对序列，且所述第二比对序列对所述第一参考序列的所述整体覆盖度达到预设第一覆盖度；其中，所述第一覆盖条件为，所述窗口覆盖度小于预设第二覆盖度，且小于所述整体覆盖度的第一预设倍数；按照所述第一比对序列与所述第一参考序列重叠长度由长至短的顺序，进行所述覆盖度计算及所述筛选；

第三比对模块，在所述第二参考序列与所述第二比对序列比对区间内，还用于执行以下步骤：

S141、使用贪婪算法计算所述第二参考序列与第二比对序列之间的比对路径，选择最优比对路径并确定其编辑距离；其中，最优比对路径的最优编辑距离为：编辑操作次数最少时所对应的编辑距离；

S142、从所述第二比对序列中，选择不大于预设编辑距离条件的第二比对序列为第三比对序列，其对应的参考序列为第三参考序列；其中，所述预设编辑距离条件为，所述比对区间内第一参考序列与第二比对序列长度之和的第三预设倍数；

S143、循环步骤S141至S142，直至第三比对序列对第三参考序列的覆盖度达到预设第六覆盖度；

步骤S141中，贪婪算法计算时约束线范围为，在所述比对路径的编辑距离的约束下，最大延伸距离减去指定距离值时最小约束线与最大约束线之间；

步骤S143中，计算所述覆盖度的区间为，根据所述最优比对路径，回溯确定第三比对序列与第三参考序列的比对区间，查找所述比对区间内的第一个预设第三碱基数完全匹配的起始位置、最后一个预设第三碱基数完全匹配的终止位置，将所述起始位置、所述终止位置确定为所述区间的起始位置、终止位置。

15.一种序列校正装置，其特征在于，所述装置包括：

序列比对模块，用于对多个测序片段之间进行比对，得到参考序列、比对序列；

序列划分模块，用于以k-mer为单位，对所述参考序列进行划分，对所述比对序列进行对应的划分，得到多个k-mer片段；

得分计算模块，用于从正序，以k-mer片段为单位，计算每个位置上所有碱基的得分；

碱基确定模块，用于从反序，确定最末碱基为最末位置上最高得分碱基，根据所述最末碱基对应的所述k-mer片段，确定倒数第二位碱基；根据所述倒数第二位碱基对应的所述k-mer片段，确定倒数第三位碱基；依次类推，获得k-mer得分链；

校正模块，用于使用k-mer得分链替换所述参考序列，获得第一校正结果；

所述序列校正装置中k值为3；

所述序列比对模块，用于使用如权利要求1至6任意一项所述的序列比对方法进行比对。