基于長讀數(shù)和多序列比對的間隙填充方法

毋 東，魏亞偉，羅軍偉，敖 山

（河南理工大學計算機科學與技術學院，河南焦作 454000）

0 概述

隨著人類基因組計劃的完成，基因組數(shù)據(jù)成為后基因組時代的重要標志，這些海量基因組數(shù)據(jù)的分析和挖掘深入到生命科學研究的各個方面，是眾多生命科學研究的基礎，更是計算機科學特別是生物信息學研究所面臨的重大挑戰(zhàn)。《國家“十三五”科技創(chuàng)新規(guī)劃》和《國家自然科學基金委“十三五”發(fā)展規(guī)劃》都明確將基因組數(shù)據(jù)分析以及生物大數(shù)據(jù)分析作為以“重大計劃”形式倡導的研究內容。

目前，隨著全基因組測序技術的發(fā)展，人們利用不同的序列組裝方法獲得越來越連續(xù)和完整的基因組序列。但是由于測序錯誤［1］、基因組重復區(qū)［2］、多態(tài)性等問題的存在，現(xiàn)有的序列組裝方法往往輸出一個包含間隙（gap）的序列集合，即scaffold［3］集合，其中每一個scaffold 是一條堿基序列，但是其包含一些gap 區(qū)域，通常由連續(xù)的“N”表示。而這些gap 區(qū)域可能對應基因編碼或調控區(qū)域，這給下游分析中的基因表達和調控、結構變異分析、物種進化等研究增加困難。

為減少scaffold 中gap 的數(shù)量，增加序列的連續(xù)性和完整性，研究人員設計了gap 填充方法。gap 填充方法利用由測序技術獲得的讀數(shù)（read/序列片段）集合，確定scaffold 集合中gap 區(qū)域的序列。當前流行的測序技術主要包括第二代和第三代測序技術。第二代測序技術具有高速度、高通量、低成本等優(yōu)點，但是第二代測序技術產(chǎn)生的讀數(shù)長度較短，在一百到五百堿基左右，越來越多的研究表明，短讀數(shù)對于復雜重復區(qū)的gap 區(qū)域的填充能力非常有限。而以太平洋生物科學公司的單分子實時測序技術［4］和牛津納米技術公司的納米孔單分子技術［5］為代表的第三代測序技術，能夠產(chǎn)生長度達到數(shù)萬堿基的長讀數(shù)。這些長讀數(shù)可以跨過基因組序列中大部分的重復區(qū)，有效彌補第二代測序技術的不足。但是第三代測序技術的測序錯誤率達到15%左右，這對gap 填充的質量和效率造成較大的影響。

本文提出一種基于長讀數(shù)和多序列比對的gap填充方法gapLM。將scaffold 集合轉化為不含gap 的序列（contig）集合，并利用比對工具把長讀數(shù)集合比對到contig 集合上，然后對獲得的比對結果進行修正。針對一個具體的gap，識別和提取覆蓋該gap 區(qū)域的長讀數(shù)序列，并將這些提取出來的序列分為3 個不同集合。最后利用多序列比對方法和新的策略對3 個不同的序列集合進行一致化，并對該gap 區(qū)域進行填充。

1 相關工作

近年來，隨著測序技術的發(fā)展，研究人員提出許多gap 填充方法。根據(jù)gap 填充方法使用的讀數(shù)類型，本文將gap 填充方法分為兩類：一類是基于雙端短讀數(shù)的gap 填充方法；另一類是基于長讀數(shù)的gap填充方法。gap 填充示例如圖1 所示。

圖1 gap 填充示例Fig.1 Example of gap filling

1）基于雙端短讀數(shù)的gap 填充方法

