合成生物学最早由Hobom于1980年提出，是建造人工生物系统的一门学科，通过人工操作生物基因分子、生物基因调控与传导等方式实现人工细胞的设计，让细胞帮助人类完成预先设想的各类任务。相比于传统的精细化工，合成生物学更多以可再生的生物质作为原料与能源，更绿色、更环保，符合未来对绿色化工、可持续化工的要求^[1]。目前常见的利用合成生物学人工设计实现合成的高附加值化合物包括抗疟疾药物青蒿素^[2]、抗癌药物紫杉醇^[3]、镇咳药物那可丁^[4]以及抗肿瘤药物人参二醇^[5]等。

代谢途径重构是合成生物学工业菌种开发流程中的重要步骤，旨在初步确定宿主内合成目标物质的转化路径，是以酶为代表的蛋白质定向进化与后期代谢流优化的基础。基于人工智能算法的生物逆合成分析是实现代谢途径重构的重要手段，其目标是实现宿主内合成路径的获取与推荐，减少实验与尝试的次数，进而降低生产目标化合物的成本。对于生物逆合成分析的路径结果，其规模会随着路径长度的增加而迅速增加，相应的路径数量也随之增加，进而需要对路径进行评估与排序，选出指定数量的最可能在宿主内实现的目标产物合成路径。路径评估的准确性决定了生物逆合成分析方法与工具的可行性与实际价值。

生物逆合成分析包括3个主要步骤：路径搜索、路径发现和路径评估。路径搜索是基于反应数据库进行的反应路径延伸。首先在反应数据库中找到反应网络中目标产物对应的化合物节点，然后逆有向边方向进行图搜索，就实现了路径延伸。路径搜索的质量取决于反应数据库的规模以及反应信息的正确性与完整性。目前，包含路径搜索的生物逆合成分析工具一般只使用单一的反应数据库，反应网络的规模与信息完整性受到限制，故ATLAS^[6]、RetroPath^[7]等工具都考虑引入路径发现方法来扩展反应网络的规模、补充缺失的反应信息。

路径发现是基于反应规则产生的虚拟反应集进行的反应路径延伸，基础是酶的底物杂泛性。通过运行反应规则可以产生一系列的虚拟反应。反应规则来自真实反应，代表了酶促反应模式，目前仍缺乏高效、自动提取反应规则的方法。路径发现的质量取决于反应规则的数量、准确性以及泛化能力。ATLAS^[6]、RetroPath^[7]、Transform-Miner^[8]等生物逆合成分析工具均引入反应规则与路径发现进行路径延伸。

路径搜索与路径发现得到的生物逆合成分析网络，一般会包含大量的预测路径，并非所有预测路径都能实现目标催化功能，不能实践的路径称为假阳性路径^[1]。路径评估即是对包含有假阳性路径的大量预测路径进行评价与排序，进而推荐最优候选路径。已有的生物逆合成分析工具的路径评估考虑了包括理论可行性和宿主适配性两方面的指标^[1]。下面以一个理论可行性指标和两个宿主适配性指标为例进行简要介绍。

1) 热力学可行性。符合热力学规律是理论可行性的基础。运行反应规则产生的虚拟反应并未被实际确认，因此需要判断反应是否具有热力学的可行性。ATLAS^[6]、XTMS^[9]在路径评估时考虑了反应的热力学可行性，使用基团贡献法计算单步反应底物和产物的热力学性质，进而计算细胞温度下反应的Gibbs自由能变化，确定细胞环境下反应的热力学可行性。

2) 化合物毒性。考虑到生物逆合成分析的结果需在宿主(细胞)中重现，一般要求路径中不出现高毒性的化合物，保证路径的宿主适配性。半数抑制浓度(IC50)可用来评价化合物的毒性。Planson等^[10]提出了一种通过训练模型基于物质结构预测化合物毒性的方法，并将预测的化合物毒性作为评估指标。

3) 理论产量。生物逆合成分析中希望能使目标产物的产量最大化，即最大化目标产物在宿主中沿某条路径的理论产量。理论产量是路径的宿主适配性的一种定量体现。GEM-Path^[11]和XTMS^[9]均考虑以最大化理论产量作为优化目标进行路径评估，利用流量平衡分析(flux balance analysis, FBA)计算理论产量。但这种方法需要大量的反应以及反应环境的相关信息，故只适用于一些发生于那些得到充分研究的模式微生物内的常见反应，目前使用并不广泛。

