基于鄰接的單面基因組片段填充問題研究進(jìn)展

李春良 宋衛(wèi)星 徐勤業(yè) 賈瀚棟 李曉峰 柳 楠

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

0 引 言

近年來(lái)，隨著基因測(cè)序技術(shù)[1-3]不斷發(fā)展，基因組測(cè)序的成本已大大降低，但僅依靠基因測(cè)序技術(shù)得到一個(gè)完整的基因組序列仍是困難的。而許多研究和應(yīng)用領(lǐng)域仍需要完整的基因組序列，如計(jì)算兩個(gè)基因組之間的最小反轉(zhuǎn)距離。因此，人們圍繞由計(jì)算機(jī)協(xié)助完成不完整基因組的片段填充問題展開了一系列的研究。

基于無(wú)重復(fù)基因的多染色體基因組，Munoz等[4]首先提出了基因組片段填充問題，給定一個(gè)完整基因組序列A和一個(gè)不完整的基因組序列B，將缺失基因插入B中得到B′，計(jì)算A和B′的距離。其中的距離包括斷點(diǎn)距離和二次切割與連接(DCJ)距離[5]。Yancopoulos等[5]提出了基于這兩種距離的線性多項(xiàng)式時(shí)間精確算法。此后，文獻(xiàn)[6]使用簡(jiǎn)單的斷點(diǎn)距離作為相似性度量，該問題也被證明為多項(xiàng)式時(shí)間可解的。

當(dāng)填充基因組中包含重復(fù)基因時(shí)，這個(gè)問題的解決變得困難，因?yàn)槿笔Щ虿迦胛恢玫倪x擇更多?；跀帱c(diǎn)距離[7-11]、基因組抽樣距離[12-13]和最小公共字符串劃分距離[14]這三種類型距離的片段填充問題已被證明是NP完全問題?；诨虺闃泳嚯x的片段填充問題甚至是不可近似的[15-16]。有重復(fù)基因的基因組片段填充問題可以設(shè)計(jì)出具有近似性能比的算法。本文重點(diǎn)討論最大化鄰接的單面基因組片段填充問題。其中，單面無(wú)重復(fù)基因組片段填充問題已經(jīng)被證明是多項(xiàng)式時(shí)間可解的[6-7,17]；當(dāng)包含了重復(fù)基因時(shí)，文獻(xiàn)[6,17]證明該問題是NP完全的，實(shí)現(xiàn)4/3-近似算法；Liu等[18]利用局部搜索等策略，將近似比提高到1.25；文獻(xiàn)[19-20]利用非盲局部搜索，將近似比進(jìn)一步提高到1.2。

伴隨片段重疊群(contig)可由越來(lái)越多的成熟工具(如Celera Assembler[21])計(jì)算獲得，許多標(biāo)準(zhǔn)基因數(shù)據(jù)的形式也多由片段重疊群序列構(gòu)成。人們也開始了對(duì)該類問題的探索[4,6]?；谄沃丿B群的基因組片段填充問題，輸入的基因片段構(gòu)成基本單位不再是單個(gè)基因，而是片段重疊群，因此缺失基因的插入位置受到限制，只能插入片段重疊群之間，從而增加了問題的復(fù)雜性。之前的基因組片段填充問題是片段填充群由單個(gè)基因構(gòu)成的情況。Zhu[22]對(duì)基于片段重疊群的問題進(jìn)行了初步研究。Liu等[23]提出了基于片段重疊群的單面片段填充問題的多項(xiàng)式時(shí)間算法。文獻(xiàn)[24-25]利用哈密頓路變換證明含有重復(fù)基因的單面該類問題是NP難的，并提出了近似性能比為2的近似算法和基因重復(fù)度為d的FPT算法[24]。隨后Bulteau等[26]針對(duì)單面該類問題，基于最大化鄰接距離和最小化斷點(diǎn)距離分別提出了k-Mer片段填充的固定參數(shù)化算法。

