串联重复序列比对的位置筛选方法

摘 要：串联重复序列是基因组构建的困难片段，由于其重复单元之间的相似性与其拷贝数的不确定性，在序列比对时容易定位到多个候选位置，如何快速而准确地筛选出正确的比对位置是一项挑战。现有方法使用种子（从测序片段中选取的短序列）来定位并扩展候选比对位置，但挑选种子时未考虑串联重复序列特性。因此，提出了一种串联重复序列比对的位置筛选方法，其通过计算稀有kmer（长度为k的子序列）序列的相似性来筛选比对结果。此外，采用合并稀有kmer的策略加速计算，并利用基于编辑距离的模糊查找以提高过滤信息密度。实验结果表明，在模拟数据集上提高比对结果的召回率与准确率的同时，该方法比现有方法快约2倍，且具有良好的并行加速性能。

关键词：串联重复； 单分子实时测序； 序列比对； 种子-扩展法

中图分类号：TP391 文献标志码：A 文章编号：1001-3695（2024）07-033-2160-05

doi：10.19734/j.issn.1001-3695.2023.12.0614

Position filtering method for tandem repeat sequence alignment

Abstract：Tandem repeat sequences are difficult part in genome construction， due to the high similarity between repeated units and the ambiguity in copy numbers，it often result in multiple candidate positions during sequence alignment. The challenge lies in rapidly and accurately filtering out the correct alignment positions. Existing methods address this issue by using seeds（short sequences selected from sequencing fragments） to locate and extend candidate alignment positions，but overlook the distinctive characteristics of tandem repeat sequences when selecting seeds. To tackle the problem， this paper proposed a position filtering method for tandem repeat sequence alignment， which filtered alignment results by calculating the similarity of rare kmer sequences. Additionally， it implemented a strategy of merging rare kmers to expedite computation， coupled with a fuzzy search based on edit distance to enhance filtering information density. Experimental results demonstrate that this approach improves both the recall and accuracy of alignment results on simulated datasets while achieving approximately a 2-fold increase in computational speed compared to existing methods， with notable parallel acceleration effects.

Key words：tandem repeat; single molecule real-time sequencing; sequence alignment; seed-and-extend method

0 引言

随着基因测序技术的高速发展，新兴的三代测序技术，如单分子实时测序技术（single molecule real-time，SMRT）［1］和牛津纳米孔测序技术（Oxford nanopore technology，ONT）［2］，连同二代测序技术（next generation sequencing，NGS）［3］一起助力了越来越多端粒到端粒T2T（telomere to telomere）基因组的诞生。其中，串联重复序列是一种由高度相似的重复单元紧密排列而成的特殊结构，构成了T2T基因组中的许多重要组件，例如着丝粒、端粒、rDNA等［4］，是T2T基因组中最复杂的结构之一。研究表明，真核生物的基因组很大部分由重复区域构成［5］，特别是哺乳动物中重复区域约占整个基因组的25%～50%［6］。然而，重复单元之间的相似性与其拷贝数的不确定性为许多T2T基因组分析问题在计算上带来了挑战［7］，其中一个关键问题就是序列比对。许多生物信息学问题的解决，例如denovo组装［8］、结构变异检测［9］与进化分析［10］等都依赖于序列比对结果的准确性［11］。

在生物信息学中，序列比对是重要的基础性问题之一。基因序列比对的问题定义为利用相关算法，将基因测序得到的reads比对到已知基因序列上，得到其相对位置关系，从而判断序列相似性与相关生物学特征。主流的序列比对算法普遍采用种子-扩展法（seed-and-extend）［12］，该方法通过种子索引来定位比对位置，将比对范围从全基因组缩小至若干候选位置，再对候选位置进行精确的比对，从而减少整体时间消耗。在面对串联重复序列时，大多数通用序列比对算法选取的种子不足以区分重复单元之间的细微差别，导致定位到多个候选比对位置，进而产生大量假阳性比对结果，形成多重比对现象［13］。例如，Minimap2［14］在将7 532条来自串联重复区域的reads比对到参考基因组时产生了多达29 106条比对结果，对筛选出正确的比对位置造成了困难。因此，设计高效的序列比对位置筛选方法，以过滤串联重复序列的假阳性比对结果，是基因序列比对面临的一大挑战。