在实际的生物逆合成分析中，一般同时考虑单步反应指标以及整体的路径指标，将这些指标通过合理的加权组合进行路径的综合评价，使其符合实际的背景反应数据库。只有评估方法与背景反应数据库相匹配，才能提高合成路径的准确与可靠性，保证最终路径数量与质量兼得。

各路径评估指标并不是彼此独立的。例如，路径长度与理论产量就高度相关，且各个针对单步反应的指标也需要进行一定的运算得到该指标对应的整体路径数值，这也与路径长度高度相关。考虑到这种相关性，合理的路径评估指标的权重选取是相对困难的。同时，不同指标可能的数值单位与分布不同，在具体的组合前首先要对各指标的数值结果进行区间统一，保证权重可以反映各指标对最终路径评估分数的影响。

总的来说，目前已有的生物逆合成分析相关工具，例如ATLAS^[6]、Transform-Miner^[8]等，可以基于逆合成分析方法进行代谢途径重构，但这些工具对路径的评估方法均较为简单，只是给出少量的指标数值，且不区分指标的权重与优先级，无法给出路径的排序与推荐。本文针对此问题，构建了一种基于理论可行性的生物合成路径评估与排序的方法，使用6个路径评估指标进行路径推荐。

1 路径评估方法

生物逆合成分析中获取的反应路径长度彼此不同，含有不同的特征，需要进行路径特征的提取与量化。现有工具^[6-8]一般仅考虑几个单一指标，指标的处理与加权组合方法也相对简单。本文的路径评估方法从6个指标维度对每条路径进行评分，并基于层次分析法，从指标重要性比较的角度，获得6个指标对应的权重，计算指标的加权和作为路径评估的结果，再基于结果对各条路径进行排序，给出最终合成路径推荐。考虑到宿主适配性指标计算困难，需要较多的数据基础，故本文的路径评估方法仅选择了理论可行性指标。

本文构建的6个路径评估指标及其基本含义如表 1所示。

1.1 路径长度分数

路径长度代表每条独立路径含有的单步反应数量，即完成指定合成转化需要的反应个数。由于每步反应都需要酶进行催化，更长的路径一方面意味着引入更多的酶，增加了成本与宿主的代谢负担；另一方面也意味着原始底物需要被多次按顺序修饰，增加了出现副产物的可能性，导致更低的选择性与理论产量，同样增加了合成成本。故在生物逆合成分析中倾向于更短的独立路径。

基于此，路径长度分数len的计算方法为

$ \text { len }=\frac{1}{n} . $

(1)

其中n代表该独立路径长度。路径越长，该路径的路径长度分数越低。在评估时倾向于更高的路径长度分数。

在描述之后的5个指标的计算方法前，为了描述方便，基于每条独立路径中的每个单步反应的来源以及真实单步反应是否连接到真实酶信息，将所有单步反应分为表 2所示的5类反应。表 2中对5类单步反应依次进行编号，并给出了其类别特性与数量划分(n₁，n₂，…，n₅)。

基于分类的完备性，应有

$ n=n_1+n_2+n_3+n_4+n_5 . $

(2)

1.2 路径反应数据库比例

为了更好地在宿主内重现路径，路径评估中更倾向于来自反应数据库内的真实单步反应。基于表 2对独立路径中单步反应的分类，路径中最终保留反应数据库内真实反应比例Rh_frac的计算方法为

$ \text { Rh_frac }=\left(n_1+n_2+n_3+n_4\right) / n. $

(3)

路径反应数据库比例是一个介于0~1之间的数。在评估时倾向于更高的路径反应数据库比例。

1.3 路径相似性分数

对于生化反应，碳链一般作为保守基团不参与反应，一定程度上保证了单步反应的底物与产物的结构相似性。基于此，考虑到路径的可行性，在路径评估中更倾向于单步相似性高的反应。对于每条独立路径，路径相似性分数path_sim的计算方法为

$ \text { path_sim }=\prod\limits_{i=1}^n \text { rxn_sim }{ }_i \text {. } $

(4)

