基于Contig的單面基因組框架填充2-近似算法

柳 楠,卞忠勇,李 洋,朱永琦

(山東建筑大學 計算機科學與技術學院,山東 濟南 250101)

0 引 言

近年來,為了高效獲得完整的生物基因組序列,基因測序和組裝技術的發展受到廣泛關注[1-2]。很多時候,使用生物測序手段無法直接獲得完整的基因組序列,通過使用計算機技術將多次測序獲取的基因組框架進行組裝,將缺失基因填充到基因組框架,能夠提高基因組序列的完整性[3-4]。

前面兩問題中基因組框架的基本單元都是單個基因,對缺失基因插入位置不做限制。但實際應用中,框架多由一系列片段重疊群(contig)組成,缺失基因只允許在contig的兩端插入[11]。對于基于contig的單面框架填充的研究[12],Jiang Haitao等人通過將哈密頓路徑問題多項式變換到本問題,證明它是NP完全問題,然后使用貪婪策略和最大匹配提出2-近似算法[13],Tan Guanlan等人指出這個算法不具一般性,并通過構造輔助圖和最大匹配設計驗證了2.57-近似算法[14]。Bulteau L等人基于最大鄰接數、最小斷點距離提出了k-Mer參數的FPT算法[15],Ma Jingjing等人提出以duo-preservations作為度量,設計實現了近似比為2的一個基本算法,并通過貪婪策略將近似比提高到1.71[16]。

該文研究的問題是基于contig的單面含重復基因組框架填充問題,研究發現該問題的2-近似算法不能避免出現冗余公共鄰接,使近似性能比無法達到2,文獻[14]提出了一種解決方案,但近似性能比只能達到2.57。該文通過構造以缺失基因、基因位點、斷點為節點的輔助圖,配合最大匹配、貪婪策略,設計了一個近似算法,證明了算法的近似性能比可達到2,并實現了算法的驗證程序。

1 相關定義

給定集合Σ,Σ包含所有用于標識基因的字符。設A是由Σ中元素構成的字符串,A中每個元素可多次出現,稱A為序列。用M(A)表示A中的字符集合,由于M(A)可能是多重集,為方便計算,用NM(A)記錄每個字符的重數。例如Σ={a,b,1,2},A=ab12a2,M(A)={a,a,b,1,2,2},NM(Α)={a:2,b:1,1:1,2:2}。若x和y在A中相鄰,則稱xy為A中的一基因對,用P(A)表示A中所有基因對的集合,P(A)也可能是多重集,用NP(A)記錄基因對重數。前面例子中P(A)={ab,b1,12,2a,a2},NP(A)={ab:1,b1:1,12:1,2a:2}。易知基因對2a與a2無區別,故后面討論中xy與yx視為相同基因對。

給定兩條序列Α=a1a2…am和B=b1b2…bn,若P(A)中的一基因對aiai+1,在P(B)可以找到bkbk+1與之對應(aiai+1=bkbk+1或aiai+1=bk+1bk),稱aiai+1和bkbk+1是公共鄰接,用a(A,B)記錄A和B全部的公共鄰接,易知DA(A,B)=P(A)-a(A,B)和DB(A,B)=P(B)-a(A,B)分別表示A和B中沒有匹配到公共鄰接的基因對,這些基因對稱為A或B中的斷點。示例如圖1所示。

圖1 公共鄰接和斷點的示例

在基于contig的情況下,兩條序列中一條是完整的基因序列G,另一條是待填充的基因框架S=?[C1]?[C2]?…?[Cn]?,S由多個contig序列Ci組成,contig兩邊的?為可插入的基因位點,稱為slot。

下面正式給出基于contig的單面基因組框架填充問題(contig-based one-sided scaffold filling, Contig-One-Sided-SF-max)的定義。

定義1:Contig-One-Sided-SF-max。

輸入:一條含重復基因的完整基因序列G=g1g2…gm和一條缺失部分基因的框架S=?[C1]?[C2]?…?[Cn]?,其中gi為單個基因,Ci是contig,且P(S)?P(G),m>n。設X=M(G)-M(S)為S中的缺失基因集合,X≠?。

問題:將X中的基因插入S中的slot得到S',使得|a(G,S')|最大。