为此，本文提出了一种基于稀有kmer的串联重复序列比对的位置筛选方法。该方法通过计算待比对序列上稀有kmer构成的序列相似性来衡量待比对序列的相似性，并用于筛选比对结果。此外，该方法将邻近的稀有kmer合并为motif以减少计算量，并采用基于编辑距离的motif模糊查找来避免测序错误与单碱基变异导致的筛选信息丢失。与现有方法相比，在模拟数据集上，本文方法的比对筛选结果具有更高的F1分数及约快两倍的运行速度，具有更强的综合性能。

1 相关研究

1.1 研究背景

序列比对中常用的种子-扩展法可分为过滤与验证两个阶段。在过滤阶段，首先从参考基因组或reads中挑选一些具有代表性的碱基序列片段作为种子，其通常表现为出现频次较低；其次将种子的位置存储在索引中以便快速查找，通常索引的方法有基于hash［15］和基于BWT（Burrows-Wheeler Transform）［16］两种；根据查找出的位置，通过一定策略定位最有可能的位置作为候选位置。在验证阶段，对候选位置进行碱基级别的精确比对，通常采用耗时的动态规划算法［17］，根据精确比对的质量决定最终的比对结果。种子-扩展法流程如图1所示。

基于hash的索引首次在BLAST［15］中被提出，后常应用于BitMapper［18］、FEM［19］等比对算法。hash索引以kmer为键，以kmer的位置为值，只需一次hash函数运算即可定位种子，适合需要频繁定位种子的找全比对问题，但存储所有种子位置的内存开销过大。

BWT在实践中常与后缀数组结合成FM索引，该索引方法对后缀数组中的特定位置进行采样，并利用BWT来求解非采样位置。FM索引仅存储采样位置的做法有效减少了空间消耗，但其定位种子的速度较慢，常应用于找最佳比对算法，如BWA［16］、bowtie2［20］等。

种子-扩展法在串联重复序列比对上的局限体现在两个阶段：

a）在过滤阶段，为避免测序错误影响（三代测序技术错误率在15%左右），种子-扩展法在选取种子时不会选取频率过低的种子，因为其可能代表了测序错误产生的碱基。而受重复影响，频率较高的种子在串联重复区域不具有代表性，导致大量候选比对位置的产生。

b）在验证阶段，为非串联重复序列设定的低比对精度不足以满足串联重复序列要求的高比对精度，因此大量候选位置被认为是正确的比对而输出，产生假阳性比对结果。

1.2 相关工作

实践中常用的是朴素的基于MAPQ的比对过滤方法，该方法考虑利用比对结果本身的信息而非额外信息来过滤比对，例如MAPQ、比对长度、完全匹配的碱基数等体现比对质量的指标。其中MAPQ的大小与比对出错的概率成反比。尽管基于MAPQ的过滤方法准确率很高，但由于串联重复序列上的比对结果的比对质量都相对较低，即使是串联重复序列的正确比对结果也容易被过滤掉，导致该方法在串联重复序列上的召回率较低。

RAfilter［21］是一个基因组重复区域假阳性比对的过滤算法。该算法的主要思想是用reads和参考基因组之间kmer（对于SMRT reads实际采用的是频次为1的uniquekmer）序列的相似性来判断一个比对是假阳性比对的概率。其采用最长递增子序列（longest increase subsequence，LIS）［22］、位操作、二分查找［23］等一系列策略提高计算精度，降低计算时空消耗。RAfilter算法可分为两个主要步骤：

a）统计参考基因组中所有kmer的频次，取出uniquekmer，并通过二进制化kmer，采用位操作遍历参考基因组与所有reads序列以构建uniquekmer的位置索引。

b）针对每一个比对结果，根据其比对区间按上一步构建的位置索引取出区间内的unique kmer向量，按照reads中uniquekmer的顺序计算参考基因组中unique kmer的LIS长度，用其计算KMAPQ以评价比对结果的可信度。KMAPQ的计算公式如下，其中lLIS代表uniquekmer序列的LIS长度，lr与lq分别代表参考基因组与reads中uniquekmer序列的长度：

RAfilter的局限在于，uniquekmer在reads和参考基因组上是稀疏的，若比对区间内缺少uniquekmer，则无法作出有效的过滤。为了提高召回率，RAfilter将这些无法过滤的比对看作是正确的，这使其准确率大幅降低。

