于長(zhǎng)永， 金建宇， 劉 鵬， 趙宇海

(東北大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院， 遼寧 沈陽(yáng) 110169)

DBG(de Bruijn圖)是生物信息學(xué)中最有用的數(shù)據(jù)結(jié)構(gòu)之一，例如，DBG最初用于組裝數(shù)十億個(gè)短基因組序列[1-2]，DBG也被用于其他領(lǐng)域，如基因組序列比對(duì)[3-6]，在人群中發(fā)現(xiàn)基因多樣性[7-8]，預(yù)測(cè)細(xì)菌的毒性[9]，以及檢測(cè)共線性塊[4].

近年來(lái)，由于DBG具有能夠處理基因組重復(fù)區(qū)域[10-11]的特性，而基因組的高度重復(fù)問(wèn)題也是目前大多數(shù)現(xiàn)有方法的瓶頸，因此DBG常用于處理基因組序列.對(duì)于基于FM-index的方法，高重復(fù)的子字符串可能會(huì)導(dǎo)致后綴的定位過(guò)程耗時(shí)較長(zhǎng)；對(duì)于基于gram的方法，由于存在過(guò)濾失敗的情況，高度重復(fù)的子字符串可能會(huì)造成候選集較大.為了避免重復(fù)子字符串的影響，提出了很多種子選擇方法，以O(shè)SS[10]為例，將頻率最小的種子定義為最優(yōu)種子，提出一種求解種子選擇問(wèn)題的動(dòng)態(tài)規(guī)劃算法.此外，在文獻(xiàn)[11]中，用DBG處理參考序列重復(fù)引起的問(wèn)題，但忽略了其他信息，如分支節(jié)點(diǎn)和由分支節(jié)點(diǎn)組成的路徑.DeBGA是一個(gè)功能強(qiáng)大的映射器，但它可能會(huì)丟失由DBG模型中的分支節(jié)點(diǎn)和路徑表示的映射位置.

DBG的一個(gè)重要特性是，參考序列的每個(gè)子字符串都是圖上的路徑，但并非所有DBG上的路徑都是參考序列子字符串.如果一個(gè)路徑是參考序列的子字符串，則稱(chēng)它為真路徑，否則為假路徑.DBG模型不僅可以對(duì)參考序列進(jìn)行索引，還可以在參考序列上定位圖的每條邊和路徑.近年來(lái)有關(guān)DBG的研究主要集中在圖的索引和快速構(gòu)建方面.例如，文獻(xiàn)[12-13]提出了有效構(gòu)造DBG的方法，文獻(xiàn)[14]提出了一種時(shí)間和空間友好的壓縮染色DBG的索引方法，文獻(xiàn)[15]提出了一個(gè)DBG的緊湊表示，它允許無(wú)限數(shù)量的節(jié)點(diǎn)和邊的增加和刪除.然而，關(guān)于在參考序列上定位DBG的節(jié)點(diǎn)、邊和路徑的方法卻很少.文獻(xiàn)[16]提出了一種新的圖數(shù)據(jù)組裝結(jié)構(gòu)，稱(chēng)為連接的de Bruijn圖(LdBG)，但該模型不能用于定位參考序列上的路徑，顯然也無(wú)法存儲(chǔ)所有邊的位置列表.

本文提出了一種DBG模型，稱(chēng)為MiniDBG，可以用于在read-mapping中索引參考序列.該模型存儲(chǔ)了最小邊集的位置列表.有了位置列表，MiniDBG可以有效地定位任何邊、節(jié)點(diǎn)和路徑.實(shí)驗(yàn)結(jié)果表明，MiniDBG能有效地在參考序列上定位邊和路徑，且存儲(chǔ)成本較低.

1 基礎(chǔ)知識(shí)

X為字母表Σ={A,C,G,T}上的基因組參考序列，A,C,G,T是生物學(xué)DNA中的四種堿基.X長(zhǎng)度用|X|來(lái)表示.X[i]表示X上的第i個(gè)字符.索引從1開(kāi)始，用X[i,j]表示X上從i到j(luò)的子字符串，X的前綴和后綴由X[*,i]和X[i,*]來(lái)表示.給定一個(gè)參考序列X和子字符串k-mer，k-mer是基因組參考序列X中長(zhǎng)度為k的字符串，用DBG(X,k)來(lái)表示X的 de Bruijn圖：