輸出:填充后的完整序列S'。

2 Contig-One-Sided-SF-max近似算法

給定實例(G,S),研究目標是設計算法將缺失基因插入S得到S',實現公共鄰接數|a(G,S')|最大化。初始公共鄰接數a(G,S)由已知可直接計算,故算法關注的重點是|a(G,S')-a(G,S)|,即新增公共鄰接數量。另外注意到a(G,S')-a(G,S)?DG(G,S),因此在算法中判斷基因插入位置需要以G中的斷點集合DG(G,S)(簡稱D)為依據,用ND記錄D中每個斷點的重數。新產生的公共鄰接與D中的元素需一一對應,避免出現冗余公共鄰接。

新產生的公共鄰接根據缺失基因的插入位置的不同可以分為兩類:

(1)外鄰接:Ci=cistart…ciend是S中的contig,x屬于S的缺失基因集合X,若xy可以構成公共鄰接,且y=cistart或y=ciend,則稱xy為外鄰接。

(2)內鄰接:除外鄰接以外的其他新增公共鄰接稱為內鄰接。

缺失基因串根據其插入S后所產生的公共鄰接數量可以分為三類:

n-typeⅠ串:由n個缺失基因組成的字符串(記做Len-n串)插入S中,新增n+1個公共鄰接,將該基因串稱為n-typeⅠ串。顯然n-typeⅠ串只能插入形如ci end]?[ci+1start的位點,該類位點記為slotⅠ。

n-typeⅡ串:Len-n串插入S中,共新增n個公共鄰接,就將該基因串稱為n-typeⅡ串。顯然n-typeⅡ串只能插入形如?[cistart或ciend]?的位點,該類位點記為slotⅡ。

n-typeⅢ串:Len-n串插入S中,共新增n-1個公共鄰接,就將該基因串稱為n-typeⅢ串。

如圖2所示,Len-5串deegg插入S,新增cd,de,ee,eg,gg,g4六個公共鄰接,故deegg是5-typeⅠ串,c]?[4為slotⅠ型位點。Len-1串a,a,b插入S,分別得到a2,a3,ab各一個公共鄰接,a2和2a記為a2,故a,a,b都是1-typeⅡ串,?[a,2]?,3]?都是slotⅡ型位點。在新增公共鄰接中a2,a3,ab,cd,g4為外鄰接,de,ee,eg,gg為內鄰接。

圖2 Contig-One-Sided-SF-max問題示例

2.1 算法核心思想框架

為方便理解,算法的大致輪廓如圖3所示。

圖3 算法核心思想

2.2 算法初始化

對實例(G,S)進行初始化,得到S中的缺失基因集合X和斷點集合D,具體算法如算法1所示。

算法1:Init(G,S)。

輸入:完整基因序列G,基因框架S

輸出:缺失基因集合X,斷點集合D

1X←?,D←?;

2 fora∈M(G) do

3 ifNM(G)[a]>NM(S)[a] anda?X then

4 添加NM(G)[a]-NM(S)[a]個atoX;

5 forab∈P(G) do

6 ifNP(G)[ab]>NP(S)[ab] andab?Dthen

7 添加NP(G)[ab]-NP(S)[ab]個abtoD;

8 returnXandD

2.3 算法實現

下面對算法每一步做出詳細介紹:

步驟一:采用貪婪策略插入1-typeⅠ串。

對于X中的元素a,從左向右掃描框架S,發現有a可插入的slotⅠ型位點,便將a插入S,將a從X中移除一個,更新S中的基因位點并將新增的兩個公共鄰接從D中刪除。具體算法如算法2所示。

算法2:Greedy-one(X,D,S)。

輸入:缺失基因集合X,斷點集合D,基因框架S=?[C1]?[C2]?…?[Cn]?

輸出:更新后的缺失基因集合X,斷點集合D,基因框架S

1 fora∈Xdo

2 forCi=cistart…ciend∈Sdo

3 ifi

4 (ciend≠ci+1startor (ciend=ci+1startand

5ND[cienda]≥2)) thenX←X-{a},

6Ci←cistart…cienda,Ci=Ci∪Ci+1,

7D←D-{cienda,aci+1start},

8ND[cienda]←ND[cienda]-1,

