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

崔雅軒，張少強

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

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

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

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

1 數據描述

1.1 10X Genomics鏈讀

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

1.2 PacBio長讀

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

表1 實驗所需PacBio長讀數據Table 1 PacBio long read data in experiment

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

2 方法和步驟

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

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

2.1 組裝長讀數據構建支架序列

對原始的PacBio測序數據，首先用第三代測序的基因組組裝算法將其拼接成支架序列集合。Wtdbg2[9]、FALCON[10]、Canu[11]是三種最新的第三代測序組裝算法。實驗分別嘗試使用這三種算法依次對表1的三組長讀數據集進行序列組裝。對組裝完成的數據進行比較分析后發現，Canu組裝后的人類基因組數據規模（約為3.42 GB）大于Wtdbg2和FALCON的數據規模（約為3 GB）。另外實驗發現Wtdbg2在三者中是運算最快的算法。例如對于Human CHM1的PacBio長讀數據，Canu的CPU用時為22 750 h，FALCON用時為68 789 h，而Wtdbg2用時只有245 h。這是因為Wtdbg2算法采用Linux的Pthreads多線程處理技術來加快算法的執行速度。由此，實驗采用Wtdbg2算法對實驗的PacBio長讀數據組裝成支架序列集合。在運行Wtdbg2時，除了記錄支架序列外，還記錄下堿基多態位點（Wtdbg2預測為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存在于支架的編號集合，據此將每組鏈讀序列集合分配給支架。如圖2示例所示，通過判斷鏈讀和支架是否共享k-mer，將四組鏈讀序列集合依次分配給兩條支架序列。其中Barcode編號為2的鏈讀組同時分配給了兩條支架。

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

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

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

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

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

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

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

2.4 糾錯和SNP預測

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

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

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

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

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

3 實驗與結果分析

3.1 實驗環境

實驗所有數據均在CPU為Intel Xeon E2244G 3.80 GHz 8 Cores，內存為512 GB的CentOS系統服務器上調試運行。

3.2 結果分析

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

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

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

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

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

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

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

表6 組裝準確率對比表Table 6 Comparison of assembly accuracy

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

4 結束語

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