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

（河南理工大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，河南焦作 454000）

0 概述

隨著人類基因組計(jì)劃的完成，基因組數(shù)據(jù)成為后基因組時(shí)代的重要標(biāo)志，這些海量基因組數(shù)據(jù)的分析和挖掘深入到生命科學(xué)研究的各個(gè)方面，是眾多生命科學(xué)研究的基礎(chǔ)，更是計(jì)算機(jī)科學(xué)特別是生物信息學(xué)研究所面臨的重大挑戰(zhàn)?！秶摇笆濉笨萍紕?chuàng)新規(guī)劃》和《國家自然科學(xué)基金委“十三五”發(fā)展規(guī)劃》都明確將基因組數(shù)據(jù)分析以及生物大數(shù)據(jù)分析作為以“重大計(jì)劃”形式倡導(dǎo)的研究內(nèi)容。

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

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

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

1 相關(guān)工作

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

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

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

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

2 本文方法

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

2.1 初始比對結(jié)果

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

2.2 比對結(jié)果修正

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

針對一條長讀數(shù)和一個(gè)contig 的比對結(jié)果，本文方法抽取它們之間的比對位置、比對方向和比對質(zhì)量信息，并將一條比對信息定義為一個(gè)六元組（Rs，Re，Cs，Ce，Ori，Q）。其中，長讀數(shù)上的［Rs，Re］區(qū)間和contig 上的［Cs，Ce］區(qū)間能夠比對上，比對方向?yàn)镺ri（0 代表反向比對，1 代表正向比對），Q代表比對質(zhì)量分?jǐn)?shù)。如果Q低于一個(gè)閾值（默認(rèn)為30），本文方法認(rèn)為該比對信息不可信，它不會用于后續(xù)的步驟。對于Q高于閾值的比對信息，本文方法計(jì)算長讀數(shù)和contig 之間比對區(qū)間的4 個(gè)偏移大小，分別是［Rs，Re］區(qū)間的左右兩端偏移大小DRs、DRe和［Cs，Ce］區(qū)間的左右兩端偏移大小DCs、DCe。如果4 個(gè)偏移量中有1 個(gè)大于α（默認(rèn)500），則認(rèn)為該比對是錯(cuò)誤的比對信息，本文方法忽略該比對信息。如果4 個(gè)偏移量均小于α，則本文方法對初始的2 個(gè)比對區(qū)間進(jìn)行修正，得到更新后的比對區(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 個(gè)比對信息進(jìn)行修正，得到最終的比對區(qū)間信息，如圖2（b）所示。

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

2.3 gap 區(qū)域?qū)?yīng)的3 類序列集合

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

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

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

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

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

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

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

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

2.4 gap 區(qū)域填充

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

針對一 個(gè)gap 和其對應(yīng)的SpanSet、LeftSet 和RightSet 集合進(jìn)行填充，具體的填充步驟如下：

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

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

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

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

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

3 時(shí)間復(fù)雜度

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

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

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

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

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

4 實(shí)驗(yàn)與結(jié)果分析

4.1 實(shí)驗(yàn)數(shù)據(jù)

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

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

4.2 評價(jià)指標(biāo)

本文使用QUAST［25］對最終的填充結(jié)果進(jìn)行評估。QUAST 將參考基因組與gap 填充后的scaffold集合進(jìn)行比對，并獲得多個(gè)評價(jià)指標(biāo)。組裝結(jié)果的好壞主要取決于長度和準(zhǔn)確度這2 個(gè)方面。本文利用QUAST 進(jìn)行相關(guān)評價(jià)時(shí)，主要列出以下4 項(xiàng)評測指標(biāo)：

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

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

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

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

4.3 結(jié)果分析

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

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

表2 gap 填充評價(jià)結(jié)果Table 2 Evaluation results of gap filling

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

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

5 結(jié)束語

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