融合原子交换特征信息的代谢路径预测

黄毅然，万志远，钟 诚

(1.广西大学 计算机与电子信息学院，广西 南宁 530004；2.广西大学 广西高校并行与分布式计算技术重点实验室，广西 南宁 530004；3.广西大学 广西多媒体通信与网络技术重点实验室，广西 南宁 530004)

0 引 言

代谢路径是生物体通过一组代谢反应将起始代谢物转化为目标代谢物的路径。代谢路径预测是生物合成路径的设计-构建-测试循环中的重要步骤之一，主要目标是通过计算方法在代谢数据库中寻找到新的有价值的合成路径[1-3]。

在代谢数据库中找到的合成代谢路径通常包含连接度很高的簇代谢物[4，5]。簇代谢物的出现可能会产生生物学上没有意义的代谢路径。原子交换是代谢路径中的代谢反应中的底物和产物之间存在原子转移的现象，研究结果表明通过利用原子交换信息搜索路径可以有效排除合成路径中的簇代谢物[6-8]。

一些基于约束的代谢路径预测方法利用代谢网络中代谢物和代谢反应具有的生化特性，构建化学计量MILP模型以寻找代谢路径，并通过在MILP模型中利用原子交换信息来寻找没有簇代谢物的替代路径以提高路径搜索质量[9]。NetFlow方法通过检测到的碳通量搜索具有生化相关性的代谢路径[10]。NFP方法[11，12]结合代谢网络中的拓扑结构和化学计量信息以构建特定的MILP模型，通过MILP模型实现从任意起始代谢物到目标代谢物的代谢路径搜索。虽然NFP方法可以搜索从任意起始代谢物到目标代谢物的代谢路径，但它没有将代谢反应中的原子转移信息融入MILP模型中以寻找生化可行的代谢路径，使得该方法获得的代谢路径的生化相关性较弱。

本文研究提出一种基于约束的代谢路径预测方法PVA，以搜索出从任意起始代谢物到目标代谢物的反应中包含特定原子交换的代谢路径。不同于搜索给定的起始代谢物和目标代谢物之间代谢路径的路径预测方法[6-8]，本文提出的方法可以搜索从任意起始代谢物到给定目标代谢物的代谢路径。此外，与NFP方法不同，本文提出的方法将代谢反应中特定原子交换信息与反应化学计量相结合，建立基于约束的MILP模型以找到含有特定原子交换的代谢路径，以进一步提高从任意起始代谢物开始搜索路径所找到的路径的生化相关性，更好地发现与已知代谢路径相近的替代路径。

1 方 法

本文方法的主要思想是：首先，从KEGG数据库中提取在输入代谢物和输出代谢物之间存在特定原子交换的代谢反应。然后，将反应的化学计量信息与包含特定原子交换信息的代谢反应结合起来，构建基于约束的MILP模型。最后，通过求解MILP模型来找到从任意起始代谢物到目标代谢物的反应中包含特定原子交换的代谢路径。

1.1 代谢路径预测模型

首先利用KEGG RPAIR数据库中输入代谢物和输出代谢物之间存在特定原子交换的代谢反应来重构代谢网络。在重构的代谢网络中，顶点代表代谢物，弧表示代谢物之间的代谢反应，弧中的反应包含输入代谢物和输出代谢物之间的原子交换特征信息。

重构的代谢网络由一组反应R和一组代谢物M组成。本文利用化学计量矩阵S、二元系数G、一组内部代谢物I、一组外部代谢物E、一组起始代谢物SM和一组可逆反应RB来刻画重构的代谢网络。在化学计量矩阵S中，行代表代谢物，列代表代谢反应。Smr表示反应r中代谢物m的化学计量值。Smr<0表示m作为反应r的输入代谢物，Smr>0表示m作为反应r的输出代谢物。内部代谢物集合I给出了代谢网络中的内部代谢物列表。外部代谢物集合E给出了需要通过预测路径中的代谢反应产生的外部代谢物列表。起始代谢物集合SM给出了可作为预测路径中的第一个代谢物的代谢物列表。可逆反应集RB给出了重构的代谢网络中的可逆反应。

依据上述模型信息，本文使用如下代谢路径预测的约束方程[11]

(1)

(2)

(3)

(4)

(5)

在上述约束方程中，二元变量aij表示输入代谢物i和输出代谢物j之间的弧，其中|A|表示代谢网络中有|A|个不同的代谢物。如果输入代谢物i到输出代谢物j的弧aij存在于预测路径中，则aij=1，否则为0，i，j=1，2，…，|A|。 约束方程(1)选择目标代谢物p为预测路径的终点[11]。约束方程(2)确保目标代谢物p不会是预测路径的起点[11]。约束方程(3)确保起始代谢物集合SM中的代谢物l或者是中间代谢物，或者是起始代谢物[11]。约束方程(4)确保起始代谢物集合以外的其它代谢物都可以作为代谢路径的中间代谢物[11]。约束方程(5)确保预测路径中的代谢物都是唯一的[11]。约束方程(1)～方程(5)确保代谢路径预测方法可以从任意起始代谢物开始搜索，其中起始代谢物可以通过约束方程(3)从起始代谢物集合SM中随机选择。

