王 昊，陳 挺,2*

(1.清華大學 計算機科學與技術(shù)系，北京 100084；2.清華大學 人工智能研究院，北京 100084)

DNA作為生物重要的遺傳物質(zhì)，一直被生物學家在不同層面進行深入研究。從發(fā)現(xiàn)DNA的結(jié)構(gòu)開始，人類通過研究復(fù)雜多樣的基因組，在健康和疾病的研究中已取得了長足的進步。為了研究DNA序列，DNA測序是必不可少的環(huán)節(jié)。由于現(xiàn)有技術(shù)僅能測得較短的DNA序列，因此將DNA測序獲得的讀段(Read)進行從頭拼接(De novo assembly)成為較長的重疊群(Contig)、進而拼裝成更長的骨架(Scaffold)、最終獲得完整的基因組(Genome)一直是生物信息學的重要課題。而隨著測序技術(shù)的發(fā)展，測序試劑和儀器不斷更新?lián)Q代，測序的速度、價格，以及測序獲得的讀段在質(zhì)量、長度、堿基偏好等方面特性也有變化，針對新的測序數(shù)據(jù)設(shè)計更適應(yīng)數(shù)據(jù)的拼裝算法也是不斷推陳出新。

1977年，由Sanger等人提出的鏈終止法測序方法開創(chuàng)了測序技術(shù)的先河[1]，這種測序方法后來通常被稱為Sanger測序法，即第一代測序技術(shù)。第一代測序技術(shù)的主要特點是獲得的讀段序列長度較長，通?？蛇_1 000 bp，準確性高，可達99.999%，但是測序的成本比較高，通量也比較低。這些特點決定了第一代測序技術(shù)通常被應(yīng)用于長度較短，或者十分重要的基因組測序任務(wù)中。

為了解決第一代測序技術(shù)具有通量低、成本高的局限性， 454焦磷酸測序[2]，Illumina[3]等測序儀器相繼誕生。這些新一代測序儀可以獲得每次運行上百萬上億次讀數(shù)的更高輸出，但是讀取長度更短，最新一代的Illumina機器現(xiàn)在輸出大約250-300 bp的高質(zhì)量讀數(shù)。這些測序技術(shù)現(xiàn)在被稱為第二代測序技術(shù)。第二代測序技術(shù)還包括Ion Torrents的Ion Proton測序儀。與Sanger，454和Illumina相比，這種技術(shù)不依賴于光學方法，可以使用購買成本較低的機器進行快速而便宜的測序。盡管第二代測序儀具有高通量、相對便宜的價格，同時具有高質(zhì)量的堿基和低錯誤率，但它們測得的讀段較短是主要問題?，F(xiàn)在，第二代測序技術(shù)可以測得包括復(fù)雜的哺乳動物在內(nèi)的絕大多數(shù)生物的基因組，并保證成本控制在可接受的范圍內(nèi)。這意味著幾乎所有生物的基因組都向人類打開了解碼的道路，并允許人們對這些基因組做進一步的分析。

隨著技術(shù)的發(fā)展，第三代測序技術(shù)應(yīng)運而生。由Pacific Biosciences(PacBio)公司開發(fā)的測序儀(如RS I，RS II和Sequel機器)使用單分子實時(SMRT)測序技術(shù)[2]，而牛津納米孔技術(shù)(ONT)開發(fā)了用于納米孔測序的裝置[4]。 與第二代技術(shù)(如454，Illumina，Ion Torrent)相比，這些第三代測序方法在測序文庫制備過程中不包括擴增步驟，因此能夠進行單分子測序，去除了擴增過程帶來的偏倚(Bias)。此外，第三代測序技術(shù)預(yù)期的讀段長度遠遠高于第二代技術(shù)，平均讀段長度可以超過6-8 kbp，最大讀段長度超過30-150 kbp[5]。使用第三代測序長讀段文庫可以定向、排序、間隔或連接基因組拼接中的重疊群，以提升這些基因組草圖的質(zhì)量。PacBio的SMRT長讀段已經(jīng)證明在解析長重復(fù)區(qū)域方面是有效的[6]，并且可能成為原核生物基因組測序的金標準。此外，SMRT長讀段已被用于解決黑猩猩[7]和人類[8]的復(fù)雜基因組區(qū)域。但是第三代測序儀的比較高的錯誤率決定了它們還不能完全取代第二代測序測序平臺。對于許多組學問題研究，比如宏基因組學，第二代技術(shù)至少在未來幾年仍將是最先進的技術(shù)[9]。