2 方法设计

2.1 问题定义

假设ε={A，C，G，T}，分别代表四种脱氧核糖核酸；r={ri|i∈［1，n］，ri∈εli}代表测序得到的n条reads，li代表第i条read的长度；R={Rj|j∈［1，m］，Rj∈εLj}代表参考基因组中的m条参考序列，Lj代表第j条参考序列的长度。若r′i是ri的子序列，R′j是Rj的子序列，则对于输入序列集合r与R，序列比对结果Z{（r′i，R′j）|i∈［1，n］，j∈［1，m］}应输出布尔集合B={bk|k∈［1，Z］}，其中，若zk是正确比对则bk=true，否则bk=false。

2.2 方法描述

串联重复序列的重复单元之间虽然高度相似，但它们也随着时间的推移积累了突变，这些突变以单碱基变异的形式出现在重复单元上，使得重复单元之间存在细微差别。这使得原本因重复而出现频次较高的kmer可能因其序列上某个单碱基变异而降低其出现频次，成为稀有kmer。因此，稀有kmer体现了单碱基变异的存在，能够特异性地区分相似的重复单元。若待比对序列间的稀有kmer能够一一对应，使其稀有kmer的最长公共子序列相对长度大于其他候选序列，那么其是真实比对的可能性将会更大。事实上经观察，在真实数据集上，若某个比对的稀有kmer最长公共子序列相对长度最长，则其恰好是正确的比对位置的概率约是92.17%。

依据上述知识与观察，本文提出的串联重复序列比对的位置筛选方法标记序列上的稀有kmer，在减少测序错误影响的同时增加筛选时的信息密度，将距离上相近的稀有kmer合并成motif来减少计算代价，并采用基于编辑距离的模糊查找来计算motif的最长公共子序列长度，用其作为相似性指标来筛选比对结果。该方法主要分为构建索引、构建motif和计算相似性三个步骤。

2.2.1 构建索引

利用Jellyfish［24］统计参考基因组上所有kmer（默认k=21）的频次，并取出频次小于4的kmer作为稀有kmer构建集合。稀有kmer比uniquekmer更加稠密，使筛选时可用的信息密度增大，提高筛选的准确性。值得注意的是，此处Jellyfish需要加入-C参数，同时统计正反义两条链上的kmer频次，避免遗漏反向互补链上的稀有kmer，使得后续方法能够对比对上反向互补链的序列进行过滤。

若令A=00，C=01，G=10，T=11，则可以将kmer从字符串转换为64位无符号整数。并行遍历reads与参考基因组序列以寻找稀有kmer，记录稀有kmer在序列上起始位置的下标与其64位无符号整数值来构建稀有kmer的有序位置索引。具体来说，在遍历过程中采用上述规则二进制化kmer，并采用滑动窗口、位操作与掩码同时遍历正反向两条链以减少转换的时间消耗。对于每一个滑动窗口中的kmer，计算出正反向kmer对应的64位无符号整数并在稀有kmer集合中查询，若查询成功则记录kmer起始位置下标与对应64位无符号整数。由于稀有kmer集合较大，reads数目较多，读取与查找稀有kmer是时间性能瓶颈。本文采用多线程分块读取稀有kmer与以每条序列为单位的多线程并行策略加速计算。索引结构如图2所示（为简化图例采用k=6）。

2.2.2 构建motif

在得到了稀有kmer的有序位置索引后，给定一个比对结果在序列上的起始位置，能够利用二分查找快速找到比对区间内包含的稀有kmer构成的向量。若一个比对结果是反向比对的（比对到反向互补链上），由于统计kmer时统计了双向链，只需将其中一个稀有kmer向量的顺序倒置即可。

然而，稀有kmer可能过于稠密，使计算单位长度序列相似性的时间代价大幅上涨。因此，考虑将二分查找得到的稀有kmer向量中起始位置间距不大于k的稀有kmer合并成一个motif，在计算序列相似性时以motif而非kmer为单位，并用构成每个motif的所有kmer的64位无符号整数值之和作为它的值。对于reads序列，构建一个向量，其元素为motif的值和构成motif的kmer值的向量。对于参考基因组序列，构建一个字典，其键为motif的值，其值为motif的起始位置和构成motif的kmer的值的向量。数据结构如图3所示。

