摘 要：串聯(lián)重復(fù)序列是基因組構(gòu)建的困難片段，由于其重復(fù)單元之間的相似性與其拷貝數(shù)的不確定性，在序列比對時容易定位到多個候選位置，如何快速而準(zhǔn)確地篩選出正確的比對位置是一項挑戰(zhàn)?，F(xiàn)有方法使用種子（從測序片段中選取的短序列）來定位并擴展候選比對位置，但挑選種子時未考慮串聯(lián)重復(fù)序列特性。因此，提出了一種串聯(lián)重復(fù)序列比對的位置篩選方法，其通過計算稀有kmer（長度為k的子序列）序列的相似性來篩選比對結(jié)果。此外，采用合并稀有kmer的策略加速計算，并利用基于編輯距離的模糊查找以提高過濾信息密度。實驗結(jié)果表明，在模擬數(shù)據(jù)集上提高比對結(jié)果的召回率與準(zhǔn)確率的同時，該方法比現(xiàn)有方法快約2倍，且具有良好的并行加速性能。

關(guān)鍵詞：串聯(lián)重復(fù)； 單分子實時測序； 序列比對； 種子-擴展法

中圖分類號：TP391 文獻標(biāo)志碼：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 引言

隨著基因測序技術(shù)的高速發(fā)展，新興的三代測序技術(shù)，如單分子實時測序技術(shù)（single molecule real-time，SMRT）［1］和牛津納米孔測序技術(shù)（Oxford nanopore technology，ONT）［2］，連同二代測序技術(shù)（next generation sequencing，NGS）［3］一起助力了越來越多端粒到端粒T2T（telomere to telomere）基因組的誕生。其中，串聯(lián)重復(fù)序列是一種由高度相似的重復(fù)單元緊密排列而成的特殊結(jié)構(gòu)，構(gòu)成了T2T基因組中的許多重要組件，例如著絲粒、端粒、rDNA等［4］，是T2T基因組中最復(fù)雜的結(jié)構(gòu)之一。研究表明，真核生物的基因組很大部分由重復(fù)區(qū)域構(gòu)成［5］，特別是哺乳動物中重復(fù)區(qū)域約占整個基因組的25%～50%［6］。然而，重復(fù)單元之間的相似性與其拷貝數(shù)的不確定性為許多T2T基因組分析問題在計算上帶來了挑戰(zhàn)［7］，其中一個關(guān)鍵問題就是序列比對。許多生物信息學(xué)問題的解決，例如denovo組裝［8］、結(jié)構(gòu)變異檢測［9］與進化分析［10］等都依賴于序列比對結(jié)果的準(zhǔn)確性［11］。

在生物信息學(xué)中，序列比對是重要的基礎(chǔ)性問題之一。基因序列比對的問題定義為利用相關(guān)算法，將基因測序得到的reads比對到已知基因序列上，得到其相對位置關(guān)系，從而判斷序列相似性與相關(guān)生物學(xué)特征。主流的序列比對算法普遍采用種子-擴展法（seed-and-extend）［12］，該方法通過種子索引來定位比對位置，將比對范圍從全基因組縮小至若干候選位置，再對候選位置進行精確的比對，從而減少整體時間消耗。在面對串聯(lián)重復(fù)序列時，大多數(shù)通用序列比對算法選取的種子不足以區(qū)分重復(fù)單元之間的細微差別，導(dǎo)致定位到多個候選比對位置，進而產(chǎn)生大量假陽性比對結(jié)果，形成多重比對現(xiàn)象［13］。例如，Minimap2［14］在將7 532條來自串聯(lián)重復(fù)區(qū)域的reads比對到參考基因組時產(chǎn)生了多達29 106條比對結(jié)果，對篩選出正確的比對位置造成了困難。因此，設(shè)計高效的序列比對位置篩選方法，以過濾串聯(lián)重復(fù)序列的假陽性比對結(jié)果，是基因序列比對面臨的一大挑戰(zhàn)。

為此，本文提出了一種基于稀有kmer的串聯(lián)重復(fù)序列比對的位置篩選方法。該方法通過計算待比對序列上稀有kmer構(gòu)成的序列相似性來衡量待比對序列的相似性，并用于篩選比對結(jié)果。此外，該方法將鄰近的稀有kmer合并為motif以減少計算量，并采用基于編輯距離的motif模糊查找來避免測序錯誤與單堿基變異導(dǎo)致的篩選信息丟失。與現(xiàn)有方法相比，在模擬數(shù)據(jù)集上，本文方法的比對篩選結(jié)果具有更高的F1分數(shù)及約快兩倍的運行速度，具有更強的綜合性能。