对路径中每个反应i计算底物与产物的分子指纹结构相似性rxn_sim_i，将n个单步结构相似性的乘积作为路径相似性分数。在评估时倾向于更高的路径相似性分数。

1.4 路径可能性分数

虚拟反应可行性建立在酶的底物杂泛性上，不同分子被同一个酶催化发生反应的概率不同。使用提取反应规则的真实反应底物与应用反应规则的底物的结构相似性comp_sim衡量单步虚拟反应的发生概率。假设化合物A与B之间的单步反应为真实反应，通过提取反应规则并将C作为反应底物，得到了化合物C转化为化合物D的虚拟反应。计算化合物A与C的结构相似性可以衡量该单步虚拟反应的发生可能性。结构相似性越高，虚拟反应发生的可能性越大。对整条路径，路径可能性分数path_prob的计算方法为

$ \text { path_prob }=\prod\limits_{i=1}^{n_5} \mathrm{e}^{\text {comp_sim }_i-1} \text {. } $

(5)

其中n₅代表独立路径中第5类单步反应的数量，即完全来自路径发现的虚拟反应的数量。在评估时倾向于更高的路径可能性分数。

1.5 路径酶比例

酶是生化反应可以在宿主内发生且高效发生的基础。为了保证路径可在宿主内重现，减少酶序列缺失导致的效率降低，在生物逆合成分析中一般要求路径内的每步反应均有可靠的酶信息。路径评估中更倾向于缺失酶信息少的路径。路径酶比例enz_frac的计算方法为

$ \text { enz_frac }=1-\frac{n_1}{n} . $

(6)

其中n₁代表路径中第1类单步反应的数量，这类反应也是唯一一类最终路径结果中依然缺乏酶信息的单步反应。在评估时倾向于更高的路径酶比例。

1.6 路径规则分数

路径发现中基于真实反应获得了一系列反应规则，不同反应规则产生虚拟反应的可行性与实现难度不同。使用路径规则分数对反应规则进行评价。每条反应规则的评估分为两部分，奖励(reward)和惩罚(penalty)。两者的计算方法如下：

$ \text { reward }=\frac{\lg n_{\text {rule }}}{\sqrt{R}}. $

(7)

$\text { penalty }=\frac{\lg C_{\mathrm{rule}}}{\sqrt{R}}. $

(8)

其中：R代表规则半径，n_rule代表该反应规则在相同酶催化下可由多少个反应数据库内的真实反应提取得到。C_rule代表与该规则关联的酶序列的类别数量，统计类别时仅考虑酶EC编号的前3位。对于奖励，n_rule越大，意味着反应规则在实际反应内出现的概率越高，反应规则的可信度越高，得到的虚拟反应的可行性越好。对于惩罚，C_rule越大，意味着反应规则提取酶促反应模式专一性差，无法体现相应酶促反应的本质，虚拟反应的实现难度大。

基于单个反应规则的奖励与惩罚，路径规则分数path_rule的计算方法为

$ \text { path_rule }=\prod\limits_{i=1}^{n_4+n_5} \mathrm{e}^{\text {reward }_i-\text { penalty }_i-1}. $

(9)

其中：n₄、n₅代表路径中第4、第5类单步反应的数量，这两类单步反应利用了反应规则的酶信息作为反应酶信息。将单步反应规则分数的乘积作为路径规则分数。在评估时倾向于更高的路径规则分数。

2 路径评估实例

本节基于生物逆合成分析获取一对指定物质之间的转化路径集合，将该集合作为路径对象来证明本文提出的路径评估方法的有效与合理性。目前还未将路径推荐结果进行实际验证，合理性判断取决于实践中对转化路径的偏好。

2.1 路径获取

基于生物逆合成分析方法，集成路径搜索与路径发现实现双向路径获取。真实反应来自Rhea数据库^[12]，虚拟反应来自运行作者基于真实反应人工提取的107条SMARTS(SMiles Arbitrary Target Specification)反应规则。考虑到本文的重点为路径评估方法，故不在本文中仔细介绍路径获取的过程与方法。路径获取中选择的化合物数据库为ChEBI (Chemical Entities of Biological Interest)数据库^[13]，用化合物在ChEBI数据库内的编号指代其本身。