基于雙端短讀數(shù)的gap 填充方法往往利用短讀數(shù)正確率高和組裝工具成熟等特點對gap 進行填充，其主要思路是基于讀數(shù)集合和scaffold 集合的比對信息確定和gap 關聯(lián)的讀數(shù)集合，然后再采用不同的策略對gap 進行填充，如圖1（a）所示。一些gap 填充方法利用讀數(shù)或者k-mer（長度為k的子序列）之間的重疊關系，從gap 區(qū)域緊鄰兩端進行擴展找到覆蓋gap 區(qū)域的序列。文獻［6］提出一種IMAGE 的gap填充方法，該方法使用比對工具將雙端讀數(shù)比對到scaffold 集合上，收集一端讀數(shù)能夠比對到gap 緊鄰區(qū)域的雙端讀數(shù)集合，然后通過Velvet［7］對讀數(shù)集合進行組裝，構建contig（更長的序列）并填充gap 區(qū)域。文獻［8］提出一種gapCloser 方法，首先確定和gap 相鄰區(qū)域有重疊關系的讀數(shù)集合，并在局部組裝時加入容錯機制［9］，從而提高填充gap 的準確性。文獻［10］提出一種gapFiller 方法，該方法將對應到gap區(qū)域的讀數(shù)轉換成k-mer 集合，然后從gap 左右兩端選擇具有重疊關系的k-mer 進行合并，反復擴展逐漸縮小gap 區(qū)域。

此外，一些啟發(fā)式方法也被用于gap 填充的過程。此類方法一般先構建讀數(shù)重疊圖或者De Bruijn圖，然后從圖中選擇一條路徑代表gap 區(qū)域的序列。文獻［11］提出一種FinIS 方法，該方法針對每個gap區(qū)域，建立一個讀數(shù)重疊圖，然后采用混合整數(shù)規(guī)劃方法和路徑評分函數(shù)確定最佳路徑。文獻［12］提出一種Sealer 方法，該方法針對一個gap，確定一個讀數(shù)集合和相對應的k-mer 集合，使用布隆過濾器數(shù)據(jù)結構存儲這些k-mer 集合，并構建De Bruijn 圖，在De Bruijn 圖中利用廣度優(yōu)先搜索尋找連接gap 區(qū)域兩側序列的路徑。文獻［13］提出一種gap2Seq 方法，該方法是在構建的De Bruijn 圖中確定最接近gap 長度的路徑，由于該過程是NP-hard 問題，因此gap2Seq方法基于一種簡單的動態(tài)規(guī)劃算法尋找近似解。文獻［14］提出一種gapReduce 方法，該方法通過迭代改變k-mer 頻次和長度閾值構建De Bruijn 圖，并基于雙端讀數(shù)的insert size 分布設計評分函數(shù)選擇路徑。

2）基于長讀數(shù)的gap 填充方法

基于長讀數(shù)的gap 填充方法的基本思想是利用比對方法將長讀數(shù)與scaffold 集合進行比對，然后選擇gap 對應的長讀數(shù)去填充，如圖1（b）所示。文獻［15］提出一種PBJelly 方法，該方法采用BLASR［16］將長讀 數(shù)比對 到scaffold 集合上，識別覆蓋gap 區(qū)域的長讀數(shù)集合，并基于OLC［17-18］組裝算法找到填充gap 區(qū)域的一致序列。文獻［19］提出一種GMcloser方法，該方法首先將scaffold 中的gap 區(qū)域刪除，使其變成不包含gap 的contig 集合，然后將這些contig 與長讀數(shù)集合進行比對。該方法基于比對概率估算模型，選擇最大可能比對到gap 區(qū)域的長讀數(shù)序列進行填充。文獻［20］提出一種LR_gapcloser方法，該方法將每個長讀數(shù)分成長度相同且有序的短序列片段，然后使用BWA-MEM［21］算法將所有的短序列片段比對到scaffold 上。該方法基于堿基覆蓋度評分函數(shù)以及短序列片段的比對方向和順序，過濾掉一些錯誤比對，進而提高填充的準確性。

2 本文方法

本文提出一種基于gapLM 的gap 填充方法，該方法以長讀數(shù)集合和scaffold 集合作為輸入，通過分析長讀數(shù)集合和scaffold 集合之間的比對信息，對gap 區(qū)域進行填充。該方法主要包括以下4 個步驟：1）刪除scaffold 中的gap 區(qū)域，獲得一個新的contig集合，然后利用比對工具BWA 將長讀數(shù)比對到contig 集合上；2）通過分析contig 和長讀數(shù)的之間比對位置的差異，對比對結果進行修正；3）根據(jù)上一個步驟獲得的比對信息，提取對應gap 區(qū)域的序列，并將這些序列分為3 種不同的類型；4）基于多序列比對方法，采取一種新的策略將一個gap 區(qū)域關聯(lián)的3 種不同類型的序列進行一致化，獲得最終的填充序列。

2.1 初始比對結果