1 相關(guān)研究

1.1 研究背景

序列比對中常用的種子-擴展法可分為過濾與驗證兩個階段。在過濾階段，首先從參考基因組或reads中挑選一些具有代表性的堿基序列片段作為種子，其通常表現(xiàn)為出現(xiàn)頻次較低；其次將種子的位置存儲在索引中以便快速查找，通常索引的方法有基于hash［15］和基于BWT（Burrows-Wheeler Transform）［16］兩種；根據(jù)查找出的位置，通過一定策略定位最有可能的位置作為候選位置。在驗證階段，對候選位置進行堿基級別的精確比對，通常采用耗時的動態(tài)規(guī)劃算法［17］，根據(jù)精確比對的質(zhì)量決定最終的比對結(jié)果。種子-擴展法流程如圖1所示。

基于hash的索引首次在BLAST［15］中被提出，后常應(yīng)用于BitMapper［18］、FEM［19］等比對算法。hash索引以kmer為鍵，以kmer的位置為值，只需一次hash函數(shù)運算即可定位種子，適合需要頻繁定位種子的找全比對問題，但存儲所有種子位置的內(nèi)存開銷過大。

BWT在實踐中常與后綴數(shù)組結(jié)合成FM索引，該索引方法對后綴數(shù)組中的特定位置進行采樣，并利用BWT來求解非采樣位置。FM索引僅存儲采樣位置的做法有效減少了空間消耗，但其定位種子的速度較慢，常應(yīng)用于找最佳比對算法，如BWA［16］、bowtie2［20］等。

種子-擴展法在串聯(lián)重復(fù)序列比對上的局限體現(xiàn)在兩個階段：

a）在過濾階段，為避免測序錯誤影響（三代測序技術(shù)錯誤率在15%左右），種子-擴展法在選取種子時不會選取頻率過低的種子，因為其可能代表了測序錯誤產(chǎn)生的堿基。而受重復(fù)影響，頻率較高的種子在串聯(lián)重復(fù)區(qū)域不具有代表性，導(dǎo)致大量候選比對位置的產(chǎn)生。

b）在驗證階段，為非串聯(lián)重復(fù)序列設(shè)定的低比對精度不足以滿足串聯(lián)重復(fù)序列要求的高比對精度，因此大量候選位置被認為是正確的比對而輸出，產(chǎn)生假陽性比對結(jié)果。

1.2 相關(guān)工作

實踐中常用的是樸素的基于MAPQ的比對過濾方法，該方法考慮利用比對結(jié)果本身的信息而非額外信息來過濾比對，例如MAPQ、比對長度、完全匹配的堿基數(shù)等體現(xiàn)比對質(zhì)量的指標(biāo)。其中MAPQ的大小與比對出錯的概率成反比。盡管基于MAPQ的過濾方法準(zhǔn)確率很高，但由于串聯(lián)重復(fù)序列上的比對結(jié)果的比對質(zhì)量都相對較低，即使是串聯(lián)重復(fù)序列的正確比對結(jié)果也容易被過濾掉，導(dǎo)致該方法在串聯(lián)重復(fù)序列上的召回率較低。

RAfilter［21］是一個基因組重復(fù)區(qū)域假陽性比對的過濾算法。該算法的主要思想是用reads和參考基因組之間kmer（對于SMRT reads實際采用的是頻次為1的uniquekmer）序列的相似性來判斷一個比對是假陽性比對的概率。其采用最長遞增子序列（longest increase subsequence，LIS）［22］、位操作、二分查找［23］等一系列策略提高計算精度，降低計算時空消耗。RAfilter算法可分為兩個主要步驟：

a）統(tǒng)計參考基因組中所有kmer的頻次，取出uniquekmer，并通過二進制化kmer，采用位操作遍歷參考基因組與所有reads序列以構(gòu)建uniquekmer的位置索引。

b）針對每一個比對結(jié)果，根據(jù)其比對區(qū)間按上一步構(gòu)建的位置索引取出區(qū)間內(nèi)的unique kmer向量，按照reads中uniquekmer的順序計算參考基因組中unique kmer的LIS長度，用其計算KMAPQ以評價比對結(jié)果的可信度。KMAPQ的計算公式如下，其中l(wèi)LIS代表uniquekmer序列的LIS長度，lr與lq分別代表參考基因組與reads中uniquekmer序列的長度：