选择苯甲酸根(ChEBI编号为16150)与2, 3-二羟基苯甲酸根(ChEBI编号为36654)为路径的原始底物与目标产物，获取两物质之间的转化路径。选择最大路径长度为9、单步底物结构相似性阈值为0.5，得到图 1的路径结果，图中的实线代表反应数据库内的真实反应，图中的虚线代表基于反应规则产生的单步虚拟反应。通过深度优先搜索，图 1共获得了12条路径，最短的路径长度为2，最长的路径长度为9。这12条路径都是潜在的可以完成原始底物到目标产物转化的途径，不同的路径经历不同的反应步数与反应类型。

图 1中出现的中间产物均不需要人工添加，只作为路径的一个节点。为了限制路径规模，设置了单步底物结构相似性阈值。每个单步反应的结构相似性在有向边上标出。图中红色节点代表路径的原始底物与目标产物，每个节点上的数字为该物质的ChEBI编号。

2.2 路径评估 2.2.1 指标计算

基于2.1节的路径获取方法与第1节的路径评估方法，对12条路径进行6个路径评估指标的计算，将结果分别列入表 3第2—7列。指标序号与表 1保持一致。

2.2.2 权重计算

使用层次分析法确定各指标权重。图 2中展示了3位专家基于1—5标度法对各指标的相对重要程度进行打分得到的判断矩阵，矩阵的二维指标序号与表 1保持一致。基于每个判断矩阵分别进行一致性检验与权重确定，结果如表 4所示。其中：λ_max代表各判断矩阵的最大特征值；w₀代表最大特征值对应的归一化特征向量；I_C代表一致性指标，衡量判断矩阵的一致性程度。

由表 4可知，所有的判断矩阵均通过了一致性检验，均可被接受。将不同判断矩阵得到的w₀在考虑I_C的情况下进行融合得到复合权重，

$ \boldsymbol{w}=\frac{\sum\limits_{i=1}^3 \boldsymbol{w}_{0, i} X\left(1-I_{\mathrm{C}, i}\right)}{\sum\limits_{i=1}^3\left(1-I_{\mathrm{C}, i}\right).} $

(10)

最终得到的复合权重w=[0.082 0.168 0.168 0.242 0.125 0.216]。

基于复合权重以及3位专家各自的权重计算各条路径评估指标的加权和，将计算结果列入表 3第8—11列。

2.2.3 路径评估结果

基于路径评估分数的计算结果(表 3)，最优的3条独立路径编号分别为4、2和3，将这3条路径展示在表 5中。表 5展示了路径涉及的单步反应产物(ChEBI编号)、酶(EC编号)以及反应类别。路径编号与表 3保持一致，反应类别与表 2保持一致。独立路径4的优势在于路径较短、全部由真实反应构成、不缺失酶信息，是最理想的路径。路径2的优势在于路径最短，同时虚拟反应帮助真实反应补充了酶信息。路径3的优势在于路径全部由真实反应构成，且结构相似性高。这3条路径在某几个维度上有明显优势，证明了本文提出的路径评估结果的合理性。特别需要说明的是，由于不同专家对各指标具有不同的主观偏好，基于每位专家权重单独计算可能会得到不同的最优路径排序结果。

3 结论

针对目前已有路径评估方法使用指标少、不区分指标之间的优先级的问题，本文基于理论可行性构建了一种生物合成路径评估方法。在指标数量方面，设计了6个反映路径理论可行性的路径评估指标，包括路径长度分数、路径反应数据库比例、路径相似性分数、路径可能性分数、路径酶比例和路径规则分数，考虑了路径长度、反应可行性与反应规则评价。在指标优先级方面，为了科学地实现指标的加权求和，邀请了3位专家对指标的相对重要程度进行两两打分。基于层次分析法，结合一致性指标，最终确定了6个指标的加权权重，得到了具有合理性的路径评估结果。苯甲酸根与2, 3-二羟基苯甲酸根之间生物合成路径评估的模拟实例测试，证明了该路径评估方法可以合理、可解释地推荐出路径结果，建立了路径获取与实验尝试的连接。本文的路径评估方法达到了预期结果，可用于应对合成生物学中数据丰富性与实践难度高的矛盾，从而降低实验的试错成本，为增加生物逆合成工具的实用性提供了基础。