通过约束方程(1)～方程(5)只是获得代谢路径中简单的代谢物序列。下面将通过约束方程(6)～方程(10)定义代谢路径中的反应通量

(6)

(7)

cr≤vr≤Max*crr=1，…，|N|

(8)

ca+cb≤1 ∀(a，b)∈RB

(9)

(10)

在约束方程(6)～方程(10)中，vr是反应r的通量，它为每个反应r分配一个非负通量，|N|表示代谢网络中有|N|个的不同代谢反应，r=1，2，…，|N|。 二元变量cr用作控制反应通量vr的开闭，cr=0使得相关通量vr为0，而cr=1使vr为非零值，r=1，2，…，|N|。Max大于或等于1。

约束方程(6)～方程(10)定义了预测路径的有效通量分布以确保找到的预测路径在生物学上是可行的，即该预测路径可以经过多个代谢反应由起始代谢物合成目标代谢物。约束方程(6)确保外部代谢物只能产生，而不能消耗[11]。约束方程(7)确保目标代谢物p作为产物存在于预测路径中[11]。约束方程(8)通过二元变量cr控制反应通量vr的切换[11]。约束方程(9)约束代谢网络中的可逆反应不能同时出现在预测路径中[11]。

本文结合约束方程(1)～方程(10)来求解MILP模型中的目标函数式(11)，以获得从任意起始代谢物到目标代谢物p的反应中含有特定原子交换的预测路径

(11)

1.2 算 法

对于给定目标代谢物T以及被追踪的原子AM，使用融合原子交换特征信息的代谢路径预测算法PVA，搜索从任意起始代谢物到目标代谢物T的反应中含有原子AM交换的代谢路径的步骤如下：

步骤1 从KEGG RPAIR数据库中提取输入代谢物和输出代谢物之间存在原子AM交换的反应，利用存在原子AM交换的代谢反应以建立MILP模型。

步骤2 对于目标代谢物T，通过MILP模型中的约束方程(1)～方程(5)找到从任意起始代谢物到目标代谢物T的代谢路径的代谢物序列。通过MILP模型中的约束方程(6)～方程(9)获取从任意起始代谢物到目标代谢物T的代谢路径的代谢反应通量。通过MILP模型中的约束方程(10)从代谢路径中筛选出存在原子AM交换的代谢反应。

步骤3 结合MILP模型约束方程(1)～方程(10)求解方程(11)的MILP模型中的目标函数，以获得从任意起始代谢物到目标代谢物T的反应中含原子AM交换的代谢路径。

算法1形式描述了本文提出的算法PVA。

算法1：PVA

输入：目标代谢物T，追踪的目标原子AM，KEGG RPAIR数据库的代谢反应集合reactions

输出：任意起始代谢物到目标代谢物T的反应中含有原子AM交换的代谢路径

Begin

(1) for allrin reactions do

RAM← 遍历reactions中每个反应，获取代谢反应r中具有原子AM交换特征信息的反应；

end for

(2)POM←利用MILP模型中的约束方程(1)～(5)找到任意起始代谢物到目标代谢物T之间的代谢物序列；

(3)RS← 利用MILP模型中的约束方程(6)～(9)找到与代谢物序列POM相匹配的代谢反应通量；

(4) for allrinRSdo

POR←通过MILP模型中的约束方程(10)从RS中筛选出代谢反应r中存在原子AM交换的反应；

end for

(5)path←结合任意起始代谢物到目标代谢物T之间的代谢物序列POM和存在原子AM交换的反应集合POR，求解MILP模型中目标函数(方程(11))以获得从任意起始代谢物到目标代谢物T的反应中含有原子AM交换的代谢路径；

(6) 输出预测路径结果path；

End.

2 实 验

2.1 实验环境

通过比较预测路径与来自KEGG的大肠杆菌数据库中的20条已知路径来评估代谢路径预测算法的有效性。

在实验中，本文方法PVA采用4种代谢路径预测模式来寻找从任意起始代谢物到目标代谢物的反应中含有特定原子交换的代谢路径。PVA-C表示通过追踪反应中的碳原子交换来搜索从任意起始代谢物到目标代谢物的代谢路径。PVA-O表示通过追踪反应中的氧原子交换来搜索从任意起始代谢物到目标代谢物的代谢路径。PVA-N表示通过追踪反应中的氮原子交换来搜索从任意起始代谢物到目标代谢物的代谢路径。PVA-P表示通过追踪反应中的磷原子交换来搜索从任意起始代谢物到目标代谢物的代谢路径。