9ND[aci+1start]←ND[aci+1start]-1;

10 returnXandDandS

對于圖2的例子,這一步的操作是將d插入到c]?[4中,與2-近似算法結果一致。

步驟二:采用最大匹配插入1-typeⅡ串。

設SⅡ={s1,s2,…,s2n}為S中slotⅡ位點集合,n等于S中contig的數量,SⅡe={e1,e2,…,e2n},ei為位點si處相鄰的基因,構造集合X,D,SⅡ之間輔助圖Γ,EXS記錄X中元素與SⅡ中元素之間的邊,ESD記錄SⅡ中元素與D中元素之間的邊。

構造Γ后,從中計算出匹配結果M1,M1是多個三元組(a,si,xy)的集合,根據M1來插入1-typeⅡ串和更新斷點集合D。具體算法如算法3所示。

算法3:Max-match(X,D,S,SⅡ,SⅡe)。

輸入:缺失基因集合X,斷點集合D,基因框架S,slotⅡ位點集合SⅡ,SⅡe

輸出:更新后的缺失基因集合X,斷點集合D,基因框架S

1EXS←?,ESD←?,V←?,E←?;

2 newSⅡ←?,newSⅡe←?;

3D′←D中由X和SⅡe元素組成的斷點;

4 fora∈Xdo

5 forsi∈SⅡ do

6 ifasi∈D'thenEXS[a][si]←1;

7 elseEXS[a][si]←0;

8 forsi∈SⅡ do

9 forxy∈D'do

10 ifsi=xandy∈X andEXS[y][si]=1

11 thenESD[y][si][xy]←1;

12 elseESD[y][si][xy]←0;

13V←S∪SⅡ∪D',E←EXS∪ESD;

14 求出Γ(V,E)上匹配結果M;

15 forM[i]∈Mdo

16 將M[i][0]插入到S中的M[i][1]位點,

17 newSⅡ←newSⅡ+M[i][1],

18 newSⅡe←newSⅡe+M[i][0],

19X←X-{M[i][0]},D←D-{M[i][2]};

20 returnXandDandS

對于圖2的例子,2-近似算法在這一步是將缺失基因a,a,b與slotⅡ位點做最大匹配,如圖4所示,經過二分匹配后將a插入到b]?,將另一個a插入到?[b,b插入到?[a,三個基因插入后得到三個基因對ab,而斷點集合中只有一個斷點ab,故經過此步操作只能得到一個新增公共鄰接ab。

圖4 基于缺失基因、基因位點的匹配算法

文中算法在這步將缺失基因a,a,b與slotⅡ位點、斷點ab,3a,2a之間做最大匹配,如圖5所示,匹配結果M1={(a,s6,3a),(a,s4,2a),(b,s1,ab)},根據M1,將a插入到3]?,另一個a插入到2]?,b插入到?[a,這樣便得到2a,3a,ab三個新增公共鄰接。

圖5 基于缺失基因、基因位點、斷點的匹配算法

步驟三:采用最大匹配,在上一步更新的位點上插入內鄰接。

與步驟二相似,這一步仍在缺失基因、位點和斷點之間做最大匹配,SⅡ只取上一步更新得到的位點即可,調用Max-match(X,D,S,newSⅡ,newSⅡe)算法完成內鄰接的插入。

對于圖2的例子,X={e,e,g,g},D={a2,a2,a3,a3,ab,b1,b5,b5,cd,de,d1,d4,ee,eg,gg,g4,61},步驟二中更新的位點為?[b,a]?,a]?,顯然D中與缺失基因存在邊的所有斷點都不能和這些位點構成邊,無法構圖進行匹配,故跳過此步。

步驟四:采用貪婪策略,在X中匹配內鄰接對,將獲得的內鄰接對作為新contig置入集合S+。

對于X中的元素a,依次掃描其后面的其他元素,若有元素b,滿足ab是D中的斷點,則將ab作為長度為2的新contig置于S+中,并從D中移除一個斷點ab。具體算法如算法4所示。

算法4:Greedy-two(X,D)。

輸入:缺失基因集合X,斷點集合D

輸出:更新后的缺失基因集合X,斷點集合D,新contig集合S+

1i←1,S+←?;

2 fora∈Xdo