构建motif减少了后续计算单位长度序列相似性的时间代价。同时，由于单碱基变异与测序错误通常会影响周围若干个kmer的频次，将稀有kmer合并成motif可以避免重复计算，并降低测序错误对序列整体相似性的影响。

2.2.3 计算相似性

在判断两个motif是否相等时，若严格地根据其值是否相等来判断，则motif中某个kmer的缺失或变异就会导致两个motif的值不相等，即使其余的kmer是完全一致的。为解决上述问题，本文通过动态规划算法计算两个motif的稀有kmer向量间的编辑距离，在查找motif时允许编辑距离存在一定误差，从而实现模糊查找。计算编辑距离的动态规划算法状态转移方程如下，其中P={pi|i∈［1，LP］}为reads序列的某个motif的稀有kmer向量，Q={qi|i∈［1，LQ］}是来自参考序列的某个motif的稀有kmer向量，LP与LQ分别代表motif的稀有kmer向量长度，x∈［0，LP］，y∈［0，LQ］。

若满足min（LP，LQ）>3∧dp［LP］［LQ］≤max（2，min（LP，LQ）×0.1），则认为两个motif相近，其中若min（LP，LQ）≤3，则两个motif相等的条件是其值完全相等。

衡量两条序列（在本文中序列的元素是motif）相似性的经典算法是求其最长公共子序列（longest common subsequence，LCS）［25］，但其O（lL）的时间复杂度不可接受，其中l是来自reads的motif序列长度，L是来自参考的motif序列长度。在实际应用中，常求解两条序列的LIS作为代替，其时间代价为O（l+L log L）。求解两条motif序列的LIS步骤如下，示例如图4所示（为简化，示例未考虑稀有kmer向量的最小长度限制）。

a）遍历reads序列的motif向量，对于每个motif，查询参考序列的字典，若其与某个motif相近或相等，则记录其在参考序列上的起始位置。

b）通过动态规划与二分查找计算上一步得到的起始位置的向量的LIS长度len，若len>min（l，L）×0.7，则认为reads与参考序列之间的比对是正确比对，否则为假阳性比对。

求解LIS长度的算法伪代码描述如算法1所示。

算法1 计算公共motif起始位置向量的LIS长度

3 实验结果与分析

为了验证本文方法的效果，将本文方法与朴素的基于比对质量分数（MAPQ）过滤的方法、RAfilter进行比较。朴素的基于MAPQ过滤的方法采用的标准是业内常用标准，正确比对需同时满足以下条件：MAPQ≥45，比对长度大于reads长度的50%，匹配的碱基数需占比对长度90%以上［26］。当运行RAfilter时，除构建位置索引时采用64个线程外，其余参数均设置为默认值。本文方法在构建索引与过滤比对时均采用64个线程。

本文实验的测试平台为：2×AMD EPYC 7H12 64-Core，16×64 GB RAM，Ubuntu 20.04.2 64位操作系统。

3.1 模拟数据实验结果

首先使用reads模拟工具PBSIM2［27］生成一定深度的SMRT reads。由于模拟reads在参考序列上的位置是已知的，所以将真实位置与比对位置进行比较，即可验证上述方法的召回率与准确率。若真实位置与比对位置有交集，且交集的长度占比对长度的99%以上，则认为比对是正确的。

模拟生成的两个串联重复序列数据集reads的覆盖度都是30x，其采样的reads为人类真实SMRT reads，其模拟的参考序列分别为人类T2T基因组的一号染色体着丝粒与二号染色体着丝粒［28］。其中一号染色体着丝粒长5 206 132b，二号染色体着丝粒长2 514 084b，两个不同规模的模拟数据集能更全面地反映方法的效果。实验结果如表1所示。

由表1可知，在两个模拟数据集上，本文方法的召回率最高。本文方法的准确率高于RAfilter，但低于基于MAPQ的方法，因为基于MAPQ的方法筛选硬指标来源于比对本身，且比较严格。综合来看，本文方法的F1分数高于另外两种方法，体现出更为全面的筛选性能。性能的提升主要有两点原因：a）稠密的稀有kmer让原本无法被筛选的reads重新被考虑，在保证召回率的同时提高准确率；b）motif的模糊查找挽回了因测序错误而错位的motif，在保证准确率的同时提高了召回率。