1 相關(guān)概念

基因通常用一個(gè)正整數(shù)來(lái)表示，如某生物中的基因序列[27-29]可表示為(1,2,3,5,8)，本文為了避免使用較大整數(shù)，也用一些英文字符來(lái)表示。假設(shè)所有的基因和基因組都是無(wú)符號(hào)的，可將結(jié)果推廣到有符號(hào)的基因組。給定一個(gè)基因集合Σ，如果Σ中元素在P中只出現(xiàn)一次，則P被稱為排列(permutation)，用c(P)表示排列P中的元素集合；如果某些基因在A中出現(xiàn)多次，A被稱為序列(sequence)，用c(A)表示A中基因的集合，是Σ的多重子集。例如，Σ={a,b,c,d,e}，P=dace，A=abcdeace，c(P)={a,c,d,e}，c(A)={a,a,b,c,c,d,e,e}。含有i個(gè)基因的子串叫做i-串，一個(gè)2-串通常也叫作一個(gè)匹配對(duì)。一個(gè)匹配對(duì)ab和一個(gè)匹配對(duì)ba是相等的。給定一個(gè)不完整的序列A=a1a2a3…an，令PA={a1a2,a2a3,…,an-1an}作為A的匹配對(duì)的集合。下面將具體給出鄰接、斷點(diǎn)、OSSF-MNSA和OSSF-max的定義。

給定兩個(gè)序列A=a1a2…an和B=b1b2…bn，如果aiai+1=bjbj+1(或者aiai+1=bj+1bj)，其中aiai+1∈PA，bj+1bj∈PB，則稱aiai+1與bjbj+1互相匹配。在PA和PB的最大匹配的對(duì)中，具有匹配關(guān)系的對(duì)aiai+1稱為A中相對(duì)于B的一個(gè)公共鄰接，不具有匹配關(guān)系的對(duì)ajaj+1稱為A中相對(duì)于B的一個(gè)斷點(diǎn)。

由上述定義可見，序列A和B包含相同的鄰接集合，但有不同的斷點(diǎn)。A(或B)中最大匹配的對(duì)形成了A和B之間公共鄰接集，用a(A,B)表示。用bA(A,B)和bB(A,B)分別表示A和B中的斷點(diǎn)集合。以上定義如圖1所示。

圖1 鄰接、斷點(diǎn)等相關(guān)定義舉例

在填充過(guò)程中，還對(duì)填充的字符串進(jìn)行分類定義。如果插入i-串新產(chǎn)生i+1個(gè)鄰接，那么這個(gè)i-串稱為i-type-1串。如果插入i-串新產(chǎn)生i個(gè)鄰接，那么這個(gè)i-串稱為i-type-2串。如果插入i-串新產(chǎn)生i-1個(gè)鄰接，那么這個(gè)i-串稱為i-type-3串。下面給出相關(guān)研究問題的定義。

定義1最大化鄰接的單面片段填充問題(OSSF-MNSA)。輸入：一個(gè)完整基因組A，一個(gè)不完整的基因組片段B，其中X=c(A)-c(B)≠?，Y=c(B)-c(A)=?。

問題：找到一組插入操作將c(A)-c(B)的基因插入B得到B*，使得|a(A,B*)|-|a(A,B)|值最大。

在研究OSSF-MNSA問題中，發(fā)現(xiàn)真實(shí)數(shù)據(jù)集中的基因組片段通常由一組連續(xù)的片段重疊群(contig)構(gòu)成。contig是作為基因中的字符串存在，這個(gè)字符串內(nèi)部不允許被添加任何的元素，一個(gè)基因組片段B是一系列連續(xù)contig組成，定義c(B)=c(C1)∪c(C2)∪…∪c(Cm)。給定一組基因X，設(shè)B+X是將X中所有元素插入B，保證B中的contig內(nèi)部不發(fā)生改變，所有可能結(jié)果排列的集合。X中的元素只能被插入到Ci的前面或者后面。為此，本文給出基于contig的最大化鄰接單面片段填充問題的定義。