盡管從頭拼接算法取得了長足的進展，但短讀段或長讀段技術(shù)的固有問題仍然阻礙了完整的基因組的構(gòu)建。一方面，由于第二代測序數(shù)據(jù)誤差小，組裝第二代測序數(shù)據(jù)集可以產(chǎn)生準確的重疊群，但由于它們讀段尺寸太小而無法識別更大的基因組重復(fù)序列。另一方面，由于第三代測序數(shù)據(jù)的長度優(yōu)勢，純第三代測序數(shù)據(jù)裝配算法(例如HGAP[7])可以輕松地解決較大的重復(fù)區(qū)域，但是為了最小化錯誤率的影響以得到較高的準確率，需要巨大的測序深度。自第一個SMRT測序平臺發(fā)布，就出現(xiàn)了利用這兩種測序方法的互補性產(chǎn)生各種組裝混合數(shù)據(jù)集的想法以利用兩者的混合拼接算法。

1 現(xiàn)有拼接算法

盡管測序手段不斷升級，但無論是二代測序數(shù)據(jù)拼接還是二三代混合拼接，對測序數(shù)據(jù)進行拼接的算法思路是整體一致的。由于基因組中絕大部分的區(qū)域是沒有互相重復(fù)的，所以通過比對讀段上的堿基序列，找到他們的重疊區(qū)域(Overlap)，然后通過將這些讀段通過重疊區(qū)域連接起來以得到更長的序列，即重疊群(Contig)。由于拼接得到的重疊群數(shù)量通常遠多于實際的染色體數(shù)量，因此需要第二階段的拼接骨架以將部分碎片拼接為更長的序列。拼接骨架(Scaffolding)是將重疊群按照正確的順序和方向進行集合并進行連接的過程，其中間可能存在一些沒有具體測序數(shù)據(jù)間隙(Gap)以N(即未知堿基)進行填充，最終得到的序列被稱為骨架(Scaffold)。通過包括對特定的一些間隙進行填充(Gap fill)等方式，對拼接得到的進一步進行基因組后處理(Genome finishing)，即可以得到相應(yīng)物種的基因組草圖。

無論是將讀段連接為重疊群，還是進一步的骨架拼接，都可以歸納為如下數(shù)學模型：在讀段/重疊群間重疊信息的約束條件下，確認所有讀段/重疊群的順序和方向。可以證明，確定所有讀段/重疊群的順序等價于最優(yōu)線性排列問題，這一問題是NP難的( NP-hard)[10]。因此，尋求最優(yōu)解在計算復(fù)雜度上目前而言是不可行的，拼接算法需要尋求某種近似解，以在合理的時間和空間復(fù)雜度下獲得準確可以率接受的解?，F(xiàn)有的拼接算法大體可以按以下分類：

1.1 貪心法

貪心法[11]是最早采用于基因拼接算法的方法(見圖1)。貪心法的策略是選取某一問題的初始解，然后選擇最小代價的步驟，轉(zhuǎn)移到下一個逼近最終目標的解，如此反復(fù)以達到最終目標或者其他終止條件。在基因組拼接問題中，貪心法的具體步驟如下：首先選取一個讀段作為初始重疊群，然后尋找與當前重疊群重疊區(qū)域相似度得分最高的讀段，與當前重疊群進行拼接，得到新的重疊群，直到?jīng)]有讀段與當前重疊群足夠相似。思路簡單、復(fù)雜度低是貪心法的主要優(yōu)勢，但該方法只能找到局部最優(yōu)解。如何從局部最優(yōu)逼近全局最優(yōu)是該方法面臨的主要挑戰(zhàn)。使用這一方法的軟件有SSPACE等。

圖1 貪心法圖解Fig.1 An example of the greedy method

