求解两物种小系统发育问题的模拟退火算法

摘要：针对复制丢失比对问题模型，提出求解复制丢失演化模型下两物种小系统发育问题（SPP）的模拟退火算法（SA2SP）。SA2SP引入比对算法用于构造问题初始解；引入标记算法用于构建问题解的目标函数，以得到问题解的进化代价；同时还引入3种智能邻域函数，利用基因序列的进化特性，指导性地产生邻域解。利用4种真实菌属的核糖体核糖核酸（rRNA）和转运核糖核酸（tRNA）基因数据对算法的性能进行测试，实验结果表明， SA2SP能够获得较伪布尔线性规划（PBLP）求解算法更小的进化代价，是求解复制丢失演化模型下两物种小系统发育问题的一种有效方法。

关键词：小系统发育问题；模拟退火；邻域函数；基因

中图分类号：TP301.6文献标志码：A英文标题

0引言

高通量基因测序技术提供了大量已测定的脱氧核糖核酸（DeoxyriboNucleic Acid， DNA）序列数据，这些数据使基因家族（gene family） 进化史的研究工作成为可能，继而帮助揭示显性的基因组学基础，研究基因的功能。例如，转运核糖核酸（transfer RiboNucleic Acid， tRNA）是合成蛋白质所必不可少的重要物质，探索物种间tRNA的进化史或许可为转录机制的研究提供新的思路[1-4]。

小系统发育问题（Small Phylogeny Problem， SPP）[1-2，5-6]是求解基因家族进化史的一个重要问题模型，由于复制（duplication）及丢失（loss）这两种基因内容修改操作在基因家族进化中起主导作用[7-10]，Holloway等[1]提出一种SPP的衍生模型，即复制丢失演化模型下的两物种小系统发育问题（TwoSpecies SPP in the DuplicationLoss model， 2SPPDL），并将该问题归约成复制丢失比对问题（DuplicationLoss Alignment problem， DLA）模型。DLA模型一般情况下是 NP 难的[1]，Holloway等[1]提出了求解该模型的伪布尔线性规划（PseudoBoolean Linear Programming， PBLP）求解算法，该算法时间复杂性是指数级别的，其性能受到问题规模的制约，当基因序列位点数增多时，算法的有效性会变差。PBLP算法为确保获得无环解，循环地加入除环约束，针对该问题，2013年，Andreotti等[2]提出基于图论理论的DupLoCut算法，该算法在求解的过程中智能地剪切掉无用边，从而提高了算法的求解效率，但该算法的时间复杂性仍然是指数级别的。

模拟退火（Simulated Annealing， SA）算法作为一种智能优化方法，是求解NP完全问题的有效手段，已被成功运用于求解基因序列比对问题，并取得较好的比对效果[11-14]。本文基于SA算法对DLA模型进行研究，提出求解2SPPDL问题的模拟退火算法（Simulated Annealing algorithm for Solving 2SPPDL， SA2SP）。主要工作包括：1）引入一种基于贪心策略的启发式比对算法，并使用该算法生成问题初始解并辅助产生邻域解；2）引入一种启发式的比对标记算法，并以此计算问题解的目标函数；3）引入基因块智能移动、相邻基因块位置互换和重新匹配基因块三种智能邻域函数， 邻域函数充分利用了基因序列的进化特性，产生具有一定指导性的邻域解。实验测试结果表明，算法SA2SP能有效求解DLA模型，并获得较算法PBLP更小的进化代价，是求解2SPPDL问题的一种有效方法。

3实验结果

本文采用真实的生物数据来比较分析算法SA2SP和PBLP的性能。SA2SP算法在一台安装了Windows XP Professional操作系统的联想工作站（Intel Core i53570 CPU 3.40GHz，内存为4GB）上运行，程序编译器为Microsoft Visual C# 2012。PBLP算法在一台安装了Ubuntu 12.04操作系统的同样配置的联想工作站上运行，程序编译器为Python 2.7.3。

3.1实验数据