由于scaffold 中的gap 區(qū)域是由“N”組成的，并且可能比較長，因此直接把長讀數(shù)比對到scaffold上，可能會使一些長讀數(shù)無法比對到gap 的臨近區(qū)域。因此，首先刪除每條scaffold 中的gap 區(qū)域，這樣一條scaffold 變成多條contig，每條contig 是不包含gap 的長序列。然后利用比對方法BWA 將長讀數(shù)比對到contig 集合上。

2.2 比對結果修正

經(jīng)過上一步驟，本文方法可以獲得一個比對結果，但是由于長讀數(shù)的測序錯誤率較高，并且基因組中往往存在一些重復區(qū)，因此現(xiàn)有的比對工具容易造成一些錯誤的比對信息，或者給出的比對區(qū)間位置和真實情況相比具有一些偏差。為了獲得更加準確的gap 區(qū)域對應的序列，本文方法對每一條比對信息進行判斷和修正。

針對一條長讀數(shù)和一個contig 的比對結果，本文方法抽取它們之間的比對位置、比對方向和比對質量信息，并將一條比對信息定義為一個六元組（Rs，Re，Cs，Ce，Ori，Q）。其中，長讀數(shù)上的［Rs，Re］區(qū)間和contig 上的［Cs，Ce］區(qū)間能夠比對上，比對方向為Ori（0 代表反向比對，1 代表正向比對），Q代表比對質量分數(shù)。如果Q低于一個閾值（默認為30），本文方法認為該比對信息不可信，它不會用于后續(xù)的步驟。對于Q高于閾值的比對信息，本文方法計算長讀數(shù)和contig 之間比對區(qū)間的4 個偏移大小，分別是［Rs，Re］區(qū)間的左右兩端偏移大小DRs、DRe和［Cs，Ce］區(qū)間的左右兩端偏移大小DCs、DCe。如果4 個偏移量中有1 個大于α（默認500），則認為該比對是錯誤的比對信息，本文方法忽略該比對信息。如果4 個偏移量均小于α，則本文方法對初始的2 個比對區(qū)間進行修正，得到更新后的比對區(qū)間，具體過程如算法1 所示。

算法1RevisePosition（Rs，Re，Cs，Ce，Ori，Len（R），Len（C），α）

比對區(qū)間信息示例如圖2 所示。一條長讀數(shù)能分別比對到2 條contig 上，其中contig 1 上的區(qū)間［Cs1，Ce1］和長讀數(shù)上的區(qū)間［Rs1，Re1］能夠正向比對上，contig 2 上的區(qū)間［Cs2，Ce2］和長讀數(shù)上的區(qū)間［Rs2，Re2］能夠正向比對上，如圖2（a）所示；利用算法1 對這2 個比對信息進行修正，得到最終的比對區(qū)間信息，如圖2（b）所示。

圖2 區(qū)間信息示例Fig.2 Sample of interval information

2.3 gap 區(qū)域對應的3 類序列集合

當對所有的比對信息進行修正后，獲得一個新的比對信息集合。針對1 個gap 和其左右相鄰的2 個contig，根據(jù)比對到這2 個contig 的長讀數(shù)位置以及方向信息，抽取長讀數(shù)中對應于該gap 的序列區(qū)域。有些長讀數(shù)可以跨過該gap 區(qū)域，并且可以同時比對到這2 個contig 上。而有些長讀數(shù)只能比對這2 個contig 中的1 個，且只能覆蓋gap 區(qū)域的左右兩側部分區(qū)域。因此，本文方法針對1 個gap，把可能覆蓋它的長讀數(shù)區(qū)域抽取出來，并分別保存到3 個不同的集合中：跨過序列集合（SpanSet），左側序列集合（LeftSet）和右側序列集合（RightSet）。

針對1 個gap，假設其左右相鄰的2 個contig 是Ci和Cj。本文方法首先確定能夠比對到Ci和Cj的所有長讀數(shù)，如果其中1 個長讀數(shù)能夠同時比對到Ci和Cj，這2 個比對 信息分別是Alignment1=（Rsi，Rei，Csi，Cei，Orii，Qi）和Alignment2=（Rsj，Rej，Csj，Cej，Orij，Qj），則本文方法用算法2 判斷該長讀數(shù)是否能夠連接Ci和Cj，并且判斷該長讀數(shù)在這2 個contig 之間的序列是否對應該gap 區(qū)域。如果算法2 返回false，表明該長讀數(shù)無法連接Ci和Cj。如果算法2 返回true，則表明該區(qū)間對應gap 區(qū)域，并把對應gap 區(qū)域的序列保存到SpanSet。其中抽取序列的長度不超過該gap 長度加β（默認200）。ReverseComplement是對一個序列進行反向互補操作。