从运行时间代价上来看，本文方法比RAfilter快约两倍，主要归功于将稀有kmer合并为motif节省了大量计算时间。从运行空间代价上来看，本文方法相比于RAfilter占用内存大幅上涨，主要由于稀有kmer相比于uniquekmer存储代价更大，且需要存储更加复杂的数据结构。

3.2 真实数据实验结果

为验证本文方法应用效果与观察结果是否一致，并分析比对位置数量对本文方法的影响，本文创建了一个真实数据集实例。将人类基因组真实的SMRT reads比对到人类T2T基因组的一号染色体着丝粒后，选择其中比对到唯一位置的reads，并记录这部分真实数据集在参考基因组上的位置。接着缩短这部分reads的长度为原来的1/4，缩短长度后的reads易比对到参考基因组的多个位置，此时可以使用之前记录的位置作为实际比对位置对筛选结果进行验证。

构建的真实数据集中，比对到唯一位置的reads有770 714条，缩短长度后有2 751 396条reads比对到了多达9 356 186个位置，其中，比对上的位置数量为1、2、3、4及4以上的reads分别有397 316（31%）、551 805（10%）、583 791（4%）、1 218 484（55%）条。在真实数据集上针对比对位置数量不同的reads，本文方法实验结果如表2所示。

由表2可知，本文方法在真实数据集上整体的召回率达91.65%，接近于2.2节提到的基于观察的概率92.17%，且不同比对位置数量上的召回率都分布于整体召回率附近，证明本文方法的实际过滤效果符合观察规律，即稀有kmer最长公共子序列相对长度最长的比对位置就是正确的比对位置。

此外，在比对位置数量不同的reads上，本文方法的过滤效果比较稳定，不存在因比对位置数量不同而导致的过滤效果偏差，且明显优于随机选择的筛选效果（例如，存在2个比对位置时随机选择到正确比对的概率是50%），证明本文方法在实际应用中是切实有效的。

3.3 并行加速实验结果

为验证本文方法并行加速效果，在上述人类一号染色体着丝粒模拟数据集上分别采用1、8、16、32、64个线程运行RAfilter与本文方法。实验结果如表3所示。

从表3可以看出，由于RAfilter只有构建索引的步骤并行化，其运行时间虽然随线程数增加而下降，但在筛选的阶段存在瓶颈。而本文方法虽然在单线程时比RAfilter慢，但随着线程数增加，其运行时间的下降幅度显著大于RAfilter，体现出其良好的并行加速效果。

4 结束语

本文针对串联重复序列存在假阳性比对的问题提出了一种基于稀有kmer的比对位置筛选方法，通过利用稀有kmer的稠密性增大信息密度，将kmer合并成motif来减少计算时间代价，采用基于编辑距离的模糊查找增强motif对错误的容忍性，计算motif序列的LIS来衡量序列相似性用于筛选比对位置，从而大幅降低比对假阳性的可能。实验结果表明，本文方法采用的稀有kmer和模糊查找能够提高比对筛选的召回率与准确率，同时以motif为计算单位能显著减少计算时间，在加速筛选的同时具有良好的并行性能。由于存储了数量更多的稀有kmer与相关数据结构，本文方法相比其他算法具有更大的内存开销，如何改进相关数据结构从而减少内存开销，是下一步的工作重点。

参考文献：

［1］Wenger A M， Peluso P， Rowell W J， et al. Accurate circular consensus long-read sequencing improves variant detection and assembly of a human genome［J］. Nature Biotechnology， 2019， 37（10）： 1155-1162.

［2］Jain M， Olsen H E， Paten B， et al. The Oxford Nanopore MinION： delivery of nanopore sequencing to the genomics community［J］. Genome Biology， 2016，17（1）： 239.

［3］Goodwin S， Mcpherson J D， Mccombie W R. Coming of age： ten years of next-generation sequencing technologies［J］. Nature Reviews Genetics， 2016， 17（6）： 333-351.

［4］Benson G. Tandem repeats finder， a program to analyze DNA［J］. Nucleic Acids Research， 1999， 27（2）： 573-580.