1.2 OLC算法

OLC算法(Overlap-Layout-Consensus)[12]如其名字所示，主要可分為以下三個步驟：(1)Overlap：對所有讀段進行兩兩比對，找到讀段間的重疊信息。但由于兩兩比對的復(fù)雜度太高，實際算法通常先采取某種方式篩選可能比較相似的讀段對，僅比較這些可能相似的讀段，其余讀段對相似度視為0。(2)Layout：根據(jù)得到的重疊信息將測序讀段視為圖上的點，讀段間的相似關(guān)系視為連線，構(gòu)建OLC圖(見圖2)。(3)Consensus：利用OLC圖上的連接信息，找到遍歷圖上所有節(jié)點的最佳路徑，即尋找圖上的哈密頓道路，這也導(dǎo)致該算法復(fù)雜度通常較高，通常用于讀段較長、總讀段數(shù)較少的拼接。使用這一方法的軟件有Celera Assembler等。

圖2 OLC法圖解Fig.2 An example of the OLC method

1.3 De Bruijn圖法

De Bruijn圖[13]是另外一種被許多流行的拼接軟件使用的圖模型算法，該算法首先被拼接軟件EULER提出。該算法首先將每個讀段分解為一個k-mer(此處使用的長度k遠遠小于OLC法中使用的長度)，并以這些k-mer為邊，(k-1)-mer為節(jié)點構(gòu)建一張圖。如果兩個(k-1)-mer 在k-mer上連續(xù)出現(xiàn)，則他們在圖上會連接一條邊，每條讀段均會形成一條鏈，將這些鏈上的點和邊的信息整合到一張圖上，就得到了De Bruijn圖(見圖3)。這種方法沒有顯式地對讀段間的重疊信息進行計算，而是通過k-mer建圖的方式進行有效的捕獲，從而節(jié)省了大量的計算時間。接下來，算法通過遍歷DBG中的無分支路徑的方式來構(gòu)造重疊群。然而，由于測序錯誤，或者存在重復(fù)序列，或者基因組的多樣性如多條染色體，都會導(dǎo)致高度相似的序列中存在局部差異，這些差異在圖中就表現(xiàn)為一些“分叉”和“氣泡”，導(dǎo)致得到的高質(zhì)量的重疊群非常短。為了解決這一問題，程序會預(yù)先在圖中搜索分叉和氣泡并檢測對應(yīng)讀段是否可能是錯誤序列，盡可能地剪除分叉和氣泡。如果路徑對應(yīng)的讀段非常少，則說明該路徑很可能是測序錯誤導(dǎo)致的，因而應(yīng)該刪去。實踐證明，如果在實際裝配過程前加入錯誤校正步驟，可以得到更好的拼接結(jié)果。該方法的主要缺點是在將讀段轉(zhuǎn)化為k-mer路徑的過程中，大量k-mer間的信息丟失，導(dǎo)致圖上路徑的正確性很容易受到其他k-mer連入邊的干擾。同時，由于該方法使用了大量基于圖的計算，導(dǎo)致難以進行并行化處理，因而在處理哺乳動物等復(fù)雜生物的基因組測序數(shù)據(jù)時不能有效拓展。除此之外，圖的存儲需要極其大量的內(nèi)存，因而限制了該方法的應(yīng)用范圍。使用這一方法的軟件有SPAdes等。

圖3 De Bruijn圖法圖解Fig.3 An example of the De Bruijn graph method

注：圖中三條讀段作為原始數(shù)據(jù)，分別先以長度為3的k-mer為節(jié)點進行分割、再用長度為4的k-mer為邊連接節(jié)點，形成一條讀段鏈；然后整合三條讀段鏈上k-mer的連接信息生成De Bruijn圖的示例。三條讀段鏈以折線標注在圖中。圖中右下的方框內(nèi)為一個“氣泡”，左上的方框內(nèi)為一個“分叉”。

1.4 其他方法