算法2DetermineSpanSet（Alignment 1，Alignment 2，SpanSet，β）

當一個長讀數(shù)能比對到Ci，但不能比對到Cj時，則該長讀數(shù)可能覆蓋該gap 區(qū)域的左側區(qū)域，本文方法用算法3 判斷該長讀數(shù)是否能覆蓋gap 區(qū)域的左側，如果能，則抽取對應gap 區(qū)域的序列，并保存到LeftSet 集合中。

算法3DetermineLeftSet（Alignment 1，LeftSet，β）

當一個長讀數(shù)能比對到Cj，但不能比對到Ci時，則該長讀數(shù)可能覆蓋該gap 區(qū)域的右側區(qū)域，本文方法用算法4 抽取對應gap 區(qū)域的序列，并保存到RightSet 集合中。

算法4DetermineRightSet（Alignment 2，RightSet，β）

根據(jù)上述步驟，本文方法獲得對應1 個gap 的3 個不同類型的序列集合，注意其中一些序列集合可能為空。

2.4 gap 區(qū)域填充

在填充1 個gap 時，關鍵是對該gap 關聯(lián)的3 個集合中的序列進行一致化，獲得一致序列。racon［22］是一個快速的序列一致算法，在該方法中，首先要確定一個目標序列（targe sequence），然后利用minimap2［23］將一個序列集合比對到該目標序列，并獲得它們之間的重疊信息（overlap），最后利用該重疊信息對目標序列進行糾錯達到序列一致化的目的。但是當目標序列和序列集合中的序列長度較短時，minimap2 往往無法獲得這些序列的重疊信息，從而無法獲得一致序列。針對racon 的缺陷，本文方法在racon 無輸出結果時，直接利用多序列比對方法mafft［24］獲得多條序列的比對信息，然后設計一種一致化策略獲得一致序列。本文方法基于序列一致算法racon 和多序列比對方法mafft，采用一種新的策略對一個gap 進行填充。

針對一 個gap 和其對應的SpanSet、LeftSet 和RightSet 集合進行填充，具體的填充步驟如下：

1）如果SpanSet 中的序列個數(shù)大于等于3，利用racon 嘗試求出一個一致序列。在此過程中，首先利用minimap2 對這些SpanSet 中的序列進行重疊檢測，統(tǒng)計每條序列有多少條其他序列能夠比對到它本身，并選擇獲得最多比對的序列作為目標序列，如果同時有多條，則選擇和gap 長度最接近的序列作為目標序列。然后以該目標序列、重疊信息和其他序列作為racon 的輸入，racon 的輸出作為一致序列。如果racon 沒有輸出任何結果，則利用mafft 對SpanSet 中的所有序列進行比對，并獲得這些序列的布局（layout），最后基于少數(shù)服從多數(shù)的原則，選擇每個堿基位置獲得最多支持的堿基，并獲得一條一致序列。通過上述步驟獲得的一致序列用SpanConsensusSequence 表示。

2）如果SpanSet 中的序列個數(shù)小于3 大于0，則首先從SpanSet 中找到一個目標序列（方法和上一步找目標序列的方法一樣），然后利用minimap2 把3 個序列集合中剩余的序列比對到目標序列上，接著用racon 進行一致化。如果racon 沒有輸出，則利用mafft 對SpanSet 和LeftSet 中的所有序列進行多序列比對，求出一個一致序列。然后再用mafft 對SpanSet 和RightSet 中的所有序列進行多序列比對，求出一個一致序列，最后對這2 條一致序列進行重疊檢測，并合并這2 個序列求出一條新的一致序列。如果LeftSet 或者RightSet 為空集，則直接把目標序列當作一致序列。通過上面步驟獲得的一致序列用SpanConsensusSequence 表示。

