基于邻接的单面基因组片段填充问题研究进展

李春良 宋卫星 徐勤业 贾瀚栋 李晓峰 柳 楠

(山东建筑大学计算机科学与技术学院 山东 济南 250101)

0 引 言

近年来，随着基因测序技术[1-3]不断发展，基因组测序的成本已大大降低，但仅依靠基因测序技术得到一个完整的基因组序列仍是困难的。而许多研究和应用领域仍需要完整的基因组序列，如计算两个基因组之间的最小反转距离。因此，人们围绕由计算机协助完成不完整基因组的片段填充问题展开了一系列的研究。

基于无重复基因的多染色体基因组，Munoz等[4]首先提出了基因组片段填充问题，给定一个完整基因组序列A和一个不完整的基因组序列B，将缺失基因插入B中得到B′，计算A和B′的距离。其中的距离包括断点距离和二次切割与连接(DCJ)距离[5]。Yancopoulos等[5]提出了基于这两种距离的线性多项式时间精确算法。此后，文献[6]使用简单的断点距离作为相似性度量，该问题也被证明为多项式时间可解的。

当填充基因组中包含重复基因时，这个问题的解决变得困难，因为缺失基因插入位置的选择更多。基于断点距离[7-11]、基因组抽样距离[12-13]和最小公共字符串划分距离[14]这三种类型距离的片段填充问题已被证明是NP完全问题。基于基因抽样距离的片段填充问题甚至是不可近似的[15-16]。有重复基因的基因组片段填充问题可以设计出具有近似性能比的算法。本文重点讨论最大化邻接的单面基因组片段填充问题。其中，单面无重复基因组片段填充问题已经被证明是多项式时间可解的[6-7,17]；当包含了重复基因时，文献[6,17]证明该问题是NP完全的，实现4/3-近似算法；Liu等[18]利用局部搜索等策略，将近似比提高到1.25；文献[19-20]利用非盲局部搜索，将近似比进一步提高到1.2。

伴随片段重叠群(contig)可由越来越多的成熟工具(如Celera Assembler[21])计算获得，许多标准基因数据的形式也多由片段重叠群序列构成。人们也开始了对该类问题的探索[4,6]。基于片段重叠群的基因组片段填充问题，输入的基因片段构成基本单位不再是单个基因，而是片段重叠群，因此缺失基因的插入位置受到限制，只能插入片段重叠群之间，从而增加了问题的复杂性。之前的基因组片段填充问题是片段填充群由单个基因构成的情况。Zhu[22]对基于片段重叠群的问题进行了初步研究。Liu等[23]提出了基于片段重叠群的单面片段填充问题的多项式时间算法。文献[24-25]利用哈密顿路变换证明含有重复基因的单面该类问题是NP难的，并提出了近似性能比为2的近似算法和基因重复度为d的FPT算法[24]。随后Bulteau等[26]针对单面该类问题，基于最大化邻接距离和最小化断点距离分别提出了k-Mer片段填充的固定参数化算法。

1 相关概念

基因通常用一个正整数来表示，如某生物中的基因序列[27-29]可表示为(1,2,3,5,8)，本文为了避免使用较大整数，也用一些英文字符来表示。假设所有的基因和基因组都是无符号的，可将结果推广到有符号的基因组。给定一个基因集合Σ，如果Σ中元素在P中只出现一次，则P被称为排列(permutation)，用c(P)表示排列P中的元素集合；如果某些基因在A中出现多次，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的最大匹配的对中，具有匹配关系的对aiai+1称为A中相对于B的一个公共邻接，不具有匹配关系的对ajaj+1称为A中相对于B的一个断点。

由上述定义可见，序列A和B包含相同的邻接集合，但有不同的断点。A(或B)中最大匹配的对形成了A和B之间公共邻接集，用a(A,B)表示。用bA(A,B)和bB(A,B)分别表示A和B中的断点集合。以上定义如图1所示。

图1 邻接、断点等相关定义举例

在填充过程中，还对填充的字符串进行分类定义。如果插入i-串新产生i+1个邻接，那么这个i-串称为i-type-1串。如果插入i-串新产生i个邻接，那么这个i-串称为i-type-2串。如果插入i-串新产生i-1个邻接，那么这个i-串称为i-type-3串。下面给出相关研究问题的定义。