3 forb∈X-{a} do

4 ifab∈Dthen

5Ci←ab,S+←S+∪Ci,

6i←i+1,D←D-{ab};

7 returnXandDandS+

對于圖2的例子X={e,e,g,g},實際計算中X中元素順序并不確定,故經過步驟四會有兩種可能:

(1)先得到內鄰接ee,然后又得到內鄰接gg,所以S+=?[ee]?[gg]?。

(2)先得到內鄰接eg,剩下缺失基因只有e和g,但是D中已經沒有斷點eg,所以S+=?[eg]?。

步驟五:采用最大匹配在S+中插入1-typeⅡ串。

調用算法Max-match(X,D,S+,S+Ⅱ,S+Ⅱe)得到填充了1-typeⅡ串的S+,下面對步驟四中兩種可能的結果繼續進行討論:

(1)若步驟四中得到的S+=?[ee]?[gg]?,則X為空,此步無任何操作。

(2)若步驟四中得到的S+=?[eg]?,調用算法Max-match(X,D,S+,S+Ⅱ,S+Ⅱe),得到的匹配結果為M2={(e,s1,ee),(g,s2,gg)},填充后S+=?e[eg]g?。

步驟六:合并S和S+為S',X中的剩余基因,在不破壞現有鄰接情況下,可插入S'中任意位點。

步驟七:輸出S'。

可以將步驟三～步驟五合并稱為匹配內鄰接。在2-近似算法中也有相關操作,該算法是將插入1-typeⅠ串后的剩余所有基因(包括二分匹配后插入的基因)構造多重圖,在多重圖上計算最大匹配,對于圖2中的例子,以{a,a,b,e,e,g,g}為頂點構造多圖,斷點有{a2,a3,b1,b5,b5,de,d1,ee,eg,gg,g4,61}。

如圖6所示,通過構造多重圖獲得匹配結果是兩對基因對eg和eg,而斷點中只有一個eg,所以2-近似算法匹配內鄰接只新增一個公共鄰接,最終得到S'=?b[adb]a?a[bc]d[452]a?[6113]a?eg?eg?,總計新增ab,cd,d4,eg共4個公共鄰接。在最優解S*=?b[adb]?[bc]deegg[452]a?[6113]a?中共有9個新增公共鄰接,近似比為9/4,出現這種情況的原因就是在最大匹配內鄰接和外鄰接時,都出現了冗余公共鄰接現象,從而影響了該算法的近似比。

圖6 2-近似算法中多重圖最大匹配

對于圖2的例子,文中算法結果S'=?b[adb]?[bc]d[452]a?[6113]a?[ee]?[gg]?或S'=?b[adb]?[bc]d[452]a?[6113]a?e[eg]g?,分別新增7個和8個公共鄰接,近似比小于2。

2.57-近似算法同樣考慮到了避免冗余公共鄰接的問題,通過在缺失基因和contig端點基因之間構造輔助圖的方式消除冗余邊,求出最大匹配來填充基因。但是該算法對于匹配結果與缺失基因和位點之間的對應關系沒有做出判斷。對于圖2的例子,依據2.57-近似算法,可以得到最大匹配結果中的(a,b)這樣一組匹配,但是依據(a,b)是將a插入?[b,還是將b插入a]?并沒有明斷判斷,若是前者將導致b無法插入框架,若是后者則a和b都可以插入。由此可見2.57-近似算法雖然做到了避免冗余公共鄰接,但是對于外鄰接處理不佳,導致近似比大于2。

2.4 算法的性質

在本節,通過對比最優解中和文中算法所求近似解中各類型串新增公共鄰接數,證明算法近似性能比為2。設OPT為最優解中的新增公共鄰接數,OPTΙ,OPTⅡ,OPTⅢ分別為最優解中由typeⅠ,typeⅡ,typeⅢ類型基因串新增公共鄰接數,可得:

OPT=OPTΙ+OPTⅡ+OPTⅢ

(1)

設APP為文中算法所求近似解中新增公共鄰接數,APPΙ,APPⅡ,APPⅢ分別為近似解中由typeⅠ,typeⅡ,typeⅢ類型基因串新增公共鄰接數,可得:

APP=APPΙ+APPⅡ+APPⅢ

(2)

