SuperLLEC：全新的鏈讀和長讀測序組裝糾錯算法

崔雅軒，張少強

天津師范大學(xué) 計算機與信息工程學(xué)院，天津 300387

1977年誕生的第一代DNA（deoxyribo nucleic acid，脫氧核糖核酸）測序技術(shù)（Sanger測序）[1]直接推動了2000年人類基因組計劃的完成[2]。Sanger測序的讀長達(dá)到1 000 bp（base pairs，堿基對），精準(zhǔn)度達(dá)到99%，但該測序技術(shù)通量低且成本高昂[2]。為了降低成本和提高通量，第二代測序技術(shù)應(yīng)運而生，第二代測序技術(shù)主要包括ABI公司的SOLID測序技術(shù)，Roche/454的GS FLX和Illumina公司的Solexa Genome Analyzer[3]。第二代測序一般產(chǎn)生長為75~500 bp的短讀（short-reads）。為了在高通量下提高測序讀長，近幾年第三代長讀（long-reads）測序技術(shù)開始大量應(yīng)用，最具代表性是太平洋生物科學(xué)公司（PacBio）的單分子實時測序（single molecule real time，SMRT）[4]和牛津納米孔技術(shù)公司（Oxford Nanopore）的納米孔測序[5]。相較于前兩代測序技術(shù)，第三代測序具有無需PCR（聚合酶鏈?zhǔn)椒磻?yīng)，polymerase chain reaction）擴增和讀長更長（一般10 000~30 000 bp）的特點，但測序準(zhǔn)確度有所降低（準(zhǔn)確率在85%左右）[6-8]。

由于第三代測序較高的錯誤率，為了提高序列組裝的準(zhǔn)確度，長讀測序的基因組組裝算法往往需要與糾錯算法進行合作。目前基于第三代測序的基因組組裝算法主要有Wtdbg2[9]、FALCON[10]和Canu[11]等；這三個算法均能夠基于單體型進行組裝，識別其中的單核酸多態(tài)性（single nucleotide polymorphism，SNP）位點。而目前第三測序的糾錯算法主要分為兩類[12]：一類是第三代數(shù)據(jù)自糾，這類最新算法主要有FLAS[13]、LoRDEC[14]等；另一類是借助第二代數(shù)據(jù)對第三代數(shù)據(jù)進行糾錯，例如基于序列比對的糾錯算法proovread[15]和CoLoRMap[16]；以及基于構(gòu)建de Bruijn圖的糾錯算法FMLRC[17]和Par-LECH[18]等。而PacBioToCA算法既能對PacBio測序數(shù)據(jù)自糾，也能用二代Illumina測序數(shù)據(jù)對PacBio數(shù)據(jù)進行糾錯[19]。

上述糾錯算法均采用先糾錯再序列組裝的策略，但2017年10xGenomics公司對Illumina二代測序方案進行升級，引入條碼（barcode）標(biāo)簽序列，相同條碼的短讀集合形成跨度在30 000~100 000 bp的鏈讀（linked-reads）信息[20]。為了充分運用序列跨度大的鏈讀數(shù)據(jù)關(guān)聯(lián)信息，和基于最新的第三代測序組裝算法性能的大幅提升，本研究放棄先對長讀糾錯再組裝的策略，擬采用先組裝再糾錯和進一步組裝的策略：先對PacBio長讀組裝得到的支架序列（scaffolds），再利用帶條碼標(biāo)簽的鏈讀測序數(shù)據(jù)設(shè)計了對支架進行糾錯和進一步組裝的算法SuperLLEC，由此得到更高質(zhì)量的基因組單體型序列。

1 數(shù)據(jù)描述

1.1 10X Genomics鏈讀

10X Genomics公司通過在序列中引入條碼（barcode）標(biāo)簽序列[20]，即一段固定長度的堿基短序列，并將其分配到不同的油滴微粒（droplet）中，利用GemCode平臺對帶條碼標(biāo)簽的長片段序列進行擴增，最后將序列打碎成合適測序短片段進行二代測序。通過條碼標(biāo)簽序列信息追蹤來自每個長片段DNA模板的多個讀序（read），稱為鏈讀（linked-reads），從而獲得該長片段的大量短讀序列。為了方便描述，本實驗將具有相同條碼標(biāo)簽的所有鏈讀序列合稱為“鏈讀序列組”。同一鏈讀序列組的所有讀序均來自相同的染色體（跨度為30 000~100 000 bp）。此信息有助于對組裝后的支架序列根據(jù)鏈讀相關(guān)性進行進一步拼接。研究所用鏈讀數(shù)據(jù)均從10X Genomics官網(wǎng)下載（https：//support.10xgenomics.com/de-novo-assembly/datasets）。