通过MILP模型搜索代谢路径的预测方法有两种。第一种方法(表示为“Topology”)在KEGG的原代谢网络中使用MILP模型搜索从任意起始代谢物到目标代谢物的代谢路径。第二种方法(表示为“Hubs”)在调整后的代谢网络中使用MILP模型搜索从任意起始代谢物到目标代谢物的代谢路径，调整后的代谢网络是指删除了高度连接的代谢物(Hubs)的代谢网络[13]。

本文分别利用方法PVA-C、PVA-O、PVA-N、PVA-P、NFP[11]、Topology和Hubs[13]在KEGG数据库中搜索20条已知代谢路径。每种方法搜索从任意起始代谢物到20条已知代谢路径中的目标代谢物的代谢路径，并通过将计算得到的替代路径与已知代谢路径进行比较来评估代谢路径预测方法的有效性。

本文采用Java和GUROBI编程实现PVA算法。基于约束的代谢路径预测方法NFP和其它两种代谢路径预测方法(“Topology”和“Hubs”)的MILP模型都使用GUROBI在Java环境中求解。所有实验均在 Intel Xeon 6130和192 GB RAM的计算机上运行。

2.2 性能指标

代谢路径预测方法得到的预测路径的生化相关性会影响预测路径在代谢工程中生物合成应用的可行性。本文通过将预测得到的替代路径与已知路径进行比较以评估预测路径的生化相关性。与同类研究方法相同，本文以代谢反应准确率[7]和代谢物准确率[7]作为评估预测路径的性能指标，具体计算方法如下：

(1)代谢反应准确率R_ACC=(R_Sn+R_PPV)/2， 其中R_S是代谢反应的灵敏度，R_Sn=R_TP/(R_TP+R_FP)，R_PPV是代谢反应的阳性预测值，R_PPV=R_TP/(R_TP+R_FN)[7]。 其中真阳性R_TP表示预测路径中存在的代谢反应同已知路径中的代谢反应相同，且在预测路径和已知路径中的代谢反应出现的相对顺序是相同的；假阳性R_FP表示存在于预测路径中，但不存在于已知路径中的反应；假阴性R_FN表示不存在于预测路径中，而是存在于已知路径中的反应[7]。代谢反应准确率R_ACC值越高说明代谢路径预测算法返回的预测路径中的代谢反应同已知路径中的代谢反应越趋于一致[7]。

(2)代谢物准确率C_ACC=(C_Sn+C_PPV)/2， 其中C_Sn是代谢物灵敏度，C_Sn=TP/(TP+FP)，C_PPV是代谢物的阳性预测值，C_PPV=TP/(TP+FN)[7]。 真阳性TP表示预测路径中存在的代谢物同已知路径中的代谢物相同，且在预测路径和已知路径中的代谢物出现的相对顺序是相同的；假阳性FP表示存在于预测路径中，但不存在于已知路径中的代谢物；假阴性FN表示代谢物不存在于预测路径中，而是存在于已知路径中[7]。代谢物准确率C_ACC值越高说明代谢路径预测方法能够更准确地发现已知路径中存在的代谢物[7]。

2.3 预测路径与已知路径中代谢反应和代谢物的比较

对于20条已知路径中的每个目标代谢物，分别采用预测方法PVA-C、PVA-O、PVA-N、PVA-P、NFP、Topo-logy和Hubs搜索从任意起始代谢物到目标代谢物的代谢路径。通过将每种方法找到的替代路径与已知路径进行比较，以评估代谢路径预测方法的性能。

表1给出了各算法预测得到的替代路径中的首条路径、前5条路径、前10条路径以及前100路径的平均代谢反应准确率R_ACC值。

表1 预测路径的平均代谢反应准确率R_ACC值

从表1可以看出，对于首条路径、前5条路径、前10条路径以及前100路径，PVA-O获得的预测路径的R_ACC值都比NFP、Topology、Hubs、PVA-C、PVA-N和PVA-P 获得的预测路径的R_ACC值更高；此外PVA-C和NFP获得的预测路径的R_ACC值与PVA-O获得的预测路径的R_ACC值相近。这些结果表明，相比其它基于约束的代谢路径预测方法，PVA-O能够更准确地恢复已知路径中的代谢反应。

从表1中还可以看到，对于首条路径、前5条路径、前10条路径以及前100路径，PVA-P获得的预测路径的R_ACC值低于PVA-C、PVA-O、PVA-N和NFP获得的预测路径的R_ACC值；PVA-P获得的预测路径的R_ACC值与Topology和Hubs获得的预测路径的R_ACC值相当。这些结果表明，在基于约束的MILP模型中追踪不同的原子交换对发现已知路径中的代谢反应有较大的影响。