3）如果SpanSet 中的序列個數(shù)為0，即為空集，則利用racon 分別對SpanSet 和RightSet 中的序列進行一致化求解，分別求出一個左側一致序列和一個右側一致序列（求解過程和上述的方法相同）。如果利用racon 無法求出其左側或者右側一致序列，則利用mafft 分別對LeftSet 和RightSet 中的序列進行多序列比對，并分別求出左側和右側一致序列，這2 條一致序列分別用LeftConsensusSequence 和RightConsensus Sequence表示，注意如果LeftSet或者RightSet 為空集，則對應的一致序列為空。

4）如果對于一個gap，通過上述步驟可以獲得SpanConsensusSequence，則直接利用該序列填充這個gap。如果對于一個gap只有LeftConsensus Sequence，則根據(jù)gap 的長度，從左往右填充這個gap，最多填充gap 長度。如果對于一個gap 只有RightConsensusSequence 序列，則根據(jù)gap 的長度，從右往左填充這個gap，最多填充gap 長度。如果對于一個gap 既有LeftConsensusSequence，又有Right ConsensusSequence，則同時從左往右，和從右往左填充gap，最多填充gap 長度。

通過上述步驟，獲得一個gap 的最終填充序列。當依次處理完所有的gap 后，則完成整個gap 填充過程。

3 時間復雜度

本文方法包括比對階段、比對修正階段、確定gap 對應的3 個序列集合階段和gap 填充階段4 個階段。

1）比對階段。本文方法調用BWA-MEM 算法將長讀數(shù)比對到scaffold 集合上，該階段的時間復雜度為O（w）+O（s），其中，w是所有長讀數(shù)長度之和，s是所有scaffold 長度之和。

2）對比對修正階段。本文方法依次調用算法1處理所有的比對信息，其時間復雜度是O（m），m是比對信息個數(shù)。

3）確定gap 對應的3 個序列集合階段。針對每一個gap，本文方法確定所有能夠比對其左右兩側contig 的比對信息。如果一個長讀數(shù)能跨過該gap，則調用算法2；如果一個長讀數(shù)只比對到該gap 的左側contig，則調用算法3；如果一個長讀數(shù)只比對到該gap 的右側contig，則調用算法4。這個步驟的時間復雜度是O（t×m），其中t是gap 的個數(shù)。

4）填充gap 階段。針對一個gap 區(qū)域及其關聯(lián)的3 個序列集合，利用多序列比對算法racon 和mafft產(chǎn)生一致序列。racon 時間復雜度是O（p×L），mafft時間復雜度是O（p2×L）+O（p×L2），p是進行多序列比對的序列個數(shù)，L是序列的長度。

4 實驗與結果分析

4.1 實驗數(shù)據(jù)

為驗證gapLM 的有效性，將gapLM 在2 個真實Pacibio 測序數(shù)據(jù)集上和其他3 種常見的gap 填充方法（GMcloser，PBJelly，LR_Gapcloser）進行比較。這2 個數(shù)據(jù)集分別來自2 個物種：釀酒酵母S288C（S.cerevisiae）和秀麗隱桿線蟲（C.elegans）。本文方法采用2 個由NGS 組裝產(chǎn)生的scaffold 集合作為原始輸入序列，使用不同的gap 填充方法對這些原始scaffold 集合進行填充。這2 個長讀數(shù)集合和scaffold 集合的詳細信息如表1 所示。數(shù)據(jù)集都來自于文獻［20］。

表1 數(shù)據(jù)集信息Table 1 Datasets information

4.2 評價指標

本文使用QUAST［25］對最終的填充結果進行評估。QUAST 將參考基因組與gap 填充后的scaffold集合進行比對，并獲得多個評價指標。組裝結果的好壞主要取決于長度和準確度這2 個方面。本文利用QUAST 進行相關評價時，主要列出以下4 項評測指標：

1）N50。contig 或scaffold 按照長度大小進行排序，并依次將其長度相加，當累加值超過所有contig（或者scaffold）長度之和的1/2 時，最后一個contig 或scaffold 的長度就是N50［26］的值。

2）NGA50。在組裝結果發(fā)生錯誤的位置斷開contig 或scaffold，然后重新計算N50 值。

3）基因組覆蓋比例（GF）。將組裝好的scaffold或者contig 與參考基因組進行比對后，計算覆蓋到參考基因組的堿基百分比。

4）組裝錯誤（MA）。在contig 或scaffold 與參考基因組進行比對后，發(fā)生在contig 或scaffold 中錯誤的個數(shù)。

4.3 結果分析