1.2 PacBio長讀

PacBio長讀序列，是太平洋生物公司的單分子實時測序技術(shù)平臺產(chǎn)生的測序讀序[21]。長讀序列的單一測序長度一般大于10 000 bp[22]。實驗用到的PacBio數(shù)據(jù)集來自三種人類基因組樣本的測序數(shù)據(jù)，具體信息見表1，數(shù)據(jù)下載地址見表2。

表1 實驗所需PacBio長讀數(shù)據(jù)Table 1 PacBio long read data in experiment

表2 實驗所需PacBio長讀數(shù)據(jù)下載地址Table 2 Download site of PacBio long read data

2 方法和步驟

算法主要分三步：首先對PacBio三代測序數(shù)據(jù)用組裝算法組裝為支架序列；然后將相同條碼標(biāo)簽的鏈讀序列組快速分配給支架集合，再用序列比對工具將所有鏈讀一一比對到對應(yīng)的支架上，并根據(jù)鏈讀信息對支架繼續(xù)組裝；最后，統(tǒng)計多態(tài)位點的堿基頻率來判斷是否為測序錯誤堿基或者是SNP。如果是測序錯誤堿基，則根據(jù)位點頻率進行糾錯以提高基因組拼接質(zhì)量。其中在鏈讀與支架比對過程中，可以對鏈讀集合進行手動拆分多個子集后多核并行化快速處理（文中用8個進程進行實驗）。圖1是算法的研究方法流程示意圖。下面就算法各步驟進行詳細(xì)闡述。

圖1 研究方法流程示意圖Fig.1 Schematic flowchart of research method

2.1 組裝長讀數(shù)據(jù)構(gòu)建支架序列

對原始的PacBio測序數(shù)據(jù)，首先用第三代測序的基因組組裝算法將其拼接成支架序列集合。Wtdbg2[9]、FALCON[10]、Canu[11]是三種最新的第三代測序組裝算法。實驗分別嘗試使用這三種算法依次對表1的三組長讀數(shù)據(jù)集進行序列組裝。對組裝完成的數(shù)據(jù)進行比較分析后發(fā)現(xiàn)，Canu組裝后的人類基因組數(shù)據(jù)規(guī)模（約為3.42 GB）大于Wtdbg2和FALCON的數(shù)據(jù)規(guī)模（約為3 GB）。另外實驗發(fā)現(xiàn)Wtdbg2在三者中是運算最快的算法。例如對于Human CHM1的PacBio長讀數(shù)據(jù)，Canu的CPU用時為22 750 h，F(xiàn)ALCON用時為68 789 h，而Wtdbg2用時只有245 h。這是因為Wtdbg2算法采用Linux的Pthreads多線程處理技術(shù)來加快算法的執(zhí)行速度。由此，實驗采用Wtdbg2算法對實驗的PacBio長讀數(shù)據(jù)組裝成支架序列集合。在運行Wtdbg2時，除了記錄支架序列外，還記錄下堿基多態(tài)位點（Wtdbg2預(yù)測為SNP）的堿基頻次（比對到該位點的長讀堿基頻次）。

2.2 鏈讀序列組分配給支架

對每個支架序列S，以k為步長，將其拆分成長度為k的短序（稱為k-mer），若支架序列S長度為L，則可以得到L/k個k-mer；并將得到的所有k-mer存儲在哈希表中，哈希表的鍵為k-mer，鍵值為該k-mer所在的支架編號S。對每組鏈讀序列中的讀序，以長度k為窗口，以步長1遍歷所有的k-mer，并通過哈希表確定該組鏈讀序列的k-mer存在于支架的編號集合，據(jù)此將每組鏈讀序列集合分配給支架。如圖2示例所示，通過判斷鏈讀和支架是否共享k-mer，將四組鏈讀序列集合依次分配給兩條支架序列。其中Barcode編號為2的鏈讀組同時分配給了兩條支架。

圖2 糾錯組裝示意圖Fig.2 Flowchart of assembly and error correction

為了加快運算速度，實驗將所有鏈讀序列組平均分成8組，調(diào)用8個線程分別進行支架分配。

2.3 鏈讀序列比對到支架并進一步組裝

將分配到每個支架的每組鏈讀短序比對到該支架序列上。由于經(jīng)典對比工具BLAST[23]采用的是近似比對算法，適用于長序列之間的快速比對，不適用于短序列比對。實驗采用Bowtie2[24]進行序列比對。Bowtie2具有比對速度快和占用內(nèi)存低的優(yōu)點。