定義2基于contig的最大化鄰接單面片段填充問題(OSSF-max)。輸入：給定一個(gè)完整的基因組A和一個(gè)片段B=〈C1,C2,…,Cm〉，其中A和contigCi都是一個(gè)基因集合Σ，多重集X=c(A)-c(B)≠?。問題：找到一個(gè)B*∈B+X，使得|a(B*,A)|最大化。

2 OSSF-MNSA問題

本文針對(duì)OSSF-MNSA問題的1.25-近似算法和1.2-近似算法進(jìn)行研究、分析和比較。

2.1 1.25-近似算法

在1.25近似算法[18]中，定義出現(xiàn)在斷點(diǎn)中的基因?yàn)閿帱c(diǎn)基因bp-gene。此外還根據(jù)斷點(diǎn)中兩個(gè)bp-gene與缺失基因集合X的隸屬關(guān)系，將斷點(diǎn)集合分為三類：①BP1(A)：斷點(diǎn)中一個(gè)bp-gene∈X，另一個(gè)bp-gene?X。②BP2(A)：斷點(diǎn)中兩個(gè)bp-gene∈X。③BP3(A)：斷點(diǎn)中兩個(gè)bp-gene?X。為了公平對(duì)待Scaffold中的每個(gè)基因，在每個(gè)輸入序列的每個(gè)端點(diǎn)處添加封閉符號(hào)“#”，可以保證原端點(diǎn)元素與中間元素都具有兩個(gè)鄰居來(lái)構(gòu)成鄰接或者斷點(diǎn)。算法的主要步驟如下：

(1) 設(shè)不完整序列B中的斷點(diǎn)集合為β(B,A)={y1z1,y2z2,…,ymzm}，缺失基因集合X={x1,x2,…,xn}，構(gòu)造二分圖G1={X,β(B,A),E}，邊(xj,yizi)∈E，當(dāng)且僅當(dāng)xj插入到y(tǒng)izi構(gòu)成一個(gè)1-type-1串時(shí)，G1中對(duì)應(yīng)節(jié)點(diǎn)間建立一條邊，在G1中求最大匹配，按照最大匹配的結(jié)果將元素插入到匹配的斷點(diǎn)處。

(2) 完成步驟(1)之后，在剩余待插入基因集合中，使用貪心算法找到2-type-1串，設(shè)找到的所有2-type-1串的集合為Q，將所有2-type-1串插入相對(duì)應(yīng)的斷點(diǎn)處，再使用局部搜索方法對(duì)插入方案進(jìn)行優(yōu)化，使得插入更多的2-type-1串，若存在下列情況，則進(jìn)行優(yōu)化：① 若xixj插入B中某斷點(diǎn)ypyp+1構(gòu)成2-type-1串，且xixj中的基因xi和xj可分別與待插入基因集合中的其他兩個(gè)基因構(gòu)成2-type-1串xixk、xjxl，則將xixj從序列和集合中刪除，將xixk和xjxl插入相應(yīng)斷點(diǎn)并添加到集合Q。② 若xixj插入I中某斷點(diǎn)ypyp+1構(gòu)成2-type-1串，同時(shí)另一個(gè)串xrxs也與ypyp+1構(gòu)成2-type-1串，且xixj中的一個(gè)基因xi與待插入集合中的基因xt構(gòu)成2-type-1串，則將xixj從序列和集合Q中刪除，將xixt插入相應(yīng)斷點(diǎn)，將xrxs插入斷點(diǎn)ypyp+1，并將xixt、xrxs加入集合Q。

(3) 完成步驟(1)和步驟(2)后，在剩余待插入基因集合中，使用貪心算法找到3-type-1串，插入相對(duì)應(yīng)的斷點(diǎn)處。

