第二、三代基因組測序數據混合拼接軟件綜述

王 昊，陳 挺,2*

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

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

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

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

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

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

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

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

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

1.1 貪心法

貪心法[11]是最早采用于基因拼接算法的方法(見圖1)。貪心法的策略是選取某一問題的初始解，然后選擇最小代價的步驟，轉移到下一個逼近最終目標的解，如此反復以達到最終目標或者其他終止條件。在基因組拼接問題中，貪心法的具體步驟如下：首先選取一個讀段作為初始重疊群，然后尋找與當前重疊群重疊區(qū)域相似度得分最高的讀段，與當前重疊群進行拼接，得到新的重疊群，直到沒有讀段與當前重疊群足夠相似。思路簡單、復雜度低是貪心法的主要優(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：對所有讀段進行兩兩比對，找到讀段間的重疊信息。但由于兩兩比對的復雜度太高，實際算法通常先采取某種方式篩選可能比較相似的讀段對，僅比較這些可能相似的讀段，其余讀段對相似度視為0。(2)Layout：根據得到的重疊信息將測序讀段視為圖上的點，讀段間的相似關系視為連線，構建OLC圖(見圖2)。(3)Consensus：利用OLC圖上的連接信息，找到遍歷圖上所有節(jié)點的最佳路徑，即尋找圖上的哈密頓道路，這也導致該算法復雜度通常較高，通常用于讀段較長、總讀段數較少的拼接。使用這一方法的軟件有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é)點構建一張圖。如果兩個(k-1)-mer 在k-mer上連續(xù)出現(xiàn)，則他們在圖上會連接一條邊，每條讀段均會形成一條鏈，將這些鏈上的點和邊的信息整合到一張圖上，就得到了De Bruijn圖(見圖3)。這種方法沒有顯式地對讀段間的重疊信息進行計算，而是通過k-mer建圖的方式進行有效的捕獲，從而節(jié)省了大量的計算時間。接下來，算法通過遍歷DBG中的無分支路徑的方式來構造重疊群。然而，由于測序錯誤，或者存在重復序列，或者基因組的多樣性如多條染色體，都會導致高度相似的序列中存在局部差異，這些差異在圖中就表現(xiàn)為一些“分叉”和“氣泡”，導致得到的高質量的重疊群非常短。為了解決這一問題，程序會預先在圖中搜索分叉和氣泡并檢測對應讀段是否可能是錯誤序列，盡可能地剪除分叉和氣泡。如果路徑對應的讀段非常少，則說明該路徑很可能是測序錯誤導致的，因而應該刪去。實踐證明，如果在實際裝配過程前加入錯誤校正步驟，可以得到更好的拼接結果。該方法的主要缺點是在將讀段轉化為k-mer路徑的過程中，大量k-mer間的信息丟失，導致圖上路徑的正確性很容易受到其他k-mer連入邊的干擾。同時，由于該方法使用了大量基于圖的計算，導致難以進行并行化處理，因而在處理哺乳動物等復雜生物的基因組測序數據時不能有效拓展。除此之外，圖的存儲需要極其大量的內存，因而限制了該方法的應用范圍。使用這一方法的軟件有SPAdes等。

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

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

1.4 其他方法

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

2 混合拼接軟件介紹

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

2.1 三代骨架拼接

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

2.2 校正-拼接

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

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

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

2.3 拼接-校正-拼接

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

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

2.4 完全混合拼接

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

3 性能比較

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

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

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

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

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

表2 Hybrid-SPAdes拼接結果Table 2 Hybrid-SPAdes assembly results

4 總 結

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

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

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