對每組具有相同條碼標(biāo)簽的鏈讀序列，用Bowtie2工具把它們比對到與之關(guān)聯(lián)的支架上。如果某組鏈讀序列之前被分配給多個支架，則根據(jù)該組鏈讀條碼信息和重疊關(guān)系，將這些支架進一步拼接成超級支架（superscaffold）序列。如圖2示例所示，A和B兩個支架序列因為Barcode 2而被拼接成一條超級支架序列。

為了加快運算速度，實驗將所有支架序列分成8組，保證每組之間沒有共享的鏈讀序列集，分別進行序列比對。同時，對于比對出現(xiàn)多態(tài)（或稱不保守）的位點，分別記錄鏈接讀序在該位點出現(xiàn)的堿基頻次。例如，如圖3示例所示，在標(biāo)出的第一列位置，堿基C在對比上的鏈讀中出現(xiàn)3次。

圖3 比對后的堿基頻率示意圖Fig.3 Example of base frequency statistics after alignment

2.4 糾錯和SNP預(yù)測

算法最后采用Fisher精確檢驗（Fisher’s exact test）來檢驗每個位點堿基是否有錯誤堿基或者SNP。Fisher精確檢驗是一種分析列聯(lián)表（contingency table）統(tǒng)計顯著性檢驗方法，用于檢驗兩個分類的關(guān)聯(lián)度。

如果某個位點只有一種堿基，則無需計算，說明該位點堿基信息一致；否則，在支架對應(yīng)位點上，統(tǒng)計比對后的鏈讀出現(xiàn)頻次最高的兩個堿基X和Y（如圖3示例，支架在堿基不一致的T/G位點，該位點鏈讀出現(xiàn)頻次最高的堿基是T和C），并記錄頻次為clong和dlong，在支架上這兩個堿基出現(xiàn)的頻次為alinked和blinked。注意如果在支架序列出現(xiàn)不是X或Y的第三種堿基Z，根據(jù)鏈讀的精確性，自動認(rèn)定堿基Z為測序錯誤的位點，不考慮Z出現(xiàn)的頻次。為了保證假設(shè)檢驗堿基樣本數(shù)量的一致性，分別計算相對頻次c=[100clong/(clong+dlong)]，d=[100dlong/(clong+dlong)]，a=[100alinked/(alinked+blinked)]，b=[100blinked/(alinked+blinked)]。如表3所示，為了比較長讀組裝和鏈讀比對的差異，構(gòu)造長讀和鏈讀堿基頻次的列聯(lián)表，然后采用Fisher精確假設(shè)檢驗。對一個位點，F(xiàn)isher精確檢驗的零假設(shè)為長讀和鏈讀在該位點的相對頻次沒有差別。運用超幾何分布計算Fisher精確檢驗的P值（P-value）：

表3 某位點在長讀和鏈讀堿基頻次列聯(lián)表Table 3 Contingency table of site located by long reads and linked-reads

為了快速計算P值，公式（1）中的階乘均可以采用公式（2）的斯特林（Stirling）公式來近似計算，由此得到取對數(shù)的組合數(shù)的近似計算公式（3），并通過公式（4）防止數(shù)據(jù)溢出。

對于某位點，如果計算的P值大于等于設(shè)定的臨界值0.05，則判定零假設(shè)成立，該位點在鏈讀中出現(xiàn)頻率最高的兩個堿基預(yù)測為SNP；如果P值小于臨界值，則判定零假設(shè)不成立，那該位點在鏈讀中出現(xiàn)頻率較小的堿基則判定為測序錯誤的堿基。

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

3.1 實驗環(huán)境

實驗所有數(shù)據(jù)均在CPU為Intel Xeon E2244G 3.80 GHz 8 Cores，內(nèi)存為512 GB的CentOS系統(tǒng)服務(wù)器上調(diào)試運行。

3.2 結(jié)果分析

針對三組人類基因組樣本PacBio測序數(shù)據(jù)集（Human NA24385、Human HG00733、Human CHM1分別記為Data1、Data2、Data3），分別用PacBioToCA和LoRDEC對這三組PacBio測序數(shù)據(jù)集進行自糾；另外也基于10X Genomics的鏈讀數(shù)據(jù)分別運行PacBioToCA和proovread對這三組PacBio數(shù)據(jù)集進行糾錯；最后對糾錯后的三組數(shù)據(jù)用Wtdbg2算法進行基因組組裝。

針對SuperLLEC算法，實驗是先用Wtdbg2算法分別對這三組PacBio測序數(shù)據(jù)組裝成支架序列（其運行時間和輸出支架總規(guī)模見表4），然后運用10X Genomics的鏈讀數(shù)據(jù)用SuperLLEC對支架進行糾錯和進一步組裝。

