基于de Bruijn圖和序列比對的長序列混合糾錯(cuò)算法

劉 剛

（廣西大學(xué)計(jì)算機(jī)與電子信息學(xué)院，南寧 530004）

0 引言

基因糾錯(cuò)對于基因工程是很重要的。序列糾錯(cuò)算法旨在識(shí)別和消除（或修復(fù)）測序序列中包含的錯(cuò)誤，從而有利于重新測序或從頭測序分析。序列糾錯(cuò)算法應(yīng)能有效地處理越來越多、越來越長的測序數(shù)據(jù)的糾錯(cuò)。

為了糾正第三代測序平臺(tái)產(chǎn)生的錯(cuò)誤率高的長序列，人們開發(fā)了混合糾錯(cuò)算法。所謂混合糾錯(cuò)算法是指利用第二代測序平臺(tái)產(chǎn)生的高準(zhǔn)確率的短序列來對第三代測序長序列進(jìn)行糾錯(cuò)的算法。長序列混合糾錯(cuò)算法分為4類：基于短序列比對的混合糾錯(cuò)算法，基于重疊序列或重疊群比對的混合糾錯(cuò)算法，基于de Bruijn 圖的混合糾錯(cuò)算法以及基于隱馬爾可夫模型的混合糾錯(cuò)算法?；诙绦蛄斜葘突谥丿B序列比對的混合糾錯(cuò)算法由于長序列的較高錯(cuò)誤率導(dǎo)致比對會(huì)出現(xiàn)偏差，因此需要改進(jìn)比對方式來提高糾錯(cuò)效果，且它們將長序列拆分成很多的長度很短的序列片段，這些序列片段會(huì)丟失長序列所攜帶的信息，這使得糾錯(cuò)效果不夠理想?；陔[馬爾可夫模型的混合糾錯(cuò)算法也需要將長序列與短序列比對，也需要改進(jìn)比對方式來提高糾錯(cuò)效果。已有的基于de Bruijn 圖的混合糾錯(cuò)算法采用固定值的k-mers（長度為的序列片段）構(gòu)建，其有時(shí)很難找到長序列的有效錨點(diǎn)，從而難以擴(kuò)展序列路徑，這使得糾錯(cuò)質(zhì)量不高。

本文研究基于值可變的de Bruijn 圖和序列比對的長序列混合糾錯(cuò)算法，以糾正第三代測序長序列中與第二代測序短序列比對的區(qū)域和未被覆蓋的區(qū)域。

1 算法

1.1 算法思想

本文的算法思想：將與長序列同物種的短序列采用短序列糾錯(cuò)算法其進(jìn)行糾錯(cuò)，從糾錯(cuò)后的短序列中提取出k-mers（長度為的序列片段），將待糾錯(cuò)的第三代測序長序列與同物種的正確率高的第二代測序短序列進(jìn)行比對，以生成種子；使用值可變的de Bruijn圖擴(kuò)展連接形成種子序列，將連接兩個(gè)相鄰種子的序列路徑覆蓋位于兩個(gè)種子之間未與長序列比對的區(qū)域；遍歷de Bruijn 圖擴(kuò)展種子序列的末端，使得種子序列能夠擴(kuò)展到待糾錯(cuò)的長序列的末端；采用短序列比對糾正長序列與短序列對準(zhǔn)的區(qū)域，并使用種子序列路徑來糾正長序列未與短序列對準(zhǔn)的區(qū)域，從而得到經(jīng)糾錯(cuò)之后的長序列。

1.2 算法描述

算法是通過將待糾錯(cuò)的長序列與同物種的短序列比對來生成種子的，種子由5 元組（id，pos，len，score，seq）組成，其中id 表示種子與長序列關(guān)聯(lián)的標(biāo)識(shí)符，pos 表示相同物種的短序列與長序列比對的開始位置，len 表示短序列與長序列比對區(qū)域的長度，score 表示比對得分，seq 表示在某個(gè)位置上短序列與長序列比對的共有序列。

本文給出的基于值可變de Bruijn圖和序列比對的長序列混合糾錯(cuò)算法（hybrid de Briujin graph errors correction algorithm，簡記為HdGEC算法）形式描述如下。

HdGEC

輸入：相同物種的長序列l(wèi)ong-read[0..n-1]和短序列short-reads[0..m-1]，參考基因組R[0..r-1]

輸出：糾正后的長序列l(wèi)ong-read[0...n-1]

Begin

1:采用QuorUM 算法糾正短序列short-reads[0..m-1］，以減少短序列中包含的錯(cuò)誤；

2: 使用KMC3 工具從糾錯(cuò)后的短序列short-reads［0..m-1］中提取出k-mers；

3:過濾掉弱值的k-mers，將強(qiáng)值的k-mers 使用PgSA索引構(gòu)造值可變的de Bruijn圖；

4: 使用BLASR 工具將經(jīng)過糾錯(cuò)的短序列shortreads[0..m-1］與同物種的長序列l(wèi)ong-read[0...n-1]比對，以獲得種子seeds（id,pos,len,score,seq）；

5: 利用序列比對將種子序列seeds（id, pos, len,score,seq）與參考基因組R[0..r-1]進(jìn)行比對，以糾正中被短序列short-reads[0..m-1］覆蓋到的區(qū)域；

6:遍歷值可變的de Bruijn 圖以找到任意相鄰兩個(gè)種子之間的序列路徑，連接這相鄰的兩個(gè)種子得到種子序列路徑，并用種子序列路徑來糾正長序列l(wèi)ong-read[0...n-1]當(dāng)中未被短序列short-reads[0..m-1]覆蓋到的區(qū)域；

7:遍歷值可變的de Bruijn 圖，以擴(kuò)展得到種子序列路徑的末端，使得種子序列路徑覆蓋到整條長序列l(wèi)ong-read[0...n-1]；

8:輸出糾正之后的長序列l(wèi)ong-read[0...n-1]；

End.

本文算法HdGEC 使用種子作為錨點(diǎn)進(jìn)行擴(kuò)展，使用貪心策略遍歷de Bruijn 圖，由于可以選擇1～階的de Bruijn 圖來連接種子序列，所以可以減少de Bruijn 圖的遍歷次數(shù)，從而縮短了序列糾錯(cuò)需要的時(shí)間。此外，本文算法HdGEC 結(jié)合k-mers 比対和值可變de Bruijn 圖來對長序列進(jìn)行糾錯(cuò)，其中通過k-mers 序列比対來糾正長序列中被短序列覆蓋到的區(qū)域，通過de Bruijn 圖糾正長序列中沒有被短序列覆蓋到的區(qū)域，這使得算法的覆蓋率更高，能夠糾正序列中更多的錯(cuò)誤堿基，獲得了更高的吞吐量（輸出的糾錯(cuò)序列的質(zhì)量）。

2 實(shí)驗(yàn)

2.1 實(shí)驗(yàn)環(huán)境與數(shù)據(jù)

實(shí)驗(yàn)使用的計(jì)算機(jī)是8 核Intel（R）Core i7-6700 k CPU@4.00 GHz處理器、內(nèi)存容量32 GB。采用C++語言編程實(shí)現(xiàn)算法。實(shí)驗(yàn)數(shù)據(jù)集采用太平洋生物科學(xué)平臺(tái)的PacBio模擬數(shù)據(jù)集和牛津納米孔平臺(tái)的ONT真實(shí)數(shù)據(jù)集。模擬數(shù)據(jù)集與真實(shí)數(shù)據(jù)集都來自于物種A.baylyi、E.coli、S.cere?visiae和C.elegans 中的數(shù)據(jù)，這4 種物種數(shù)據(jù)集和參考基因組來源于基因測序中心https://www.genoscope.cns.fr/externe/NaS/datasets。模擬的長序列數(shù)據(jù)集采用SimLord模擬器生成。實(shí)驗(yàn)使用的參考基因組、模擬數(shù)據(jù)集與真實(shí)數(shù)據(jù)集的信息如表1所示。

表1 實(shí)驗(yàn)數(shù)據(jù)集

2.2 模擬數(shù)據(jù)集上的長序列糾錯(cuò)實(shí)驗(yàn)

對于模擬數(shù)據(jù)集上的實(shí)驗(yàn)，采用糾正后的長序列的錯(cuò)誤率（%）、吞吐量、split reads比率、運(yùn)行時(shí)間這4個(gè)指標(biāo)測試長序列混合糾錯(cuò)算法的性能，其中、和的值越小越好，值越大越好。使用LRCStats 軟件測量得到算法CoLoRMap、HALC、Jabba、LoRDEC、NaS和HdGEC 在物種A.baylyi和E.coli 模擬數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果，分別如表2和表3所示。

從表2 與表3 的實(shí)驗(yàn)結(jié)果可以看到：在糾正后的長序列的錯(cuò)誤率上，6 種混合糾錯(cuò)算法的值均低于0.3%，本文算法HdGEC 對物種A.baylyi 糾正后的長序列的錯(cuò)誤率最低，對物種E.coli 糾正后的長序列的錯(cuò)誤率是第二低的，它比Jabbaa 算法糾正后的長序列的錯(cuò)誤率0.0462%僅高了0.0134%。在吞吐量上，除了NaS算法，其他5種混合糾錯(cuò)算法的指標(biāo)值都較高，并且本文算法HdGEC 的吞吐量高于其他5 種混合糾錯(cuò)算法。在split reads 比率上，算法CoLoRMap、HALC和LoRDEC 的srr值較高，值高意味著長序列的大部分區(qū)域沒有被糾正，本文算法HdGEC 的值較低，表明參考基因組與短序列比對的區(qū)域有很大的機(jī)會(huì)被發(fā)現(xiàn)，長序列被短序列覆蓋到的更多區(qū)域?qū)?huì)被糾正；在運(yùn)行時(shí)間方面，糾錯(cuò)算法Jabba 所需時(shí)間最少，而糾錯(cuò)算法NaS所需時(shí)間最多。

表2 長序列糾錯(cuò)算法在物種A.baylyi上模擬數(shù)據(jù)集的實(shí)驗(yàn)結(jié)果

表3 糾錯(cuò)算法在物種E.coli上模擬數(shù)據(jù)集的實(shí)驗(yàn)結(jié)果

這表明，對于物種E.coli和A.baylyi 的模擬數(shù)據(jù)集的長序列糾錯(cuò)，本文算法HdGEC 的吞吐量高于其他算法，并且錯(cuò)誤率較低，產(chǎn)生的split reads 的比例也較小，能糾錯(cuò)長序列的大部分區(qū)域。

2.3 真實(shí)數(shù)據(jù)集上的長序列糾錯(cuò)實(shí)驗(yàn)

對于在牛津納米孔真實(shí)數(shù)據(jù)集上的實(shí)驗(yàn)，采用糾正的長序列數(shù)量、比率（%）、糾正的長序列平均長度、糾正的堿基數(shù)量、平均同一性、覆蓋率、運(yùn)行時(shí)間這7個(gè)指標(biāo)測試長序列混合糾錯(cuò)算法的性能，其中、、、和的 值 越 大越好，和值越小越好。長序列糾錯(cuò)算法HdGEC、 CoLoRMap、 HALC、 Jabba、LoRDEC和NaS在小型物種A.baylyi和E.coli、中型物種S.cerevisiae、大型物種C.elegans 真實(shí)數(shù)據(jù)集上運(yùn)行的實(shí)驗(yàn)結(jié)果分別如表4～表7所示。

表7 算法在大型物種C.elegans真實(shí)數(shù)據(jù)集運(yùn)行的實(shí)驗(yàn)結(jié)果

從表4和表5可以看到，對于小型物種真實(shí)數(shù)據(jù)集的長序列糾錯(cuò)，在6個(gè)算法中，本文算法HdGEC 糾正的長序列的平均長度是最長；本文算法HdGEC的split reads比率第二低，僅比算法NaS高了一點(diǎn)點(diǎn)；本文算法HdGEC、CoLoRMap、HALC、LoRDEC和NaS 的覆蓋率均達(dá)到100%；算法HALC 能糾正的長序列數(shù)量最多，本文算法HdGEC 能糾正的長序列數(shù)量第二多；在糾正的堿基數(shù)量方面，本文算法HdGEC和NaS 是較多的；在運(yùn)行時(shí)間方面，算法Jabba 所需時(shí)間最少，本文算法所需時(shí)間是第三少。

表4 算法在小型物種A.baylyi真實(shí)數(shù)據(jù)集的實(shí)驗(yàn)結(jié)果

表5 算法在小型物種E.coli真實(shí)數(shù)據(jù)集的實(shí)驗(yàn)結(jié)果

從表6可以看到，對于中型物種真實(shí)數(shù)據(jù)集的長序列糾錯(cuò)，在6個(gè)算法中，本文算法HdGEC 能糾正的長序列平均長度最長、能糾正的堿基數(shù)量最多、split reads 比率最低、覆蓋率最高；在能糾正的長序列數(shù)量方面，算法HALC最多，本文算法HdGEC 第二多；在平均同一性方面，算法Jabba 最高，本文算法HdGEC 第三高；在運(yùn)行時(shí)間方面，Jabba 算法所需時(shí)間最少，本文算法HdGEC 所需時(shí)間是第四少；而算法NaS運(yùn)行時(shí)間超過了16天尚未計(jì)算得到結(jié)果。

表6 算法在物種S.cerevisiae真實(shí)數(shù)據(jù)集的實(shí)驗(yàn)結(jié)果

從表7可以看到，對于大型物種真實(shí)數(shù)據(jù)集的長序列糾錯(cuò)，在6個(gè)算法中，本文算法HdGEC 能糾正的長序列數(shù)量最多、能糾正的長序列平均長度最長、能糾正的堿基數(shù)量最多、split reads 比率最低、覆蓋率最高；在平均同一性方面，算法Jabba 最高，本文算法HdGEC 第二高；在運(yùn)行時(shí)間方面，Jabba 算法所需時(shí)間最少，本文算法HdGEC 所需時(shí)間是第三少；而算法HALC和NaS 運(yùn)行時(shí)間超過了16 天尚未計(jì)算得到結(jié)果。

上述結(jié)果表明：對于中型和大型物種真實(shí)數(shù)據(jù)集的長序列糾錯(cuò)，本文算法HdGEC 的整體糾錯(cuò)質(zhì)量最好；對于小型物種真實(shí)數(shù)據(jù)集的長序列糾錯(cuò)，本文算法HdGEC 的整體糾錯(cuò)質(zhì)量也較好。

3 結(jié)語

為了有效糾正第三代測序平臺(tái)產(chǎn)生的錯(cuò)誤率高的長序列，本文給出一個(gè)基于值可變de Bruijn圖和序列比對的混合糾錯(cuò)算法。在模擬數(shù)據(jù)集和真實(shí)數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果表明，與已有的長序列混合糾錯(cuò)算法相比，本文算法獲得了整體較好的糾錯(cuò)質(zhì)量。本文算法為了獲得較高質(zhì)量的長序列糾錯(cuò)結(jié)果，是利用de Bruijn 圖來尋找種子序列路徑，這需要花費(fèi)較時(shí)間長。下一步工作是研究將本文的算法并行化，以在維持獲得較高質(zhì)量的長序列糾錯(cuò)結(jié)果的同時(shí)，顯著減少糾錯(cuò)過程所需的時(shí)間。