除了以上基本方法之外，還有一些工具通過結(jié)合以上方法的優(yōu)勢，提出組合式的新算法。例如Pacific Biosciences公司開發(fā)的軟件FALCON提出的序列圖(String graph)算法[14]。該方法首先使用OLC中的方法計算讀段間的重疊信息，然后根據(jù)重疊區(qū)域的邊緣將序列切割成長度不等子序列，并以這些子序列為節(jié)點，在同一讀段中的連接順序為邊建立類似DBG的圖結(jié)構(gòu)。然后在該圖上遍歷無分支路徑獲得重疊群。

2 混合拼接軟件介紹

自第三代測序問世以來，已經(jīng)有許多混合拼接工具問世。這些工具除了在拼接算法主要思路的選取上有所不同，如何結(jié)合使用二、三代測序數(shù)據(jù)的方式也各有異同。大體而言，可以主要分為以下幾種：

2.1 三代骨架拼接

用二代數(shù)據(jù)拼接為重疊群，然后結(jié)合三代數(shù)據(jù)進行骨架拼接(見圖4)。這一方式的主要創(chuàng)新點在骨架拼接階段。該方式的優(yōu)點在于可以利用之前二代測序單獨拼接、骨架拼接的相關(guān)成果。但是如果由于重復(fù)區(qū)域等原因，二代拼接結(jié)果有較多錯誤，難以利用三代測序數(shù)據(jù)對二代拼接結(jié)果進行直接校正，這是該方法的不足。采用這種方式的軟件有：AHA[15]，Cerulean[16]和SSPACE-LongRead[17]等。

2.2 校正-拼接

將二代測序數(shù)據(jù)比對到三代測序數(shù)據(jù)上，對三代測序數(shù)據(jù)進行校正，再使用校正后的三代數(shù)據(jù)進行拼接(見圖5)。這種方法的創(chuàng)新點主要在比對和校正三代測序數(shù)據(jù)上。校正后的流程與前代拼接流程需求基本相同，但由于讀段長度變化極大，拼接方法的選擇相對比較受限，Celera Assembler (runCA)[18]是通常的選擇。采用這種方式的軟件有：PBcR pipeline[19]，此外還有專門的糾錯軟件Proovread[20]。

圖4 三代骨架拼接Fig.4 Third-generation sequencing scaffolding

圖5 校正-拼接Fig.5 Correction-assembly

2.3 拼接-校正-拼接

將二代測序數(shù)據(jù)進行預(yù)處理，再將重疊群比對到三代測序數(shù)據(jù)上，對三代測序數(shù)據(jù)進行校正，再使用校正后的三代數(shù)據(jù)進行拼接。由于第二代測序的讀段較短，而第三代讀段錯誤率較高，直接將第二代測序數(shù)據(jù)比對到第三代測序數(shù)據(jù)通常會得到較多可能的正確位置，直接選出正確位置是一個相當具有挑戰(zhàn)性的問題。為了優(yōu)化這一問題，先對第二代測序數(shù)據(jù)進行預(yù)處理，如拼接成重疊群、構(gòu)建de-Bruijn圖、將相鄰的相同堿基壓縮等方式。以拼接成重疊群為例，這相當于增加了讀段的長度，可以減少可能比對到的位置，提高了比對和校正的準確性。但是如何矯正可能的二代數(shù)據(jù)錯誤拼接是這一方法新引入的問題。ECTools[21]采取的方式是將二代數(shù)據(jù)拼接為重疊群，LoRDEC[22]采取的方式是先用二代測序數(shù)據(jù)構(gòu)建de-Bruijn圖，LSC[23]采取的方式是將相鄰的相同堿基壓縮為一個以應(yīng)對三代數(shù)據(jù)插入、刪除錯誤較多的問題。拼接流程圖見圖6。

圖6 拼接-校正-拼接Fig.6 Assembly-correction-assembly

2.4 完全混合拼接

即直接將第二代、第三代測序數(shù)據(jù)混合輸入，設(shè)計綜合使用兩種數(shù)據(jù)的拼接算法。這種方案理論上可以更好地應(yīng)用兩種數(shù)據(jù)的優(yōu)勢，但算法通常復(fù)雜。采取這種方式的有SPAdes-Hybrid[24]的混合拼接，其主要思路為先用二代測序數(shù)據(jù)構(gòu)建de-Bruijn圖，然后利用三代測序數(shù)據(jù)在圖上的比對結(jié)果計算得到重疊群。