(4) 完成步驟(1)-步驟(3)后，對(duì)于X中剩余的缺失基因沒有嚴(yán)格的插入方法，只要保證每個(gè)插入的基因至少產(chǎn)生一個(gè)新鄰接就可以。

4/3-近似算法沒有對(duì)長(zhǎng)度大于2的type-1串的插入進(jìn)行研究，文獻(xiàn)[17]首先根據(jù)算法插入的type-1串的數(shù)量與產(chǎn)生新鄰接的數(shù)量分析得到更好的近似下界：

首先分析1-type-1串的數(shù)量，該算法插入的1-type-1串?dāng)?shù)量不少于任意最優(yōu)解中的1-type-1串?dāng)?shù)量。算法得到的1-type-1串?dāng)?shù)量為：

b′1=b′11+b′12+b′13+…+b′19+|M∩W|

然后分析2-type-1串的數(shù)量，涉及二分圖最大匹配以及局部搜索優(yōu)化。如圖2所示，在二分圖中，每個(gè)頂點(diǎn)的度不超過(guò)2，圖中包含路徑、圈和孤立頂點(diǎn)。

圖2 二分圖舉例

最終得到的2-type-1串?dāng)?shù)量為：b′2=b′21+b′22+…+b′29+|Q∩P|。

分析3-type-1串的數(shù)量，經(jīng)過(guò)貪心算法插入3-type-1串，得到3-type-1串?dāng)?shù)量為：b′3=b′31+b′32+b′33+b′34+b′35+|Z|。

從而得到1.25-近似算法。該算法中時(shí)間復(fù)雜度最高的操作是最大匹配，其運(yùn)行時(shí)間為O(n2.5)，所以該算法時(shí)間復(fù)雜度為O(n2.5)。

2.2 1.2-近似算法

文獻(xiàn)[20]通過(guò)非盲的局部搜索技術(shù)來(lái)提高近似性能比，該技術(shù)使用一種新的目標(biāo)函數(shù)，由一個(gè)、兩個(gè)、三個(gè)和四個(gè)基因缺失的基因串的權(quán)值和而不是直觀的鄰接數(shù)表示。算法首先定義了好串和不好串，如果從B′中移除這個(gè)i-串，使B′變成B+，若|a(B′,A)|-|a(B+,A)|=i+1，那么B′中的這個(gè)缺失i-串是好串；若|a(B′,A)|-|a(B+,A)|=i，則此串是不好串。

局部搜索技術(shù)需要一個(gè)目標(biāo)函數(shù)來(lái)量化其解的最優(yōu)程度，因此在i=1,2,3,4的情況下，用i-串?dāng)?shù)目的權(quán)值和來(lái)設(shè)置局部搜索目標(biāo)函數(shù)。算法的主要步驟如下：

(1) 先插入串，如果插入一個(gè)串，該串包含X-Y中的所有基因，這樣使D(B′)增加一個(gè)正數(shù)，最后，沒有其他串可以被插入來(lái)增加D(B′)。

(2) 如果從B′中刪除一個(gè)好串，使得B′變成B+，則將2-串插入B+得到B++，其中D(B++)>D(B′)。之后沒有2-串可以替代B′中的一個(gè)好串從而增加D(B′)。

(3) 添加較短的好串。如果從B′中刪除一個(gè)好串，使得B′變成B+，插入到B′中的其他串，短于從B′中刪除的串，使得B+變成B++，其中D(B++)>D(B′)。至少用1-串去替換1-串，可以允許1-串去替換1-串、2-串、3-串和4-串。最后，沒有1-串可以替換B′中的1-串、2-串、3-串和4-串來(lái)增加D(B′)。

(4) 一個(gè)串可以替換長(zhǎng)度相同的好串，把B′變?yōu)榕cD(B′)相同目標(biāo)函數(shù)值的基因組片段。另一個(gè)對(duì)B′的替換完成，只接受可采納的字符串替換，每個(gè)包含兩個(gè)子替換，第一個(gè)發(fā)生在長(zhǎng)度相同的好串上，第二個(gè)必須如步驟(3)的情形。