V={vi|vi是k-mer并且vi是X上的字符串}；

E={eij|vi[k-1,*]=vj[*,k-1],|pij|≠0}；

P={pij|pij=list(eij)}.

其中V,E,P分別為圖的頂點(diǎn)集合、邊集合和位置列表集合.集合V由序列X上出現(xiàn)的一系列長(zhǎng)度為k的字符串組成.eij是一條從vi到vj并且滿足vi[k-1,*]=vj[*,k-1]的邊，并且|pij|≠0.pij是eij的位置列表(將eij看作X上長(zhǎng)度為k+1的子字符串)，能夠按照從前到后的順序保存eij在X上出現(xiàn)的位置.|pij|表示pij的長(zhǎng)度.

給定參考序列X和整數(shù)k，找到邊的位置列表的集合，用Pmin表示，并需要滿足以下條件：

①DBG′(X,k)=(V,E,Pmin)可以通過(guò)Pmin的位置列表的并集定位每一條邊；

②對(duì)于任意的P′?Pmin,DBG″(X,k)=(V,E,P′)不能滿足條件1.

DBG′(X,k)=(V,E,Pmin)為參考序列X上字符串長(zhǎng)度為k的MiniDBG.

2 路徑定位算法

2.1 MiniDBG模型

MiniDBG是傳統(tǒng)DBG的一種壓縮存儲(chǔ)形態(tài)，通過(guò)對(duì)分支、非分支節(jié)點(diǎn)進(jìn)行區(qū)分，找到圖中的超級(jí)邊，以及超級(jí)邊所對(duì)應(yīng)的最小邊，最終生成最小邊的位置列表.對(duì)最小邊及最小邊位置列表進(jìn)行存儲(chǔ)，從而實(shí)現(xiàn)傳統(tǒng)DBG的壓縮存儲(chǔ).

定義1(非分支節(jié)點(diǎn)、單路徑)：

如果DBG(X,k)上的一個(gè)頂點(diǎn)vi滿足：入度與出度相等且為1，則vi被稱(chēng)為非分枝節(jié)點(diǎn)，否則為分支節(jié)點(diǎn).對(duì)于一個(gè)DBG(X,k)上的一條路徑h=vi1vi2…vik，如果其每個(gè)節(jié)點(diǎn)都是非分支節(jié)點(diǎn)，則稱(chēng)之為單路徑.如果vi的入度不為1，則稱(chēng)vi為入分支節(jié)點(diǎn)；如果vi的出度不為1，則稱(chēng)vi為出分支節(jié)點(diǎn).對(duì)于一條邊eij，如果vi是出分支節(jié)點(diǎn)，則eij為一條出邊.

假設(shè)邊a和b是相鄰的，a為節(jié)點(diǎn)的入邊,b為節(jié)點(diǎn)的出邊，相鄰邊位置的關(guān)系分為4種：

①a是非分支入邊，b是非分支出邊；

②a是非分支入邊，b是分支出邊；

③a是分支入邊，b是非分支出邊；

④a是分支入邊，b是分支出邊.

用a=b,a>b,a

定義2(超級(jí)邊)：

在移位操作下，超級(jí)邊內(nèi)各邊的位置列表相同，因此將超級(jí)邊的位置列表定義為其第一條邊的位置列表：

在某些情況下，邊的超級(jí)邊仍然是它本身.因?yàn)槌?jí)邊上的邊的位置列表彼此是等價(jià)的，超級(jí)邊的位置列表可以表示超級(jí)邊上所有的邊.為了消除超級(jí)邊中等價(jià)位置列表的影響，DBG可以看作是一個(gè)新的圖，稱(chēng)為super-DBG，它由分支節(jié)點(diǎn)和超級(jí)邊組成.

表1 四種鄰接邊的關(guān)系

算法1給出了超級(jí)邊的生成過(guò)程.T表示用于保存結(jié)果的節(jié)點(diǎn)列表，第2～6行為向右擴(kuò)展邊緣的過(guò)程，直至到達(dá)分支節(jié)點(diǎn).同樣，第7～11行向左擴(kuò)展邊緣，直至到達(dá)分支節(jié)點(diǎn).T在擴(kuò)容過(guò)程中按順序依次保存節(jié)點(diǎn).