3 性能比較

為了比較不同混合拼接軟件的拼接效果、拼接性能等指標，Hsin-Hung Lin等人[25]使用Escherichia coli K12 MG1655 (E. coli), 151-bp paired end 二代真實測序數(shù)據(jù)，S. Koren 等人研究中發(fā)布的第三代SMRT測序數(shù)據(jù)數(shù)據(jù)(SMRT1)[19]，Pacific Biosciences’ DevNet網(wǎng)站(http://pacificbiosciences.github.io/DevNet/)發(fā)布的第三代SMRT E. coli測序數(shù)據(jù)(SMRT2) 進行測試。其中，SMRT1讀段平均長度為2.0 Kbp，測序深度為16X；SMRT2讀段的平均長度為2.6 Kbp，測序深度23X。測試中，二代測序數(shù)據(jù)拼接使用軟件Abyss[26]先進行二代數(shù)據(jù)單獨從頭拼接，為需要二代測序數(shù)據(jù)拼接的軟件提供支持。拼接結(jié)果使用QUAST[27]軟件與參考基因組 (NCBI reference sequence NC000913.2)進行評估，該軟件是被廣泛使用對組裝結(jié)果的質(zhì)量進行評估的軟件。評估的主要指標為NGA50。NGA50是將所有拼接得到的重疊群在其錯誤拼接處打斷，然后統(tǒng)計這些打斷后的重疊群按長度從大到小排列，然后依次相加，當相加結(jié)果恰好超過參考基因組總長度的50%時，取該重疊群的長度作為NGA50的值。使用該指標是評估拼接結(jié)果準確性和連續(xù)性的指標。對該文章中的測試結(jié)果進行整合后，繪制出圖7。 其中部分軟件運行結(jié)果給出了更詳盡的數(shù)據(jù)，整理后列在表1中。

圖7 以NGA50度量的不同拼接軟件結(jié)果對比Fig.7 Comparison of NGA50 results of different assembly software tools

從圖7中，可以看出，在本次實驗中，AHA，Cerulean和SSPACE-LongRead這三款三代骨架拼接的軟件中，AHA和SSPACE-LongRead的效果穩(wěn)定，受二代拼接結(jié)果和三代測序讀段情況影響較小，但是結(jié)果最差；而Cerulean拼接結(jié)果受二代拼接結(jié)果和三代測序讀段情況的影響較大，拼接性能較為不穩(wěn)定，在測序深度和平均長度較小的SMRT1的數(shù)據(jù)中表現(xiàn)較好，這可能是由于在骨架拼接中，長讀段的測序深度不需要太高，過高的測序深度反而引入更多的噪音影響軟件性能。校正-拼接軟件PBcR pipeline受三代測序讀段的情況影響非常大，結(jié)果很不穩(wěn)定，好的結(jié)果超過前三者，但是差的結(jié)果甚至不如二代數(shù)據(jù)單獨拼接。PBcR在三代測序深度更高的SMRT2上效果更好，這可能是由于其需要使用矯正后的三代序列進行單獨拼接，因此更高的三代測序深度可以得到更好的結(jié)果。運行時間相對較短是該方法的優(yōu)點。ECTools + runCA等拼接-校正-拼接的軟件組合，都可以得到較好的結(jié)果，相較之下ECTools + runCA和Proovread + runCA的結(jié)果好于其他兩個組合。此外，這ECTools + runCA也在覆蓋深度較高的SMRT2數(shù)據(jù)上結(jié)果較好，可以說明長讀段測序深度的增加這類軟件組合最終拼接質(zhì)量。相比而言，這類算法運行時間整體偏長。

表1 不同長讀段校正工具參與的混合拼接結(jié)果比較Table 1 Comparison of hybrid assembly results with different assembly strategy softwares

為了對完全混合拼接方法的軟件Hybrid-SPAdes進行性能測試，本文在大腸桿菌數(shù)據(jù)集(E.colistr.K12 )上對其進行了測試。使用的大腸桿菌數(shù)據(jù)集包括：2 代 WGS 配對數(shù)據(jù)(測序儀器Illumina HiSeq 2000)，讀段長度 150 bp，數(shù)據(jù)總量 37.5 Gbp，下載自 NCBI SRA 數(shù)據(jù)庫，編號 SRR8365224；3 代 SMRT WGS 數(shù)據(jù)(測序儀器 PacBio RS II)，讀段平均長度 1 027 bp，測序深度22X，數(shù)據(jù)總量 1.2 Gbp，下載自 NCBI SRA 數(shù)據(jù)庫，編號SRR2063112。此外本文還使用了PBSIM[]模擬器生成了不同深度的3 代 SMRT 測序數(shù)據(jù)以測試。拼接參數(shù)選擇默認參數(shù)(16線程并行)，最后拼接得到結(jié)果以QUAST軟件評估，具體結(jié)果如表二所示。從模擬數(shù)據(jù)的結(jié)果上，可以看到NGA50指標變化較小，說明Hybrid-SPAdes對三代測序深度的魯棒性較好，但綜合其他指標，還是可以看到拼接結(jié)果隨著測序深度增加而變好。另外真實數(shù)據(jù)上的結(jié)果遠遠好于同等測序深度下的拼接結(jié)果，關(guān)鍵指標NGA50比深度更深(40X)的模擬數(shù)據(jù)結(jié)果，說明現(xiàn)有模擬器生成的模擬數(shù)據(jù)與實際測序數(shù)據(jù)仍然有較大偏差，比較應(yīng)以真實數(shù)據(jù)為主。真實數(shù)據(jù)上，Hybrid-SPAdes的NGA50指標為492.2 kbp，拼接結(jié)果質(zhì)量較為一般，與PBcR pipeline的NGA50較為接近。內(nèi)存消耗方面，Hybrid-SPAdes內(nèi)存消耗相對較小且非常穩(wěn)定，但在運算時間方面有一定波動，但整體用時大約在3 h左右，相對較少。結(jié)果見表2。

表2 Hybrid-SPAdes拼接結(jié)果Table 2 Hybrid-SPAdes assembly results

4 總 結(jié)

本文分析了在二代、三代測序數(shù)據(jù)混合拼接領(lǐng)域的現(xiàn)有軟件，并比較了他們測試性能，希望為需要使用二代、三代測序數(shù)據(jù)進行相關(guān)生物學研究的團隊提供參考，方便選出更合適的拼接軟件。通過文中的比較，拼接-校正-拼接的軟件組合性能較好，可以優(yōu)先考慮使用，其中以ECTools + runCA和Proovread + runCA的組合最為推薦。此外，如果三代測序數(shù)據(jù)深度較低，三代骨架拼接的Abyss+Cerulean的軟件組合也是可以考慮嘗試的。如果三代測序深度較高質(zhì)量較好，則可以考慮嘗試校正-拼接的軟件PBcR pipeline。而無論三代數(shù)據(jù)測序深度如何，完全混合拼接軟件Hybrid-SPAdes都是值得嘗試的。不同混合拼接策略的比較整理在表3中。

表3 不同拼接策略混合拼接軟件比較Table 3 Comparison of hybrid assembly softwares with different assembly strategies

另外，對于致力于通過優(yōu)化拼接算法，得到更好的拼接結(jié)果的研究者，希望通過本文的分析，能幫助有志于優(yōu)化混合拼接算法的團隊提供參考和比較的方向，以便更好地設(shè)計出更加高效、高質(zhì)量的混合拼接軟件。同時，也為需要使用二代、三代測序數(shù)據(jù)進行相關(guān)生物學研究的團隊提供參考，方便選出更合適的拼接軟件。隨著測序技術(shù)的發(fā)展，測序數(shù)據(jù)通量更高、價格更便宜、測序深度更長，這些對基因組拼接軟件即是數(shù)據(jù)優(yōu)勢，也是計算挑戰(zhàn)。通過優(yōu)化算法、提高并行度等方式，提高計算效率、降低計算復(fù)雜度，應(yīng)該是未來一段時間二三代數(shù)據(jù)混合拼接算法的發(fā)展方向。