表2给出了各算法预测得到的替代路径中的首条路径、前5条路径、前10条路径以和前100路径的平均代谢物准确率C_ACC值。

表2 预测路径的平均代谢物准确率C_ACC值

从表2中可以观察到，对于首条路径、前5条路径、前10条路径以及前100路径，PVA-O获得的预测路径的C_ACC值也都高于NFP、Topology、Hubs、PVA-C、PVA-N和PVA-P获得的预测路径的C_ACC值。这些结果表明PVA-O能够更准确地找到已知路径中的代谢物。另一方面，从表2还可以看到，对于首条路径、前5条路径、前10条路径以及前100路径，与PVA的其它追踪模式相比，PVA-N和PVA-P获得的预测路径的C_ACC值比较接近，但是都低于PVA-C、PVA-O获得的预测路径的C_ACC值。这些结果表明，在基于约束的MILP模型中通过追踪磷原子和氮原子可能较难获得生化相关较好的代谢路径。

表1和表2的结果表明，融合原子交换特征信息的代谢路径预测算法PVA能够有效地找到与已知路径更接近的预测路径，其中PVA-O的预测路径的反应准确率和代谢物准确率都较高，说明已知路径反应中的氧原子交换特征更显著，人们可以通过追踪氧原子在代谢路径反应中的交换更容易地发现生化相关性较好的替代路径。

2.4 实例分析：UDP-葡萄糖代谢路径的搜索

第2.3节的实验结果表明融合原子交换特征信息的代谢路径预测算法PVA可以有效发现生化相关性较好的代谢路径。下面将通过一个实例来研究PVA所找到的路径的生化特点。

UDP-葡萄糖(UDP-glucose)是一种核苷酸糖，在许多糖基化反应中充当葡萄糖残基的供体，它可以通过以α-D-葡萄糖(alpha-D-Glucose)为起始代谢物经过多个中间代谢物由多个代谢反应合成获得[14，15]。同时，UDP-葡萄糖在一系列代谢物的糖基化中必不可少[15]。

图1(a)是KEGG数据库中从起始代谢物alpha-D-Glucose到目标代谢产物UDP-glucose的已知代谢路径，该已知代谢路径存在于大肠杆菌和酵母菌中。图1(b)是PVA-C搜索从任意起始代谢物到目标代谢产物UDP-glucose所获得的前5条预测路径，图1(c)是PVA-O搜索从任意起始代谢物到目标代谢产物UDP-glucose所获得的前5条预测路径，图1(d)是PVA-P搜索从任意起始代谢物到目标代谢产物UDP-glucose所获得的前5条预测路径。在图1中，预测路径中与已知路径相同的代谢反应和代谢物分别用虚线矩形和虚线椭圆绘制。

图1 UDP-glucose合成路径

从图1中可以发现PVA-C、PVA-O和PVA-P选取了6种不同的代谢物作为起始代谢物搜索代谢路径，其中就包含起始代谢物alpha-D-Glucose。

从图1(b)和图1(c)可以看到，PVA-C和PVA-O分别发现了从起始代谢物alpha-D-Glucose到目标代谢物UDP-glucose的预测路径b4和d4，这两条代谢路径都是大肠杆菌和酵母菌中常用于合成目标代谢物UDP-glucose的代谢路径[14]。对比图1(b)和图1(c)还可看出，PVA-C和追踪氧原子交换的方法获得的前5条预测路径是相同，只是返回的路径排序不同。这些结果表明通过追踪碳原子和氧原子在代谢路径反应中的交换能较好地寻找到与合成代谢物UDP-glucose的已知路径相似的替代路径，在一定程度上可以帮助代谢工程研究人员更方便地设计合成目标代谢物UDP-glucose的替代路径。

从图1(b)、图1(c)和图1(d)还可以看到，虽然PVA-P没有返回已知路径，但是PVA-P找到的从起始代谢物alpha-D-Glucose6-phosphate到目标代谢物UDP-glucose的路径f1有很大一部分与已知路径相同。这些结果表明，通过追踪磷原子在代谢路径反应中的交换也可以为研究人员发现合成目标代谢物UDP-glucose提供潜在有用的替代路径。

3 结束语

本文提出的代谢路径预测算法的特色是：从KEGG代谢数据库中获取含有原子交换信息的代谢反应和代谢物，并结合反应化学计量信息构建基于约束的MILP模型，通过求解该MILP模型以获取不限定起始代谢物的生化相关性较好的代谢路径。实验结果表明，通过在基于约束的代谢路径预测模型中追踪原子交换，能够获得从任意起始代谢物到目标代谢物的反应中含有原子交换的生化相关性更好的替代路径，为合成路径的设计提供启发和参考。