崔雅軒，張少?gòu)?qiáng)

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

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

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

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

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

1.1 10X Genomics鏈讀

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

1.2 PacBio長(zhǎng)讀

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

表1 實(shí)驗(yàn)所需PacBio長(zhǎng)讀數(shù)據(jù)Table 1 PacBio long read data in experiment

表2 實(shí)驗(yàn)所需PacBio長(zhǎng)讀數(shù)據(jù)下載地址Table 2 Download site of PacBio long read data

2 方法和步驟

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

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

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

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

2.2 鏈讀序列組分配給支架

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

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

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

2.3 鏈讀序列比對(duì)到支架并進(jìn)一步組裝

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

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

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

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

2.4 糾錯(cuò)和SNP預(yù)測(cè)

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

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

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

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

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

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

3.1 實(shí)驗(yàn)環(huán)境

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

3.2 結(jié)果分析

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

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

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

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

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

表5 NG50長(zhǎng)度對(duì)比表Table 5 Comparison of NG50 length

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

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

表7 SNP覆蓋率對(duì)比表Table 7 Comparison of SNP coverage

4 結(jié)束語(yǔ)

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