算法1 生成給定邊的超級(jí)邊Gen_super_edge(ei0j0)：

輸入: 圖DBG(X,k)=(V,E), 一條邊ei0j0;

1:T=?;k=0;l=0;

2: while 1

3: insertvjkintoTat the tail termination;

4: ifvjkis branching

5: break;

6:vjk+1=vjk.successor();k++;

7:while 1

8: insertvilintoTat the head termination;

9: ifvilis branching

10: break;

11:vil+1=vik.precursor();l++;

12: ReturnT.

定義3(最小超級(jí)邊)：

算法2描述了基于MiniDBG模型計(jì)算一條邊的位置列表的過(guò)程.第3～4行對(duì)應(yīng)于輸入邊的超級(jí)邊是最小超級(jí)邊的情況(Δ為偏移量)，超級(jí)邊中第一條邊的位置列表保存在MiniDBG模型中.第5～18行對(duì)應(yīng)輸入邊的超級(jí)邊不是最小超級(jí)邊的情況；這里分兩種情況：①超級(jí)邊從入分支節(jié)點(diǎn)開(kāi)始；②超級(jí)邊在出分支節(jié)點(diǎn)結(jié)束.其中El-表示從頂點(diǎn)l出發(fā)的邊的集合，E-k表示在頂點(diǎn)k結(jié)束的邊的集合.

顯然，超級(jí)邊可以從入分支節(jié)點(diǎn)開(kāi)始，同時(shí)又在出分支節(jié)點(diǎn)結(jié)束.在這種情況下，位置列表可以通過(guò)情況①或②計(jì)算.

算法2 計(jì)算每條邊的位置列表Gen_list(eij)：

輸入:圖DBG(X,k)=(V,E,Pmin), 一條邊eij;

輸出:list(eij)，其中l(wèi)ist(eij)=Gen_list(eij);

1:list(eij)=?;

5:else ifvlis out-branching

6: foreluinEl-

8: ifvyis in-branching

9:

10: else

11: list(eij)=list(eij)∪(Gen_list(elu)+Δ);

12:else ifvkis in-branching

13: foreukinE-k

15: ifvyis out-branching

17: else

18: list(eij)=list(eij)∪(Gen_list(euk)+Δ);

19: Return list(eij).

2.2 MiniDBG路徑定位算法

定理1DBG中的閉環(huán)路徑示例如圖1所示，對(duì)于一個(gè)閉環(huán)路徑，在循環(huán)路徑中至少存在一個(gè)入分支節(jié)點(diǎn)和一個(gè)出分支節(jié)點(diǎn).

圖1 DBG中的閉環(huán)路徑

定理2將整條路徑h=ei0i1ei1i2…ein-1in進(jìn)行分割，得到路徑片段hseg=eikik+1eik+1ik+2…eik+l-1ik+l，則有以下結(jié)論：

①?q,k≤q≤k+l-1滿足list(hseg)=list(eiqiq+1)-(q-k).

②eiqiq+1是最小超級(jí)邊.

③hseg中除eiqiq+1之外沒(méi)有任何邊是最小超級(jí)邊.

定理3 圖DBG(X,k)上的路徑h=ei0i1ei1i2…ein-1in一定滿足以下結(jié)論：

①h可以用一個(gè)超級(jí)邊序列等價(jià)地表示

②路徑上的邊有相同的位置列表list(h)=list(h′);

③

根據(jù)定理3可以推得基于MiniDBG模型的路徑位置列表計(jì)算方法，如算法3所示.第1行，L=Ω表示結(jié)果位置列表，并初始化為一般集合.Rpre表示初始的邊與當(dāng)前的超級(jí)邊之間的關(guān)系.類(lèi)似地，Rcur表示當(dāng)前超級(jí)邊和下一條邊之間的關(guān)系.l表示當(dāng)前超級(jí)邊的長(zhǎng)度.Get_relationship(*,*)表示根據(jù)連接兩條邊的節(jié)點(diǎn)的類(lèi)型獲取兩條相鄰邊之間的關(guān)系.第4～5行，算法跳過(guò)“=”關(guān)系，如果某一步驟后的交集為空，則算法返回的路徑為假路徑.