RAfilter的局限在于，uniquekmer在reads和參考基因組上是稀疏的，若比對區(qū)間內(nèi)缺少uniquekmer，則無法作出有效的過濾。為了提高召回率，RAfilter將這些無法過濾的比對看作是正確的，這使其準(zhǔn)確率大幅降低。

2 方法設(shè)計

2.1 問題定義

假設(shè)ε={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，序列比對結(jié)果Z{（r′i，R′j）|i∈［1，n］，j∈［1，m］}應(yīng)輸出布爾集合B={bk|k∈［1，Z］}，其中，若zk是正確比對則bk=true，否則bk=false。

2.2 方法描述

串聯(lián)重復(fù)序列的重復(fù)單元之間雖然高度相似，但它們也隨著時間的推移積累了突變，這些突變以單堿基變異的形式出現(xiàn)在重復(fù)單元上，使得重復(fù)單元之間存在細微差別。這使得原本因重復(fù)而出現(xiàn)頻次較高的kmer可能因其序列上某個單堿基變異而降低其出現(xiàn)頻次，成為稀有kmer。因此，稀有kmer體現(xiàn)了單堿基變異的存在，能夠特異性地區(qū)分相似的重復(fù)單元。若待比對序列間的稀有kmer能夠一一對應(yīng)，使其稀有kmer的最長公共子序列相對長度大于其他候選序列，那么其是真實比對的可能性將會更大。事實上經(jīng)觀察，在真實數(shù)據(jù)集上，若某個比對的稀有kmer最長公共子序列相對長度最長，則其恰好是正確的比對位置的概率約是92.17%。

依據(jù)上述知識與觀察，本文提出的串聯(lián)重復(fù)序列比對的位置篩選方法標(biāo)記序列上的稀有kmer，在減少測序錯誤影響的同時增加篩選時的信息密度，將距離上相近的稀有kmer合并成motif來減少計算代價，并采用基于編輯距離的模糊查找來計算motif的最長公共子序列長度，用其作為相似性指標(biāo)來篩選比對結(jié)果。該方法主要分為構(gòu)建索引、構(gòu)建motif和計算相似性三個步驟。

2.2.1 構(gòu)建索引

利用Jellyfish［24］統(tǒng)計參考基因組上所有kmer（默認k=21）的頻次，并取出頻次小于4的kmer作為稀有kmer構(gòu)建集合。稀有kmer比uniquekmer更加稠密，使篩選時可用的信息密度增大，提高篩選的準(zhǔn)確性。值得注意的是，此處Jellyfish需要加入-C參數(shù)，同時統(tǒng)計正反義兩條鏈上的kmer頻次，避免遺漏反向互補鏈上的稀有kmer，使得后續(xù)方法能夠?qū)Ρ葘ι戏聪蚧パa鏈的序列進行過濾。

若令A(yù)=00，C=01，G=10，T=11，則可以將kmer從字符串轉(zhuǎn)換為64位無符號整數(shù)。并行遍歷reads與參考基因組序列以尋找稀有kmer，記錄稀有kmer在序列上起始位置的下標(biāo)與其64位無符號整數(shù)值來構(gòu)建稀有kmer的有序位置索引。具體來說，在遍歷過程中采用上述規(guī)則二進制化kmer，并采用滑動窗口、位操作與掩碼同時遍歷正反向兩條鏈以減少轉(zhuǎn)換的時間消耗。對于每一個滑動窗口中的kmer，計算出正反向kmer對應(yīng)的64位無符號整數(shù)并在稀有kmer集合中查詢，若查詢成功則記錄kmer起始位置下標(biāo)與對應(yīng)64位無符號整數(shù)。由于稀有kmer集合較大，reads數(shù)目較多，讀取與查找稀有kmer是時間性能瓶頸。本文采用多線程分塊讀取稀有kmer與以每條序列為單位的多線程并行策略加速計算。索引結(jié)構(gòu)如圖2所示（為簡化圖例采用k=6）。

2.2.2 構(gòu)建motif

在得到了稀有kmer的有序位置索引后，給定一個比對結(jié)果在序列上的起始位置，能夠利用二分查找快速找到比對區(qū)間內(nèi)包含的稀有kmer構(gòu)成的向量。若一個比對結(jié)果是反向比對的（比對到反向互補鏈上），由于統(tǒng)計kmer時統(tǒng)計了雙向鏈，只需將其中一個稀有kmer向量的順序倒置即可。