文中提到的“新增公共鄰接比”指一條或多條基因串在最優解中新增公共鄰接數與在文中算法中新增公共鄰接數之比。設某最優解中的Len-k串經過步驟四、五得到的內鄰接數為nl(k)。

證明:假設k個缺失基因在最優解中可以組成長度為k的連續基因串,則這些缺失基因通過最大匹配可得到?k/2」個匹配,即?k/2」個內鄰接。在步驟四中,使用的是貪婪策略,在最壞情況下,最大匹配中兩個匹配對錯位匹配,會使得兩個基因無法匹配,如:斷點有ac,ab,cd,缺失基因有a,b,c,d,最大匹配為{ab,cd},若ac匹配到一起,則b,d無法匹配,稱ac為錯位匹配。

因此,在步驟四使用貪婪策略可得到m個內鄰接,m的取值范圍是:

k=1和2時顯然成立。

綜上所述,引理1成立。

證明:步驟一中使用貪婪策略插入1-typeⅠ串,在最壞情況下,一個錯位的1-typeⅠ串s1,會導致一個u-typeⅠ串su變為typeⅢ,一個v-typeⅠ串sv和一個w-typeⅠ串sw變為typeⅡ,u,v,w表示各類型串的長度。

用OPT(s),APP(s)分別表示基因s在最優解和近似解中獲得的公共鄰接數,則s1,su,sv,sw新增公共鄰接比為:

u=1時,su在最優解中可以得到兩個外鄰接,在本算法中變成1-typeⅢ串無新增公共鄰接,所以

2≤u時,su在步驟四、五得到nl(u)個內鄰接,所以

1≤v,w時,sv和sw在步驟二均可獲得一個外鄰接。2≤v,w≤3時,sv和sw在步驟二獲得外鄰接后,在步驟三又均獲得一個內鄰接。

故1≤v,w≤3時,sv,sw新增公共鄰接比如下:

4≤v,w時,sv和sw在步驟三結束后都得到兩個公共鄰接,剩余長度為v-2和w-2的串在步驟四、五獲得的公共鄰接數分別為nl(v-2),nl(w-2),所以

最優解中類型為r-typeⅠ串且沒有受到錯位插入影響的基因串sr,在步驟二后均可以獲得兩個外鄰接,若r≥4在步驟三后又可以獲得兩個內鄰接。

易知r≤7時,新增公共鄰接比≤2成立。

當r≥8時,新增公共鄰接比

綜上所述,引理2得證。

證明:最優解中l-typeⅡ串,在步驟二獲得一個外鄰接,若l≥2在步驟三又可以獲得一個內鄰接。

易知,1≤l≤3時,新增公共鄰接比≤2顯然成立。

當l≥4時,新增公共鄰接比

綜上所述,引理3得證。

證明:由引理1易證得。

證明:

(3)

文中算法的時間復雜度主要取決于在X,SⅡ,D集合元素之間求最大匹配,最壞情況下集合的X每個元素,要考慮每條與SⅡ中元素構成的無向邊,接著還要考慮每條由相連的SⅡ中元素與D中元素構成的無向邊,這需要O(|X||SΙΙ||D|)時間。所以,算法時間復雜度為O(nml),其中n,m,l分別指缺失基因個數、兩倍的contig數量以及斷點個數。

3 算法的程序實現

完成Contig-One-sided-SF-max問題的多項式時間近似算法后,用Python語言實現了該算法,并編寫了具有可視化界面的Contig-One-sided-SF-max問題的程序,增強人機交互體驗。

程序主界面如圖7所示。

用戶在程序輸入框中輸入G和S,點擊確定,程序將根據文中設計的算法自動完成基因組框架填充,并給出填充結果和新增公共鄰接總數。

程序的執行如圖8所示。

圖8 程序運行界面

4 結束語

該文重新審視了基于contig的單面含重復基因的基因組框架填充問題。在設計算法時以缺失基因、基因位點、斷點的對應關系為依據(而不像以往研究中,只關注于其中兩個之間的對應),解決了2-近似算法的冗余公共鄰接問題,又解決了2.57-近似算法精確度較低的問題。對于近似比的優化方面,還需要進一步研究,希望能夠吸引更多人共同學習來完善這一問題。