算法3 計(jì)算給定路徑的位置列表：

輸入:DBG(X,k)=(V,E,Pmin), 一路徑

ei0i1ei1i2…ein-1in;

輸出: list(ei0i1ei1i2…ein-1in);

1:L=Ω;Rpre=′′;Rcur=′?′;l=1;

2: for(k=0;k

3:Rpre=Rcur;k=k+l-1;l=1;

4: while(Rcur=Get_relationship(eik+lik+l+1,eik+l+1ik+l+2)=′=′)

5:l++;

6: ifRpre=′?′orRpre=′>′

7: ifRcur=′?′orRcur=′<′

8:L=L∩(list(eikik+1)-k);

9: ifL=?

10: returnL;

11:k=k+l-1;

12: ifRcur=′>′orRcur=′?′

13:L=L∩(list(ein-1in)-n+1);

14: ReturnL.

為了定位一個(gè)點(diǎn)vi，MiniDBG方法必須定位從vi開(kāi)始或在vi結(jié)束的邊.然后結(jié)合所有邊的位置生成節(jié)點(diǎn)vi的位置列表.定位節(jié)點(diǎn)比定位邊花費(fèi)更多的時(shí)間.然而，k為n+1的圖上的節(jié)點(diǎn)列表相當(dāng)于k取n的圖上的邊列表.因此，MiniDBG可以通過(guò)設(shè)置較小的k來(lái)減少定位所花費(fèi)的時(shí)間.

對(duì)于邊的定位，如果是一條最小邊，MiniDBG將返回邊的位置列表；對(duì)于非最小邊，MiniDBG會(huì)生成一組最小邊，并結(jié)合最小邊位置列表得到原始邊的位置列表.這個(gè)過(guò)程遞歸地調(diào)用定位邊的方法，直到所有邊都最小為止.對(duì)于路徑，MiniDBG首先檢查所有的節(jié)點(diǎn)和邊是否在圖上：如果不是，則返回空列表；否則，MiniDBG將使用定理3中的公式計(jì)算路徑的位置列表.在此過(guò)程中，一旦獲取的位置列表的交集為空，算法返回并且路徑的位置列表為空.所有這些操作都是為了避免在定位假路徑上浪費(fèi)時(shí)間.

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

所有實(shí)驗(yàn)在一臺(tái)搭載了4個(gè)Intel(R)Xeon(R)E5-2640 2.40 GHz CPU的服務(wù)器上完成，每個(gè)CPU有10個(gè)核心，服務(wù)器內(nèi)存空間為160 GB ，并搭載了4臺(tái) NVDIA Tesla K40 顯卡.操作系統(tǒng)為 Ubuntu 16.04.SHDex算法在eclipse中實(shí)現(xiàn)，并使用 C 和C++ 編程.

本文使用了以下數(shù)據(jù)集：從NCBI RefSeq數(shù)據(jù)庫(kù)下載了所有的62E大腸桿菌菌株數(shù)據(jù)[17](http://www. ncbi. nlm. nih. gov/refseq/)，選取其中之一稱(chēng)為 D1數(shù)據(jù)集；62E中全部數(shù)據(jù)集稱(chēng)為D2；文獻(xiàn)[12]中引用的人類(lèi)基因組數(shù)據(jù)在本文中被用作數(shù)據(jù)樣本D3.

k-mer長(zhǎng)度為k的頂點(diǎn)對(duì)應(yīng)的圖和k-mer長(zhǎng)度為k-1的邊對(duì)應(yīng)的圖是一樣的；因此，MiniDBG可以設(shè)置較小的k-mer長(zhǎng)度.對(duì)于FM-index，當(dāng)后綴數(shù)組間隔為空時(shí)，路徑為假路徑，定位過(guò)程結(jié)束.

測(cè)試MiniDBG在增加查詢長(zhǎng)度的同時(shí)定位字符串的性能.當(dāng)字符串長(zhǎng)度L增加時(shí)，查找子字符串的時(shí)間開(kāi)銷(xiāo)Cme增加.結(jié)果如圖2所示.MiniDBG的內(nèi)存開(kāi)銷(xiāo)非常小.隨著字符串長(zhǎng)度L的增加，MiniDBG耗費(fèi)了可接受的時(shí)間來(lái)定位這些字符串.在圖2c中，隨著查詢子字符串長(zhǎng)度L的增加，時(shí)間成本增加得非常緩慢.

圖3～圖5為在數(shù)據(jù)集D1，D2和D3上定位子字符串的時(shí)間開(kāi)銷(xiāo)，圖例中F和M分別代表FM-index和MiniDBG.查找子字符串是由參考序列上的滑動(dòng)窗口生成的.查詢總數(shù)為5 521 840，圖中顯示了定位查找子字符串的總時(shí)間開(kāi)銷(xiāo).對(duì)于FM-index，查找的子字符中長(zhǎng)度分別從9/10/15到17/20/23.對(duì)于MiniDBG，查找的子字符串長(zhǎng)度從k+1到k+9，參數(shù)k分別為9/10/19 到 14/15/21.

從實(shí)驗(yàn)結(jié)果中得到如下結(jié)論：

1) 對(duì)于非頻繁查詢(特別是在參考序列上只出現(xiàn)一次的查詢)，F(xiàn)M-index比MiniDBG效率更高，如圖3所示.對(duì)于在參考序列上出現(xiàn)超過(guò)10次的頻繁查詢，MiniDBG的效率更高，如圖4所示，這是因?yàn)镕M-index和MiniDBG的工作方式不同.FM-index首先計(jì)算查找字符串的后綴數(shù)組間隔，然后針對(duì)間隔內(nèi)的每個(gè)索引計(jì)算后綴數(shù)組的值；因此，間隔越大，花費(fèi)的時(shí)間就越多.對(duì)于MiniDBG，它可以一次性獲得所有位置，而不像FM-index需要依次計(jì)算每個(gè)后綴數(shù)組.