定义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问题中，发现真实数据集中的基因组片段通常由一组连续的片段重叠群(contig)构成。contig是作为基因中的字符串存在，这个字符串内部不允许被添加任何的元素，一个基因组片段B是一系列连续contig组成，定义c(B)=c(C1)∪c(C2)∪…∪c(Cm)。给定一组基因X，设B+X是将X中所有元素插入B，保证B中的contig内部不发生改变，所有可能结果排列的集合。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]中，定义出现在断点中的基因为断点基因bp-gene。此外还根据断点中两个bp-gene与缺失基因集合X的隶属关系，将断点集合分为三类：①BP1(A)：断点中一个bp-gene∈X，另一个bp-gene∉X。②BP2(A)：断点中两个bp-gene∈X。③BP3(A)：断点中两个bp-gene∉X。为了公平对待Scaffold中的每个基因，在每个输入序列的每个端点处添加封闭符号“#”，可以保证原端点元素与中间元素都具有两个邻居来构成邻接或者断点。算法的主要步骤如下：

(1) 设不完整序列B中的断点集合为β(B,A)={y1z1,y2z2,…,ymzm}，缺失基因集合X={x1,x2,…,xn}，构造二分图G1={X,β(B,A),E}，边(xj,yizi)∈E，当且仅当xj插入到yizi构成一个1-type-1串时，G1中对应节点间建立一条边，在G1中求最大匹配，按照最大匹配的结果将元素插入到匹配的断点处。

(2) 完成步骤(1)之后，在剩余待插入基因集合中，使用贪心算法找到2-type-1串，设找到的所有2-type-1串的集合为Q，将所有2-type-1串插入相对应的断点处，再使用局部搜索方法对插入方案进行优化，使得插入更多的2-type-1串，若存在下列情况，则进行优化：① 若xixj插入B中某断点ypyp+1构成2-type-1串，且xixj中的基因xi和xj可分别与待插入基因集合中的其他两个基因构成2-type-1串xixk、xjxl，则将xixj从序列和集合中删除，将xixk和xjxl插入相应断点并添加到集合Q。② 若xixj插入I中某断点ypyp+1构成2-type-1串，同时另一个串xrxs也与ypyp+1构成2-type-1串，且xixj中的一个基因xi与待插入集合中的基因xt构成2-type-1串，则将xixj从序列和集合Q中删除，将xixt插入相应断点，将xrxs插入断点ypyp+1，并将xixt、xrxs加入集合Q。

(3) 完成步骤(1)和步骤(2)后，在剩余待插入基因集合中，使用贪心算法找到3-type-1串，插入相对应的断点处。

(4) 完成步骤(1)-步骤(3)后，对于X中剩余的缺失基因没有严格的插入方法，只要保证每个插入的基因至少产生一个新邻接就可以。

4/3-近似算法没有对长度大于2的type-1串的插入进行研究，文献[17]首先根据算法插入的type-1串的数量与产生新邻接的数量分析得到更好的近似下界：

首先分析1-type-1串的数量，该算法插入的1-type-1串数量不少于任意最优解中的1-type-1串数量。算法得到的1-type-1串数量为：

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

然后分析2-type-1串的数量，涉及二分图最大匹配以及局部搜索优化。如图2所示，在二分图中，每个顶点的度不超过2，图中包含路径、圈和孤立顶点。

图2 二分图举例

最终得到的2-type-1串数量为：b′2=b′21+b′22+…+b′29+|Q∩P|。

分析3-type-1串的数量，经过贪心算法插入3-type-1串，得到3-type-1串数量为：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]通过非盲的局部搜索技术来提高近似性能比，该技术使用一种新的目标函数，由一个、两个、三个和四个基因缺失的基因串的权值和而不是直观的邻接数表示。算法首先定义了好串和不好串，如果从B′中移除这个i-串，使B′变成B+，若|a(B′,A)|-|a(B+,A)|=i+1，那么B′中的这个缺失i-串是好串；若|a(B′,A)|-|a(B+,A)|=i，则此串是不好串。

局部搜索技术需要一个目标函数来量化其解的最优程度，因此在i=1,2,3,4的情况下，用i-串数目的权值和来设置局部搜索目标函数。算法的主要步骤如下：

(1) 先插入串，如果插入一个串，该串包含X-Y中的所有基因，这样使D(B′)增加一个正数，最后，没有其他串可以被插入来增加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′变为与D(B′)相同目标函数值的基因组片段。另一个对B′的替换完成，只接受可采纳的字符串替换，每个包含两个子替换，第一个发生在长度相同的好串上，第二个必须如步骤(3)的情形。

