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

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

(山東建筑大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,山東 濟(jì)南 250101)

0 引 言

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

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

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

1 相關(guān)定義

給定集合Σ,Σ包含所有用于標(biāo)識基因的字符。設(shè)A是由Σ中元素構(gòu)成的字符串,A中每個元素可多次出現(xiàn),稱A為序列。用M(A)表示A中的字符集合,由于M(A)可能是多重集,為方便計(jì)算,用NM(A)記錄每個字符的重?cái)?shù)。例如Σ={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中的一基因?qū)?用P(A)表示A中所有基因?qū)Φ募?P(A)也可能是多重集,用NP(A)記錄基因?qū)χ財(cái)?shù)。前面例子中P(A)={ab,b1,12,2a,a2},NP(A)={ab:1,b1:1,12:1,2a:2}。易知基因?qū)?a與a2無區(qū)別,故后面討論中xy與yx視為相同基因?qū)Α?/p>

給定兩條序列Α=a1a2…am和B=b1b2…bn,若P(A)中的一基因?qū)iai+1,在P(B)可以找到bkbk+1與之對應(yīng)(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中沒有匹配到公共鄰接的基因?qū)?這些基因?qū)ΨQ為A或B中的斷點(diǎn)。示例如圖1所示。

圖1 公共鄰接和斷點(diǎn)的示例

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

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

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

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

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

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

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

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

新產(chǎn)生的公共鄰接根據(jù)缺失基因的插入位置的不同可以分為兩類:

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

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

缺失基因串根據(jù)其插入S后所產(chǎn)生的公共鄰接數(shù)量可以分為三類:

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

n-typeⅡ串:Len-n串插入S中,共新增n個公共鄰接,就將該基因串稱為n-typeⅡ串。顯然n-typeⅡ串只能插入形如?[cistart或ciend]?的位點(diǎn),該類位點(diǎn)記為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Ⅰ型位點(diǎn)。Len-1串a(chǎn),a,b插入S,分別得到a2,a3,ab各一個公共鄰接,a2和2a記為a2,故a,a,b都是1-typeⅡ串,?[a,2]?,3]?都是slotⅡ型位點(diǎn)。在新增公共鄰接中a2,a3,ab,cd,g4為外鄰接,de,ee,eg,gg為內(nèi)鄰接。

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

2.1 算法核心思想框架

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

圖3 算法核心思想

2.2 算法初始化

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

算法1:Init(G,S)。

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

輸出:缺失基因集合X,斷點(diǎn)集合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 算法實(shí)現(xiàn)

下面對算法每一步做出詳細(xì)介紹:

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

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

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

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

輸出:更新后的缺失基因集合X,斷點(diǎn)集合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-近似算法結(jié)果一致。

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

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

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

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

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

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

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

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

3D′←D中由X和SⅡe元素組成的斷點(diǎn);

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)上匹配結(jié)果M;

15 forM[i]∈Mdo

16 將M[i][0]插入到S中的M[i][1]位點(diǎn),

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Ⅱ位點(diǎn)做最大匹配,如圖4所示,經(jīng)過二分匹配后將a插入到b]?,將另一個a插入到?[b,b插入到?[a,三個基因插入后得到三個基因?qū)b,而斷點(diǎn)集合中只有一個斷點(diǎn)ab,故經(jīng)過此步操作只能得到一個新增公共鄰接ab。

圖4 基于缺失基因、基因位點(diǎn)的匹配算法

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

圖5 基于缺失基因、基因位點(diǎn)、斷點(diǎn)的匹配算法

步驟三:采用最大匹配,在上一步更新的位點(diǎn)上插入內(nèi)鄰接。

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

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

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

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

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

輸入:缺失基因集合X,斷點(diǎn)集合D

輸出:更新后的缺失基因集合X,斷點(diǎn)集合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},實(shí)際計(jì)算中X中元素順序并不確定,故經(jīng)過步驟四會有兩種可能:

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

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

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

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

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

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

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

步驟七:輸出S'。

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

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

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

對于圖2的例子,文中算法結(jié)果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端點(diǎn)基因之間構(gòu)造輔助圖的方式消除冗余邊,求出最大匹配來填充基因。但是該算法對于匹配結(jié)果與缺失基因和位點(diǎn)之間的對應(yīng)關(guān)系沒有做出判斷。對于圖2的例子,依據(jù)2.57-近似算法,可以得到最大匹配結(jié)果中的(a,b)這樣一組匹配,但是依據(jù)(a,b)是將a插入?[b,還是將b插入a]?并沒有明斷判斷,若是前者將導(dǎo)致b無法插入框架,若是后者則a和b都可以插入。由此可見2.57-近似算法雖然做到了避免冗余公共鄰接,但是對于外鄰接處理不佳,導(dǎo)致近似比大于2。

2.4 算法的性質(zhì)

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

OPT=OPTΙ+OPTⅡ+OPTⅢ

(1)

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

APP=APPΙ+APPⅡ+APPⅢ

(2)

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

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

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

k=1和2時顯然成立。

綜上所述,引理1成立。

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

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

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

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

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

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

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

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

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

當(dāng)r≥8時,新增公共鄰接比

綜上所述,引理2得證。

證明:最優(yōu)解中l(wèi)-typeⅡ串,在步驟二獲得一個外鄰接,若l≥2在步驟三又可以獲得一個內(nèi)鄰接。

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

當(dāng)l≥4時,新增公共鄰接比

綜上所述,引理3得證。

證明:由引理1易證得。

證明:

(3)

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

3 算法的程序?qū)崿F(xiàn)

完成Contig-One-sided-SF-max問題的多項(xiàng)式時間近似算法后,用Python語言實(shí)現(xiàn)了該算法,并編寫了具有可視化界面的Contig-One-sided-SF-max問題的程序,增強(qiáng)人機(jī)交互體驗(yàn)。

程序主界面如圖7所示。

用戶在程序輸入框中輸入G和S,點(diǎn)擊確定,程序?qū)⒏鶕?jù)文中設(shè)計(jì)的算法自動完成基因組框架填充,并給出填充結(jié)果和新增公共鄰接總數(shù)。

程序的執(zhí)行如圖8所示。

圖8 程序運(yùn)行界面

4 結(jié)束語

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