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

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

0 引 言

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

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

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

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

1 相關(guān)概念

基因通常用一個正整數(shù)來表示，如某生物中的基因序列[27-29]可表示為(1,2,3,5,8)，本文為了避免使用較大整數(shù)，也用一些英文字符來表示。假設(shè)所有的基因和基因組都是無符號的，可將結(jié)果推廣到有符號的基因組。給定一個基因集合Σ，如果Σ中元素在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個基因的子串叫做i-串，一個2-串通常也叫作一個匹配對。一個匹配對ab和一個匹配對ba是相等的。給定一個不完整的序列A=a1a2a3…an，令PA={a1a2,a2a3,…,an-1an}作為A的匹配對的集合。下面將具體給出鄰接、斷點、OSSF-MNSA和OSSF-max的定義。

給定兩個序列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的最大匹配的對中，具有匹配關(guān)系的對aiai+1稱為A中相對于B的一個公共鄰接，不具有匹配關(guān)系的對ajaj+1稱為A中相對于B的一個斷點。

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

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

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

定義1最大化鄰接的單面片段填充問題(OSSF-MNSA)。輸入：一個完整基因組A，一個不完整的基因組片段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ù)據(jù)集中的基因組片段通常由一組連續(xù)的片段重疊群(contig)構(gòu)成。contig是作為基因中的字符串存在，這個字符串內(nèi)部不允許被添加任何的元素，一個基因組片段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)。輸入：給定一個完整的基因組A和一個片段B=〈C1,C2,…,Cm〉，其中A和contigCi都是一個基因集合Σ，多重集X=c(A)-c(B)≠?。問題：找到一個B*∈B+X，使得|a(B*,A)|最大化。

2 OSSF-MNSA問題

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

2.1 1.25-近似算法

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

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

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

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

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

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

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

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

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

圖2 二分圖舉例

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

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

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

2.2 1.2-近似算法

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

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

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

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

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

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

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

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

3 OSSF-max問題

本文除了研究OSSF-MNSA問題，還對OSSF-max問題中不含重復基因的多項式時間可解算法和含重復基因的2-近似算法進行研究、分析和比較。

3.1 不含重復基因的OSSF-max問題多項式時間可解算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 語

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

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

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

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