然而，稀有kmer可能過于稠密，使計算單位長度序列相似性的時間代價大幅上漲。因此，考慮將二分查找得到的稀有kmer向量中起始位置間距不大于k的稀有kmer合并成一個motif，在計算序列相似性時以motif而非kmer為單位，并用構(gòu)成每個motif的所有kmer的64位無符號整數(shù)值之和作為它的值。對于reads序列，構(gòu)建一個向量，其元素為motif的值和構(gòu)成motif的kmer值的向量。對于參考基因組序列，構(gòu)建一個字典，其鍵為motif的值，其值為motif的起始位置和構(gòu)成motif的kmer的值的向量。數(shù)據(jù)結(jié)構(gòu)如圖3所示。

構(gòu)建motif減少了后續(xù)計算單位長度序列相似性的時間代價。同時，由于單堿基變異與測序錯誤通常會影響周圍若干個kmer的頻次，將稀有kmer合并成motif可以避免重復(fù)計算，并降低測序錯誤對序列整體相似性的影響。

2.2.3 計算相似性

在判斷兩個motif是否相等時，若嚴格地根據(jù)其值是否相等來判斷，則motif中某個kmer的缺失或變異就會導(dǎo)致兩個motif的值不相等，即使其余的kmer是完全一致的。為解決上述問題，本文通過動態(tài)規(guī)劃算法計算兩個motif的稀有kmer向量間的編輯距離，在查找motif時允許編輯距離存在一定誤差，從而實現(xiàn)模糊查找。計算編輯距離的動態(tài)規(guī)劃算法狀態(tài)轉(zhuǎn)移方程如下，其中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）相似性的經(jīng)典算法是求其最長公共子序列（longest common subsequence，LCS）［25］，但其O（lL）的時間復(fù)雜度不可接受，其中l(wèi)是來自reads的motif序列長度，L是來自參考的motif序列長度。在實際應(yīng)用中，常求解兩條序列的LIS作為代替，其時間代價為O（l+L log L）。求解兩條motif序列的LIS步驟如下，示例如圖4所示（為簡化，示例未考慮稀有kmer向量的最小長度限制）。

a）遍歷reads序列的motif向量，對于每個motif，查詢參考序列的字典，若其與某個motif相近或相等，則記錄其在參考序列上的起始位置。

b）通過動態(tài)規(guī)劃與二分查找計算上一步得到的起始位置的向量的LIS長度len，若len>min（l，L）×0.7，則認為reads與參考序列之間的比對是正確比對，否則為假陽性比對。

求解LIS長度的算法偽代碼描述如算法1所示。

算法1 計算公共motif起始位置向量的LIS長度

3 實驗結(jié)果與分析

為了驗證本文方法的效果，將本文方法與樸素的基于比對質(zhì)量分數(shù)（MAPQ）過濾的方法、RAfilter進行比較。樸素的基于MAPQ過濾的方法采用的標(biāo)準(zhǔn)是業(yè)內(nèi)常用標(biāo)準(zhǔn)，正確比對需同時滿足以下條件：MAPQ≥45，比對長度大于reads長度的50%，匹配的堿基數(shù)需占比對長度90%以上［26］。當(dāng)運行RAfilter時，除構(gòu)建位置索引時采用64個線程外，其余參數(shù)均設(shè)置為默認值。本文方法在構(gòu)建索引與過濾比對時均采用64個線程。

本文實驗的測試平臺為：2×AMD EPYC 7H12 64-Core，16×64 GB RAM，Ubuntu 20.04.2 64位操作系統(tǒng)。

3.1 模擬數(shù)據(jù)實驗結(jié)果

首先使用reads模擬工具PBSIM2［27］生成一定深度的SMRT reads。由于模擬reads在參考序列上的位置是已知的，所以將真實位置與比對位置進行比較，即可驗證上述方法的召回率與準(zhǔn)確率。若真實位置與比對位置有交集，且交集的長度占比對長度的99%以上，則認為比對是正確的。

模擬生成的兩個串聯(lián)重復(fù)序列數(shù)據(jù)集reads的覆蓋度都是30x，其采樣的reads為人類真實SMRT reads，其模擬的參考序列分別為人類T2T基因組的一號染色體著絲粒與二號染色體著絲粒［28］。其中一號染色體著絲粒長5 206 132b，二號染色體著絲粒長2 514 084b，兩個不同規(guī)模的模擬數(shù)據(jù)集能更全面地反映方法的效果。實驗結(jié)果如表1所示。