［5］Rice E S，Green R E. New approaches for genome assembly and scaffolding［J］.Annual Review of Animal Biosciences， 2019， 7： 17-40.

［6］Kazazian H H. Mobile elements： drivers of genome evolution［J］. Science， 2004， 303（5664）： 1626-1632.

［7］Bzikadze A V， Pevzner P A. Automated assembly of centromeres from ultra-long error-prone reads［J］. Nature Biotechnology， 2020， 38（11）： 1309-1316.

［8］Cheng Haoyu， Concepcion G T， Feng Xiaowen， et al. Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm［J］. Nature Methods， 2021， 18（2）： 170-175.

［9］Poplin R， Chang Pichuan， Alexander D， et al. A universal SNP and small-indel variant caller using deep neural networks［J］. Nature Biotechnology， 2018， 36（10）： 983-987.

［10］Wlodzimierz P， Rabanal F A， Burns R， et al. Cycles of satellite and transposon evolution in Arabidopsis centromeres［J］. Nature， 2023，618：557-565.

［11］Sedlazeck F J， Lee H， Darby C A， et al. Piercing the dark matter： bioinformatics of long-range sequencing and mapping［J］. Nature Reviews Genetics， 2018， 19（6）： 329-346.

［12］Ahmed N， Bertels K， Al-Ars Z. A comparison of seed-and-extend techniques in modern DNA read alignment algorithms［C］//Proc of IEEE International Conference on Bioinformatics and Biomedicine. Piscataway， NJ： IEEE Press， 2016： 1421-1428.

［13］Bzikadze A V， Pevzner P A. UniAligner： a parameter-free framework for fast sequence alignment［J］. Nature Methods， 2023， 20（9）： 1346-1354.

［14］Li Heng. Minimap2： pairwise alignment for nucleotide sequences［J］. Bioinformatics， 2018， 34（18）： 3094-3100.

［15］Altschul S G W， Miller W， Myers E W， et al. Basic local alignment search tool［J］. Journal of Molecular Biology， 1990， 215（3）： 403-410.

［16］Li Heng， Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform［J］. Bioinformatics， 2009， 25（14）： 1754-1760.

［17］Smith T W M. Identification of common molecular subsequences［J］. Journal of Molecular Biology， 1981， 147（1）： 195-197.

［18］Cheng Haoyu， Jiang Huaipan， Yang Jiaoyun， et al. BitMapper： an efficient all-mapper based on bit-vector computing［J］. BMC Bioinformatics， 2015， 16（1）： 192.

［19］Zhang Haowen， Chan Yuandong， Fan Kaichao， et al. Fast and efficient short read mapping based on a succinct hash index［J］. BMC Bioinformatics， 2018， 19（1）： 92.

［20］Langmead B， Salzberg S L. Fast gapped-read alignment with Bowtie 2［J］. Nature Methods， 2012， 9（4）： 357-359.

［21］Yang Jinbao， Zhao Xianjia， Jiang Heling， et al. RAfilter： an algorithm for detecting and filtering false-positive alignments in repetitive genomic regions［J］. Horticulture Research， 2023， 10（1）： 288.

［22］Fredman M L. On computing the length of longest increasing subsequences［J］. Discrete Mathematics， 1975， 11（1）： 29-35.

［23］Bentley J L. Multidimensional binary search trees used for associative searching［J］. Communications of the ACM， 1975， 18（9）： 509-517.

［24］Marcais G， Kingsford C. A fast， lock-free approach for efficient parallel counting of occurrences of kmers［J］. Bioinformatics， 2011， 27（6）： 764-770.

［25］Levenshtein V. Binary codes capable of correcting deletions， insertions， and reversals［J］. Soviet Physics Doklady， 1966， 10（8）： 707-710.

［26］Jarvis E D， Formenti G， Rhie A， et al. Semi-automated assembly of high-quality diploid human reference genomes［J］. Nature， 2022，611：519-531.

［27］Ono Y， Asai K， Hamada M. PBSIM2： a simulator for long-read sequencers with a novel generative model of quality scores［J］. Bioinformatics， 2021，37（5）： 589-595.

［28］Nurk S， Koren S， Rhie A， et al. The complete sequence of a human genome［J］. Science， 2022， 376（6588）： 44-53.