本文采用表1中4种菌属的真实tRNA和rRNA数据对SA2SP和PBLP两种算法进行算法性能比较：5个芽孢杆菌（Bacillus）、5个耶尔森氏菌（Yersinia）、10个链球菌（Streptococcus）和12个不动杆菌（Acinetobacter）。为方便描述，本文对各菌属的基因进行编号，详情参见表1，其中NC开头的基因序列来自于National Center of Biotechnology Information （NCBI） 网站（ http：//www.ncbi.nlm.nih.gov/），其余基因序列来自于PathoSystems Resource Integration Center （PATRIC）网站（http：//patricbrc.vbi.vt.edu/）。

3.2性能评价

本文通过进化代价和运行时间两项评价指标对算法SA2SP和PBLP进行比较分析。SA2SP算法参数设置如下：初始温度T0=8000，终止温度Tf=0，降温系数a=0.90，内循环次数N=25。表2中的基因序列组合采用表1中的基因编号组合代替。例如，耶尔森氏菌的基因编号组合{1，2}表示基因Yersinia_pestis_Pestoides_F（94）和基因Yersinia_pestis_Pestoides_A（87）的组合。

表2中，给出芽孢杆菌10种基因序列组合的实验测试结果，稳定RNA家族数为42左右，每个基因序列稳定RNAs数目在96～149。从表2可看出，算法SA2SP在各种序列组合下，均能获得较算法PBLP更小的进化代价。例如，对芽孢杆菌NC_004722 （编号为1）与芽孢杆菌NC_006274（编号为2）的基因序列组合，算法SA2SP推算其祖先基因序列所对应的进化代价为9，而PBLP推算其祖先基因序列所对应的进化代价为17，基于最大节约原则[1-2]，算法SA2SP能推算出更优的祖先基因序列。

图5～6分别通过10种链球菌的45种基因序列组合和12种不动杆菌的66种基因序列组合对算法SA2SP和PBLP的进化代价进行比较。链球菌的基因组有35个左右稳定RNA家族，每个基因组的稳定RNA数目在80～101；不动杆菌的基因组有36个左右稳定RNA家族，每个基因组的稳定RNA数目在96～115。以下各图中，x和y坐标轴为基因序列编号，z坐标轴为进化代价，其中灰色小圆点代表PBLP算法的进化代价，黑色小圆点代表SA2SP算法的进化代价。从图5～6可看出：在两种菌属的各种基因序列编号组合下，算法SA2SP均能获得较PBLP更小的进化代价，故其能推算出较PBLP更优的祖先基因序列。

表3给出在相同计算机配置下，算法SA2SP及PBLP运行各菌属基因序列所有组合的平均运行时间。虽然SA2SP算法的运行时间较PBLP算法要长，但最长时间仅为232.6s。因此，在实际应用中是切实可行的。

s

菌属名称PBLPSA2SP菌属名称PBLPSA2SP芽孢杆菌58.2111.4链球菌58.8176.1耶尔森氏菌11.046.9不动杆菌128.7232.6

上述实验数据显示，SA2SP算法能获得较PBLP算法更小的进化代价，故其能推算出较PBLP算法更优的祖先基因序列。SA2SP算法具有两个特点：1）SA2SP算法通过使用比对算法来获得初始解，初始解质量的提高有助于在合理的时间范围内获得较好的解；2）利用具有指导性的智能邻域函数产生邻域解，在保持邻域解多样性的基础上，提高了算法的寻优能力。

4结语

本文提出一种求解复制丢失演化模型下两物种小系统发育问题的模拟退火算法SA2SP。算法SA2SP引入了比对算法、标记算法及3种智能邻域函数。智能邻域函数在保证邻域解多样性的前提下，充分利用基因序列的进化特性，使算法能够获得较好的邻域解，从而快速地收敛到最优解。利用4种菌属的稳定RNA数据进行实验测试，对算法SA2SP和PBLP的进化代价和时间进行对比。实验结果表明，SA2SP算法能够获得较PBLP算法更小的进化代价，且其运行时间在实际应用中是切实可行的。目前，SA2SP算法仅解决了包含复制和丢失两种操作的演化模型下的两物种小系统发育问题，由于替换（substitution）操作也是一种常见的基因内容修改操作[1]，在今后的工作中将进一步对扩充的问题模型开展深入研究，提高求解方法的实用性。此外，其他智能优化算法对该问题的求解思路也是今后进一步的研究方向。