按照步驟(1)-步驟(4)重復(fù)查找并實(shí)現(xiàn)一個(gè)字符串替換，只有當(dāng)i沒有替換時(shí)，才可以用i+1來(lái)替換改善B′。

B′中一個(gè)好的i串最多破壞B*中i+1個(gè)好串，B*中一個(gè)好的i-串可以被B′中最多i+1個(gè)好串破壞。為了比較B′和B*中好串的數(shù)量，建立二分圖G=(L,R,E)，并且為了不等式的證明，嘗試刪除G中的一些邊，簡(jiǎn)化為R中每個(gè)頂點(diǎn)最多有兩條邊，即頂點(diǎn)的度小于等于2。經(jīng)過(guò)推導(dǎo)，證得近似性能比為1.2。該算法主要數(shù)據(jù)結(jié)構(gòu)是最大匹配，其運(yùn)行時(shí)間是O(n2.5)，故該算法時(shí)間復(fù)雜度是O(n2.5)。

3 OSSF-max問題

本文除了研究OSSF-MNSA問題，還對(duì)OSSF-max問題中不含重復(fù)基因的多項(xiàng)式時(shí)間可解算法和含重復(fù)基因的2-近似算法進(jìn)行研究、分析和比較。

3.1 不含重復(fù)基因的OSSF-max問題多項(xiàng)式時(shí)間可解算法

(1) 在A串中，將所有缺失基因標(biāo)記為紅色，識(shí)別所有紅色子串。

(2) 將所有type-1紅色子串插入Ci和Ci+1之間正確的槽里，鎖上槽，不允許其他基因插入此槽。

(3) 對(duì)于槽〈βi,αi+1〉，如果βiαi+1=〈j,j+1〉已經(jīng)有一個(gè)鄰接，j-1，j+2∈X，在Ci最后附加j-1，在Ci+1開始加上j+2，更新X和相應(yīng)的紅色子串，如果j-1、j+2最多其中一個(gè)在X中，那么鎖上槽〈βi,αi+1〉，不允許其他基因插入此槽。

(4) 為所有剩余紅色子串和可能未鎖定的槽建立二分圖G，它們可以被定位到type-2。對(duì)于G的每個(gè)連通分量，使用標(biāo)準(zhǔn)方法去計(jì)算最大匹配并插入相應(yīng)的type-2紅色子串。

(5) 在不破壞現(xiàn)有鄰接的情況下，將剩余的type-3紅色子串插到最左邊或最右邊的槽中。

該算法第1步時(shí)間復(fù)雜度為O(n2)，掃描A和B基因，A-B基因涂成紅色，確定B中αi和βi的關(guān)系，可以識(shí)別A中type-1串。在步驟(2)、步驟(3)、步驟(5)中是線性時(shí)間可完成。步驟(4)需要花費(fèi)O(n2.5)時(shí)間來(lái)計(jì)算最大匹配。所以得到以下結(jié)果：不含重復(fù)基因的OSSF-max問題可以在O(n2.5)時(shí)間內(nèi)解決。

3.2 含重復(fù)基因的OSSF-max問題2-近似算法

(1) 使用貪心算法從左到右掃描B中元素，元素x被插入到槽中形成1-type-1串，生成2個(gè)新的外部鄰接，鎖定此槽。

(2) 用二分圖最大匹配識(shí)別所有1-type-2串(或者外部公共鄰接)xz，將x插到槽里并將此槽進(jìn)行更新，如果x被插到槽z°中(z之后)，將此槽更新為槽x°；如果x被插到槽°z中(z之前)，將此槽更新為槽°x。