由表1可知，在兩個模擬數(shù)據(jù)集上，本文方法的召回率最高。本文方法的準(zhǔn)確率高于RAfilter，但低于基于MAPQ的方法，因為基于MAPQ的方法篩選硬指標(biāo)來源于比對本身，且比較嚴格。綜合來看，本文方法的F1分數(shù)高于另外兩種方法，體現(xiàn)出更為全面的篩選性能。性能的提升主要有兩點原因：a）稠密的稀有kmer讓原本無法被篩選的reads重新被考慮，在保證召回率的同時提高準(zhǔn)確率；b）motif的模糊查找挽回了因測序錯誤而錯位的motif，在保證準(zhǔn)確率的同時提高了召回率。

從運行時間代價上來看，本文方法比RAfilter快約兩倍，主要歸功于將稀有kmer合并為motif節(jié)省了大量計算時間。從運行空間代價上來看，本文方法相比于RAfilter占用內(nèi)存大幅上漲，主要由于稀有kmer相比于uniquekmer存儲代價更大，且需要存儲更加復(fù)雜的數(shù)據(jù)結(jié)構(gòu)。

3.2 真實數(shù)據(jù)實驗結(jié)果

為驗證本文方法應(yīng)用效果與觀察結(jié)果是否一致，并分析比對位置數(shù)量對本文方法的影響，本文創(chuàng)建了一個真實數(shù)據(jù)集實例。將人類基因組真實的SMRT reads比對到人類T2T基因組的一號染色體著絲粒后，選擇其中比對到唯一位置的reads，并記錄這部分真實數(shù)據(jù)集在參考基因組上的位置。接著縮短這部分reads的長度為原來的1/4，縮短長度后的reads易比對到參考基因組的多個位置，此時可以使用之前記錄的位置作為實際比對位置對篩選結(jié)果進行驗證。

構(gòu)建的真實數(shù)據(jù)集中，比對到唯一位置的reads有770 714條，縮短長度后有2 751 396條reads比對到了多達9 356 186個位置，其中，比對上的位置數(shù)量為1、2、3、4及4以上的reads分別有397 316（31%）、551 805（10%）、583 791（4%）、1 218 484（55%）條。在真實數(shù)據(jù)集上針對比對位置數(shù)量不同的reads，本文方法實驗結(jié)果如表2所示。

由表2可知，本文方法在真實數(shù)據(jù)集上整體的召回率達91.65%，接近于2.2節(jié)提到的基于觀察的概率92.17%，且不同比對位置數(shù)量上的召回率都分布于整體召回率附近，證明本文方法的實際過濾效果符合觀察規(guī)律，即稀有kmer最長公共子序列相對長度最長的比對位置就是正確的比對位置。

此外，在比對位置數(shù)量不同的reads上，本文方法的過濾效果比較穩(wěn)定，不存在因比對位置數(shù)量不同而導(dǎo)致的過濾效果偏差，且明顯優(yōu)于隨機選擇的篩選效果（例如，存在2個比對位置時隨機選擇到正確比對的概率是50%），證明本文方法在實際應(yīng)用中是切實有效的。

3.3 并行加速實驗結(jié)果

為驗證本文方法并行加速效果，在上述人類一號染色體著絲粒模擬數(shù)據(jù)集上分別采用1、8、16、32、64個線程運行RAfilter與本文方法。實驗結(jié)果如表3所示。

從表3可以看出，由于RAfilter只有構(gòu)建索引的步驟并行化，其運行時間雖然隨線程數(shù)增加而下降，但在篩選的階段存在瓶頸。而本文方法雖然在單線程時比RAfilter慢，但隨著線程數(shù)增加，其運行時間的下降幅度顯著大于RAfilter，體現(xiàn)出其良好的并行加速效果。

4 結(jié)束語

本文針對串聯(lián)重復(fù)序列存在假陽性比對的問題提出了一種基于稀有kmer的比對位置篩選方法，通過利用稀有kmer的稠密性增大信息密度，將kmer合并成motif來減少計算時間代價，采用基于編輯距離的模糊查找增強motif對錯誤的容忍性，計算motif序列的LIS來衡量序列相似性用于篩選比對位置，從而大幅降低比對假陽性的可能。實驗結(jié)果表明，本文方法采用的稀有kmer和模糊查找能夠提高比對篩選的召回率與準(zhǔn)確率，同時以motif為計算單位能顯著減少計算時間，在加速篩選的同時具有良好的并行性能。由于存儲了數(shù)量更多的稀有kmer與相關(guān)數(shù)據(jù)結(jié)構(gòu)，本文方法相比其他算法具有更大的內(nèi)存開銷，如何改進相關(guān)數(shù)據(jù)結(jié)構(gòu)從而減少內(nèi)存開銷，是下一步的工作重點。

參考文獻：

［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.