表4 Wtdbg2運行時間和輸出規(guī)模Table 4 Running times and result sizes of Wtdbg2 on three datasets

支架NG50[7]是在基因組組裝中衡量算法優(yōu)劣的重要指標(biāo)之一，其代表含義是算法所組裝得到的支架序列的中位長度。如表4所示，在沒有運用任何糾錯算法的情況下，Wtdbg2對三個數(shù)據(jù)集進行組裝，其支架NG50長度分別為25.8 MB、12.7 MB和15.4 MB。而如表5所示，在運用PacBioToCA、LoRDEC和proovread進行糾錯后，再用Wtdbg2組裝，如表5所示，這三個數(shù)據(jù)集的支架NG50增長不明顯，三個工具中最大的分別為26.1 MB、13.1 MB和16.2 MB。除了PacBio數(shù)據(jù)自糾外，這三種算法均用到了鏈讀數(shù)據(jù)，但是這三種算法只把鏈讀單純作為獨立的短讀，沒有用到鏈讀條碼關(guān)聯(lián)信息。但是從表5可以看出SupperLLEC組裝得到的支架NG50長度分別為28.3 MB、16.1 MB和17.5 MB，遠(yuǎn)遠(yuǎn)大于其他幾個先糾錯再組裝的算法實驗。進一步合并三個數(shù)據(jù)集的結(jié)果，圖4分別畫出Wtdbg2、SuperLLEC和其他所有工具的支架長度分布，發(fā)現(xiàn)SuperLLEC與Wtdbg2和其他工具之間有顯著差異，其大尺寸的支架占比更大。因為SuperLLEC充分利用鏈讀數(shù)據(jù)的關(guān)聯(lián)關(guān)系進行了進一步的支架組裝，而其他糾錯算法無法運用鏈讀數(shù)據(jù)的條碼信息，所以其他糾錯算法無法達(dá)到SuperLLEC算法的組裝長度。

圖4 數(shù)據(jù)集糾錯前后支架序列長度分布比較Fig.4 Length distributions of scaffolds before and after error corrections

表5 NG50長度對比表Table 5 Comparison of NG50 length

另外，通過和NCBI RefSeq真實基因數(shù)據(jù)集比較，如表6所示，在組裝準(zhǔn)確率上，SuperLLEC在三個數(shù)據(jù)集上均明顯優(yōu)于其他糾錯算法，達(dá)到了非常高的準(zhǔn)確率。同時，因為其他糾錯算法是先糾錯后讀序組裝，沒有單獨SNP的預(yù)測，其SNP預(yù)測均來自Wtdbg2的預(yù)測，因此在SNP覆蓋率（即召回率，預(yù)測正確的SNP數(shù)除以SNP總預(yù)測數(shù)）上其他糾錯算法沒有明顯變化。如表7所示，在三個數(shù)據(jù)集上，SuperLLEC的SNP覆蓋率均超過80%，大大提升了Wtdbg2的召回率。因此SuprLLEC可以與Wtdbg2結(jié)合使用來提高組裝的準(zhǔn)確率和SNP預(yù)測精度。

表6 組裝準(zhǔn)確率對比表Table 6 Comparison of assembly accuracy

表7 SNP覆蓋率對比表Table 7 Comparison of SNP coverage

4 結(jié)束語

SuperLLEC算法是借用跨度大的鏈讀數(shù)據(jù)集對PacBio三代測序數(shù)據(jù)進行糾錯和組裝。SuperLLEC算法有別于傳統(tǒng)的糾錯算法：首先，放棄直接糾錯，而采用先組裝后糾錯，并且在糾錯過程中進一步組裝。其次，根據(jù)鏈讀及其條碼信息與PacBio數(shù)據(jù)結(jié)合進行組裝和糾錯，目前主要算法均為第三代數(shù)據(jù)自糾或運用第二代測序短讀數(shù)據(jù)進行糾錯，因此只能借用短讀的雙末端（pair-end）信息，而SuperLLEC首次根據(jù)鏈讀條碼信息進行糾錯和組裝，組裝和糾錯效果更好。再次，算法采用的Fisher精確檢驗比經(jīng)典的卡方檢驗更適合離散分布的檢驗，能夠更好地辨別單核酸多態(tài)性堿基和測序錯誤堿基。最后，算法步驟中多次運用并行計算技巧加快算法的運算速度。目前由于第三代測序技術(shù)還存在測序準(zhǔn)確率偏低的問題，因此該糾錯組裝算法能夠很好地提高基因組組裝的準(zhǔn)確率。