(3) 對(duì)于第1步之后X的所有剩余元素(包括步驟(2)插入的元素)，計(jì)算一個(gè)多圖Q，其頂點(diǎn)是X中元素(步驟(1)之后的元素)，如果xy是A中一個(gè)潛在的內(nèi)部鄰接(忽略已經(jīng)與步驟(1)和步驟(2)計(jì)算相匹配的元素)，則在所有x∈X和y∈X之間存在一條邊。在Q中計(jì)算最大匹配M。對(duì)于M中所有的匹配對(duì)xy，x是步驟(2)插到一端的元素，在x之前或者x之后相應(yīng)地插入y。對(duì)于M中其余的匹配對(duì)，在不破壞現(xiàn)有鄰接的前提下，插在任意沒有鎖定的槽里。

(4) 不破壞現(xiàn)有鄰接的情況下，在B中任意未鎖定的槽里插入X中剩余元素。

在步驟(1)之后，1-type-1串可將最優(yōu)解中的1-type-1串改變?yōu)閠ype-3，可以改變兩個(gè)type-1串為type-2串，因此可得近似性能比為2。

在步驟(3)，一般圖最大匹配的大小滿足：

b′11是在步驟1獲得的1-type-1串?dāng)?shù)量，然后通過(guò)步驟3產(chǎn)生的鄰接數(shù)。滿足：

該算法運(yùn)用貪心算法、二分圖最大匹配和一般圖最大匹配對(duì)缺失基因進(jìn)行填充，得到近似性能比為2。算法運(yùn)行時(shí)間主要是頂點(diǎn)數(shù)為n的二分圖最大匹配和一般圖最大匹配，均為O(n2.5)時(shí)間復(fù)雜度，因此該算法運(yùn)行時(shí)間最少為O(n2.5)。

4 結(jié) 語(yǔ)

本文對(duì)最大化鄰接的單面基因組片段填充算法進(jìn)行研究、分析和比較，將結(jié)果總結(jié)如表1所示。

表1 OSSF-MNSA和OSSF-max問題比較結(jié)果

隨著快速基因測(cè)序技術(shù)的發(fā)展，基因組片段填充將會(huì)更加廣泛、科學(xué)、有效并充分地運(yùn)用于計(jì)算基因組學(xué)及相關(guān)領(lǐng)域中，保持基因組序列的完整性和準(zhǔn)確性，從而節(jié)省全基因組生物測(cè)序成本。基因組片段填充為生物信息問題的分析和研究做了鋪墊，并已經(jīng)被大量的生物學(xué)[30-31]實(shí)驗(yàn)數(shù)據(jù)所驗(yàn)證。本文首先介紹了基因組片段填充技術(shù)的發(fā)展歷程；然后對(duì)基于普通序列的最大化鄰接距離的片段填充問題和基于片段重疊群的該問題進(jìn)行了深入研究，給出了相關(guān)的基本概念，闡述了算法的主要思想；最后對(duì)已有算法及算法復(fù)雜性進(jìn)行了詳細(xì)的分析與比較。

基于本文的綜合研究、分析和比較，可以發(fā)現(xiàn)以下研究趨勢(shì)：(1) 0SSF-MNSA問題的近似性能比目前已經(jīng)達(dá)到6/5，通過(guò)研究改善繼續(xù)提高這一比例。(2) 0SSF-max問題的近似性能比目前已經(jīng)達(dá)到2，通過(guò)研究改善繼續(xù)提高這一比例。(3) 不含重復(fù)基因基于contig的雙面基因組片段填充問題，證明是多項(xiàng)式時(shí)間可解問題和NP完全問題，并完成相應(yīng)算法的設(shè)計(jì)。(4) 含重復(fù)基因基于contig的雙面基因組片段填充問題，證明是多項(xiàng)式時(shí)間可解問題和NP完全問題，并完成相應(yīng)算法的設(shè)計(jì)。(5) 對(duì)性能較好的算法進(jìn)行實(shí)現(xiàn)，設(shè)計(jì)并開發(fā)相應(yīng)的實(shí)驗(yàn)平臺(tái)。