當gapLM、PBJelly、GMcloser 和LR_gapcloser分別對scaffold 集合進行填充后，為了更好地評價填充結果，本文將填充后的scaffold 在gap 區(qū)域處打斷，產(chǎn)生不包含gap 區(qū)域的contig 集合。這樣填充越多越正確的gap 填充方法產(chǎn)生的contig 數(shù)目也就越少，并且準確性也高。然后利用評價工具QUAST 將這些contig 集合和基因組參考序列進行比對，并獲得相應的評價指標。

4 種gap 填充方法最終結果的評價如表2 所示，對于S.cerevisiae 數(shù)據(jù)集，其具有相對較小的基因組。與參考數(shù)據(jù)集相比，原始的輸入數(shù)據(jù)集包含29 個錯誤組裝（MA）。gapLM 對可能覆蓋gap 區(qū)域的序列進行劃分，并分別抽取有關序列進行多序列比對產(chǎn)生一致序列。在第1 組數(shù)據(jù)集上，與其他方法相比，gapLM 獲得最大的NGA50值。雖然MA值比GMcloser 要高，但是GMcloser 的GF 值是最低的，填充效果不明顯。而對較大的C.elegans 數(shù)據(jù)集，本文方法獲得的gap 填充結果的N50 和NGA50 值都是最優(yōu)的，并且本文方法的MA 值在4 種方法中排名第2，這說明本文方法能夠獲得較好的填充結果。

表2 gap 填充評價結果Table 2 Evaluation results of gap filling

在運行時間上，所用gap 填充方法主要消耗在比對階段，這和選擇的比對工具相關。PBjelly 使用BLASR 進行比對，GMcloser 使用MUMcer 進行比對，LR_Gapcloser 和gapLM 都使用BWA-MEM進行比對。但是，LR_gapcloser 把長讀數(shù)打斷成一些小片段，然后將這些小片段比對到scaffold 集合上，這樣減少耗費在比對上的時間。gap 填充評價結果如表2 所示，其中，contigs≥1 000 bp。從表2 可以看出，LR_Gapcloser 在2 個數(shù)據(jù)集上都運行最快，并且消耗的內存最小，而GapLM 在運行時間和最大內存消耗上比GMcloser 和PBjelly 更優(yōu)。

與其他3 個常見的gap 填充方法相比，gapLM 的N50 和NGA50 值相對較高，并且gapLM 的錯誤率也較低，這表明填充區(qū)域序列的準確度比其他方法具有一定優(yōu)勢。因此，gapLM 產(chǎn)生更加準確和有效的填充結果。雖然GMcloser 方法利用長讀數(shù)進行填充，但是它填充后的scaffold 集合的覆蓋度和原始scaffold 集合的覆蓋度一樣，說明其沒有進行有效填充。PBjelly 方法能夠利用長讀數(shù)對gap 進行填充，但是其錯誤率顯著增加。LR_gapcloser 方法在填充結果上取得不錯的效果并且其運行效率最優(yōu)，但是由于其直接利用長讀數(shù)上的序列去填充對應的gap區(qū)域，因此其結果仍然存在一定的錯誤率，特別是在大基因組數(shù)據(jù)上，其填充結果的N50 值較低。

5 結束語

本文提出一種基于長讀數(shù)和多序列比對的gap填充方法GapLM。將scaffold 集合中的每條contig在gap 區(qū)域分割成不包含gap 的contig 集合，然后把長讀數(shù)比對到contig 上，并對比對結果進行修正以獲得更加準確的比對信息?；诒葘ξ恢煤头较蛐畔ⅲ@得可能覆蓋一個gap 的3 個序列集合。利用多序列比對方法和一種新的策略，對這3 個序列集合進行一致化，獲得一致序列以填充gap 區(qū)域。將gapLM 與目前常見的3 種其他gap 填充方法在2 個真實數(shù)據(jù)集上進行比較，實驗結果表明，gapLM 能夠產(chǎn)生更加準確的填充結果。但該方法仍然存在運行時間稍長等不足，因此設計一種快速的針對gap 填充的比對方法是下一步的工作重點。另外，為提高gap填充的完整性和準確性，充分挖掘長讀數(shù)和雙端短讀數(shù)的互補優(yōu)勢，同時利用長讀數(shù)和雙端短讀數(shù)進行gap 填充也是未來研究的方向。