按照步骤(1)-步骤(4)重复查找并实现一个字符串替换，只有当i没有替换时，才可以用i+1来替换改善B′。

B′中一个好的i串最多破坏B*中i+1个好串，B*中一个好的i-串可以被B′中最多i+1个好串破坏。为了比较B′和B*中好串的数量，建立二分图G=(L,R,E)，并且为了不等式的证明，尝试删除G中的一些边，简化为R中每个顶点最多有两条边，即顶点的度小于等于2。经过推导，证得近似性能比为1.2。该算法主要数据结构是最大匹配，其运行时间是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-1，j+2∈X，在Ci最后附加j-1，在Ci+1开始加上j+2，更新X和相应的红色子串，如果j-1、j+2最多其中一个在X中，那么锁上槽〈βi,αi+1〉，不允许其他基因插入此槽。

(4) 为所有剩余红色子串和可能未锁定的槽建立二分图G，它们可以被定位到type-2。对于G的每个连通分量，使用标准方法去计算最大匹配并插入相应的type-2红色子串。

(5) 在不破坏现有邻接的情况下，将剩余的type-3红色子串插到最左边或最右边的槽中。

该算法第1步时间复杂度为O(n2)，扫描A和B基因，A-B基因涂成红色，确定B中αi和βi的关系，可以识别A中type-1串。在步骤(2)、步骤(3)、步骤(5)中是线性时间可完成。步骤(4)需要花费O(n2.5)时间来计算最大匹配。所以得到以下结果：不含重复基因的OSSF-max问题可以在O(n2.5)时间内解决。

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中一个潜在的内部邻接(忽略已经与步骤(1)和步骤(2)计算相匹配的元素)，则在所有x∈X和y∈X之间存在一条边。在Q中计算最大匹配M。对于M中所有的匹配对xy，x是步骤(2)插到一端的元素，在x之前或者x之后相应地插入y。对于M中其余的匹配对，在不破坏现有邻接的前提下，插在任意没有锁定的槽里。

(4) 不破坏现有邻接的情况下，在B中任意未锁定的槽里插入X中剩余元素。

在步骤(1)之后，1-type-1串可将最优解中的1-type-1串改变为type-3，可以改变两个type-1串为type-2串，因此可得近似性能比为2。

在步骤(3)，一般图最大匹配的大小满足：

b′11是在步骤1获得的1-type-1串数量，然后通过步骤3产生的邻接数。满足：

该算法运用贪心算法、二分图最大匹配和一般图最大匹配对缺失基因进行填充，得到近似性能比为2。算法运行时间主要是顶点数为n的二分图最大匹配和一般图最大匹配，均为O(n2.5)时间复杂度，因此该算法运行时间最少为O(n2.5)。

4 结 语

本文对最大化邻接的单面基因组片段填充算法进行研究、分析和比较，将结果总结如表1所示。

表1 OSSF-MNSA和OSSF-max问题比较结果

随着快速基因测序技术的发展，基因组片段填充将会更加广泛、科学、有效并充分地运用于计算基因组学及相关领域中，保持基因组序列的完整性和准确性，从而节省全基因组生物测序成本。基因组片段填充为生物信息问题的分析和研究做了铺垫，并已经被大量的生物学[30-31]实验数据所验证。本文首先介绍了基因组片段填充技术的发展历程；然后对基于普通序列的最大化邻接距离的片段填充问题和基于片段重叠群的该问题进行了深入研究，给出了相关的基本概念，阐述了算法的主要思想；最后对已有算法及算法复杂性进行了详细的分析与比较。

基于本文的综合研究、分析和比较，可以发现以下研究趋势：(1) 0SSF-MNSA问题的近似性能比目前已经达到6/5，通过研究改善继续提高这一比例。(2) 0SSF-max问题的近似性能比目前已经达到2，通过研究改善继续提高这一比例。(3) 不含重复基因基于contig的双面基因组片段填充问题，证明是多项式时间可解问题和NP完全问题，并完成相应算法的设计。(4) 含重复基因基于contig的双面基因组片段填充问题，证明是多项式时间可解问题和NP完全问题，并完成相应算法的设计。(5) 对性能较好的算法进行实现，设计并开发相应的实验平台。