2) MiniDBG在字符串定位時(shí)比FM-index更靈活.在原子串上增加一個(gè)字符的定位問(wèn)題中，F(xiàn)M-index 需要重新計(jì)算后綴數(shù)組的所有值.例如，給定“CGT”的位置，只添加一個(gè)基數(shù)的“ACGT”或“CGTA”的位置必須重新計(jì)算得出；而對(duì)于MiniDBG，只需要添加一條邊，可以通過(guò)將新邊的列表與原始列表合并來(lái)計(jì)算.

圖2 子字符串長(zhǎng)度增加時(shí)查找子字符串的時(shí)間成本

3) MiniDBG查找字符串的時(shí)間開(kāi)銷(xiāo)比FM-index更低，即使它比FM-index占用更少的內(nèi)存.如圖4所示，在數(shù)據(jù)集D2上，F(xiàn)M-index在需要566 MB內(nèi)存的情況下，耗費(fèi)104，94和 94 s 來(lái)分別定位長(zhǎng)度為14，15和16的字符串；而MiniDBG能夠分別使用519，221和221 MB內(nèi)存進(jìn)行同樣的查詢操作，分別只需要41，76和93 s.如圖5所示，在數(shù)據(jù)集D3上，F(xiàn)M-index使用5 GB的內(nèi)存，在2 150，1 715，1 396和1 132 s內(nèi)定位長(zhǎng)度為20，21，22，23的字符串；而MiniDBG分別使用3.6，3，3和3 GB的內(nèi)存進(jìn)行相同的查詢操作，僅需515，955，900和902 s.

圖3 在數(shù)據(jù)集D1上定位子字符串的時(shí)間開(kāi)銷(xiāo)

圖4 在數(shù)據(jù)集D2上定位子字符串的時(shí)間開(kāi)銷(xiāo)

圖5 在數(shù)據(jù)集D3上定位子字符串的時(shí)間開(kāi)銷(xiāo)

4 結(jié) 語(yǔ)

本文提出了MiniDBG模型，用于存儲(chǔ)最小邊位置列表，以便在參考序列上定位圖的邊和路徑，并給出了基于MiniDBG模型的路徑定位算法.討論了k-mer長(zhǎng)度增加時(shí)位置列表所占內(nèi)存的大小.通過(guò)實(shí)驗(yàn)驗(yàn)證了模型和算法的有效性.