参考文献：

[1]HOLLOWAY P， SWENSON K， ARDELL D， et al. Ancestral genome organization：an alignment approach [J]. Journal of Computational Biology， 2013，20（4）：280-295.

[2]ANDREOTTI S， KNUT R， CANZAR S. The duplicationloss small phylogeny problem： from cherries to trees [J]. Journal of Computational Biology， 2013，20（9）：643-659.

[3]DONG H， NILSSON L， KURLAND C G. Covariation of tRNA abundance and codon usage in escherichia coli at different growth rates [J]. Journal of Molecular Biology， 1996， 260（5）： 649-663.

[4]KANAYA S， YAMADA Y， KUDO Y， et al. Studies of codon usage and tRNA genes of 18 unicellular organisms and quantification of Bacillus subtilis tRNAs： gene expression level and speciesspecific diversity of codon usage based on multivariate analysis [J]. Gene， 1999，238（1）：143-155.

[5]KOVC J， BREJOV B， VINAR T. A new approach to the small phylogeny problem [EB/OL]. [20150225]. http：///abs/1012.0935.

[6]GUPTA A， MANUCH J， STACHO L， et al. Small phylogeny problem： character evolution trees [M]// SAHINALP S C， MUTHUKRISHNAN S， DOGRUSOZ U. Combinatorial Pattern Matching， LNCS 3109. Berlin： Springer， 2004：230-243.

[7]BLOMME T， VANDEPOELE K， DE BODT S， et al. The gain and loss of genes during 600 million years of vertebrate evolution [J]. Genome Biology， 2006，7（5）： R43.BLOMME T， van DEPOELE K， DE BODT S， et al. The gain and loss of genes during 600 million years of vertebrate evolution [EB/OL]. [20150206]. http：//pmcc.webt.cisti.nrc.ca/picrender.cgi？accid=PMC1779523&blobtype= pdf.

[8]COTTON J A， PAGE R D. Rates and patterns of gene duplication and loss in the human genome [C]// Proceedings of the 2005 Royal Society B： Biological Sciences. London： The Royal Society， 2005，272：277-283.COTTON J A， PAGE R D. Rates and patterns of gene duplication and loss in the human genome [J]. Proceedings of the Royal Society， 2005，272（1560）：277-283.

[9]WAPINSKI I， PFEFFER A， FRIEDMAN N， et al. Natural history and evolutionary principles of gene duplication in fungi [J]. Nature， 2007，449（7158）：54-61.

[10]DEMUTH J P， DE BIE T， STAJICH J E， et al. The evolution of mammalian gene families [J]. PloS One， 2006，1（1）： E85.https：///publication/6616169_The_evolution_of_mammalian_gene_families_PLoS_ONE_11_e85

[11]KEITH J M， ADAMS P， BRYANT D， et al. A simulated annealing algorithm for finding consensus sequences [J]. Bioinformatics， 2002，18（11）：1494-1499.

[12]ZOLA J， TRYSTRAM D， TCHERNYKH A， et al. Parallel multiple sequence alignment with local phylogeny search by simulated annealing [C]// Proceedings of the 2006 20th International Conference on Parallel and Distributed ProcessingInternational Parallel and Distributed Processing Symposium. Washington， DC： IEEE Computer Society， 2006：253.

[13]CHEN S M， LIN C H. Multiple DNA sequence alignment based on genetic simulated annealing techniques [J]. Information and Management SciencesInternational Journal of Information and Management Sciences， 2007，18（2）：97-111.

[14]SARIYER O S， GVEN C. Sequence alignment using simulated annealing [J]. Physica A： Statistical Mechanics and its Applications， 2010，389（15）：3007-3012.

[15]DONDI R， ELMABROUK N. On the complexity of minimum labeling alignment of two genomes [EB/OL]. （20120608） [20150522]. http：///abs/1206.1877v1.