来源：Evaluation of machine learning-assisted directed evolution across diverse combinatorial landscapes（Cell Systems）| 编译：DNA Lab Space | 原文许可：elsevier-subscription
导读
定向进化是蛋白质工程的核心手段，但耗时费力。机器学习辅助定向进化（MLDE）虽已被证明更高效，其适用条件却始终模糊。本研究系统评估了MLDE、主动学习及聚焦训练等策略在16个蛋白质适应度景观上的表现，通过量化景观可导航性，揭示了MLDE在越崎岖的景观上优势越明显。研究发现，结合零样本预测器的聚焦训练能稳定超越随机采样，为湿实验策略选择提供了实用指南。
研究背景
工程化蛋白质在疾病治疗、作物改良和绿色催化等领域不可或缺。定向进化（DE）通过反复突变和功能筛选累积有益突变，本质上是高维适应度景观上的经验性贪心爬坡过程。然而，DE耗时且资源密集——筛选昂贵，且可能需要多轮突变和筛选才能获得理想改进。当适应度景观富含氨基酸替换的上位效应（非加性效应）时，景观变得更加崎岖、难以穿越。上位效应常见于结构邻近的突变之间，且在结合界面或酶活性位点因残基、底物和辅因子间的直接相互作用而富集。蛋白质工程师常针对这些互作位点进行突变，通常使用同时饱和突变（SSM）构建文库。
尽管多种MLDE策略已被证明能比传统DE更高效地鉴定高适应度变体，但人们对影响MLDE跨蛋白质表现的因素了解有限，阻碍了湿实验策略的最优选择。
研究方法
关键资源
本研究使用的数据包括：细菌毒素-抗毒素（ParD-ParE）适应度景观（Lite et al.，https://doi.org/10.7554/eLife.60924）及其结构（PDB: 6X0A; PDB: 5CEG）；Protein G B1结构域（GB1）适应度景观（Wu et al.，https://doi.org/10.7554/eLife.16965）及其结构（PDB: 2GI9）；二氢叶酸还原酶（DHFR）适应度景观（Papkou et al.，https://doi.org/10.1126/science.adh3860）及其结构（PDB: 6XG5）；T7 RNA聚合酶适应度景观（Tu et al.，https://doi.org/10.1101/2022.03.09.483646）及其结构（PDB: 1CEZ）；TEV蛋白酶适应度景观（Tu et al.，https://doi.org/10.1101/2022.03.09.483646）及其结构（PDB: 1LVM）；色氨酸合酶β亚基（TrpB）适应度景观（Johnston et al.，https://doi.org/10.1073/pnas.2400439121）及其结构（PDB: 8VHH）。编译数据与结果见 https://doi.org/10.5281/zenodo.13910505。
使用的软件与算法包括：SSMuLA代码及conda环境（https://github.com/fhalab/SSMuLA; https://doi.org/10.5281/zenodo.15694110）；ALDE4SSMuLA（https://github.com/fhalab/alde4ssmula; https://doi.org/10.5281/zenodo.15694151）；EVmutation网络服务器（Hopf et al.，https://v2.evcouplings.org/）及GitHub仓库（https://github.com/debbiemarkslab/EVmutation）；ESM-2（Rives et al.; Meier et al.; Lin et al.，https://github.com/facebookresearch/esm）；ESM逆折叠（ESM-IF）（Hsu et al.，https://github.com/facebookresearch/esm）；基于结构的组合变体效应（CoVES）（Ding et al.，https://github.com/ddingding/CoVES）；Triad（Protabit, Pasadena, CA, USA，https://triad.protabit.com/）；MLDE（Wittmann et al.，https://github.com/fhalab/MLDE）；微调蛋白质语言模型（Schmirler et al.，https://github.com/RSchmirler/data-repo_plm-finetune-eval）。
景观准备
为避免偏差和误述（尤其对推导景观属性而言），本研究选择了基本完整的数据集，未进行任何插补。假设缺失值遵循与现有数据相同的分布，因此不影响属性计算。为减少噪声数据或可靠性较低的景观属性计算带来的偏差（这些偏差可能导致不可推广的结论），正文聚焦于活性变体至少占1%的景观。然而，为确保比较不同方法的结论在所有景观中均有效，我们在补充信息中提供了对活性变体少于1%的景观的广泛分析，并在适用处讨论了差异。1%阈值的推导基于传统DE筛选中最大景观出现一个活性变体的预期发生率，计算为(1 / (4 × 20)) × 100%。所有变体适应度值在每个景观内归一化，使最大适应度变体赋值为1：ω′ = ω / ω_max，其中ω为变体原始适应度值，ω_max为景观内观测到的最大适应度，ω′为所有分析中使用的归一化适应度值。
景观属性
为全面理解模型结果，本研究分析了两组属性：1）适应度统计，包括活性变体百分比和简单统计建模导出的参数；2）崎岖度，包括成对上位效应和局部最优数。所有值均为经验推导和计算得出。未插补缺失值，假设其遵循与现有数据相同的分布，因此不影响属性计算。所有值可在数据存储中找到，所有实现可在SSMuLA代码库中找到。
活性变体定义
对于包含终止密码子变体适应度数据的景观，"活性"变体定义为高于所有含终止密码子序列（预期无活性）平均适应度1.96个标准差的变体。对于GB1、T7和TEV，遵循作者基于其适应度测量系统检测限设定的截断值。
适应度统计
使用SciPy Python包中的"统计函数"（`scipy.stats`）和信号（`scipy.signal`）模块计算峰度、估计柯西峰位置并确定核密度估计（KDE）峰数。具体而言，峰度使用`stats`模块中默认设置的`kurtosis`函数计算。柯西峰位置使用`stats`模块中`cauchy`分布对象的`fit`方法估计。KDE峰数通过`stats`模块中的`gaussian_kde`函数估计概率密度函数，然后使用信号模块中的`argrelextrema`函数识别局部最优值来确定。
成对上位效应计算
将成对上位效应分为三类：幅度、符号和互惠符号。对每个活性变体，为所选位点的每个可能的双替换分配上位效应类型。然后计算每个起始变体在景观中每种上位效应类型的比例。为增强与DE可导航性的相关性，将加性相互作用纳入幅度上位效应，将符号和互惠符号上位效应合并为非幅度上位效应。缺失值被省略。对每个独特的起始变体ab，设ω_ab为起始变体的适应度值，ω_Ab为仅在A位点发生氨基酸替换的变体适应度值，ω_aB为仅在B位点发生氨基酸替换的变体适应度值，ω_AB为在A和B位点均发生氨基酸替换的变体适应度值。
幅度上位效应：两个氨基酸替换的联合效应大于或等于各自突变方向上的加性效应。这可通过单步或重组型DE方法导航。幅度上位效应的比例通过简单计数归类为幅度上位效应的实例数除以成对相互作用总数来计算。
非幅度上位效应：包括符号上位效应和互惠符号上位效应。符号上位效应中，一个氨基酸替换在另一个替换存在时的效应方向改变，使得替换顺序影响单步DE可导航性。互惠符号上位效应中，联合效应在彼此存在时改变两个突变的方向。这无法通过本质上是贪心上坡行走的单步DE实现。非幅度上位效应的比例简单计算为1减去幅度上位效应的比例。
成对上位效应与C-α距离的相关性
突变残基的成对C-α距离基于每个亲本结构（PDB: 6X0A, 5CEG, 2GI9, 6XG5, 1CEZ, 1LVM和8VHH）计算，然后对每个景观取平均值。随后计算景观中残基间平均C-α距离与成对非幅度上位效应比例之间的Spearman相关性。对景观内每个残基对(i, j)，C-α距离d_ij计算为：d_ij = ‖r_i^Cα − r_j^Cα‖，其中r_i^Cα和r_j^Cα为残基i和j的C-α原子的笛卡尔坐标，‖·‖表示欧几里得范数。
【数据】活性变体阈值：1.96个标准差；归一化公式：ω′ = ω/ω_max；1%阈值来源：(1/(4×20))×100%；结构PDB编号：6X0A, 5CEG, 2GI9, 6XG5, 1CEZ, 1LVM, 8VHH；上位效应分类：幅度、符号、互惠符号
研究结果
景观概览
本研究选取了16个实验组合景观，涵盖一系列结合相互作用和酶活性。所有景观的突变均位于结合相互作用点、活性位点或先前显示可调节适应度的位置，这些位置通常是工程改造的目标。结合方面，考察了两个三位点细菌毒素-抗毒素ParD-ParE景观和用于免疫球蛋白结合的GB1景观。酶活性方面，分析了一个三位点二氢叶酸还原酶（DHFR）景观、一个三位点T7 RNA聚合酶景观、一个四位点TEV蛋白酶景观，以及十个三位点或四位点热稳定色氨酸合酶β亚基（TrpB）景观。
首先对景观进行表征，为后续评估提供解释。为最小化理论景观建模与实验应用之间的错位，选择了两组直接从数据集导出的经验性可解释属性：（1）适应度统计（不包含序列信息）和（2）景观崎岖度（涉及将序列映射到适应度）。使用以下统计量间接推断适应度景观的复杂性：活性变体百分比、柯西峰位置对应的适应度值、适应度分布的峰度（"尾重程度"）和核密度估计（KDE）峰数。柯西分布以其重尾特征著称。使用其峰位置对应的适应度作为景观属性，以捕获大多数变体的适应度。KDE是一种非参数方法，用于估计适应度分布的概率密度函数。KDE不假设适应度具有任何特定潜在分布，有助于平滑噪声。KDE峰数反映适应度分布的模态性，可作为潜在景观可导航性的代理指标，影响DE的结果。
为量化对DE构成导航挑战的崎岖度，纳入了局部最优数和成对非幅度上位效应百分比。局部最优定义为适应度高于所有相差单个氨基酸替换的活性邻居的变体。近期研究也强调了上位效应对DE的影响，并指出大多数上位效应是成对的。因此，也将非幅度成对上位效应的量作为相关景观属性纳入。
所有MLDE策略一致优于DE，尤其在景观可导航性降低时
首先评估景观属性如何影响蛋白质工程策略的效果。具体而言，使用两个指标评估蛋白质工程活动的成果：（1）"平均最大适应度"，即每种方法平均达到的最终变体适应度；（2）"达到全局最优的比例"，衡量达到真实最大适应度的频率。在MLDE、ALDE、聚焦训练和三种不同DE策略中考察了这些指标。
DE策略包括："重组"（recomb），即每个位点最佳SSM变体的重组（n_total = 19 × n_site + 2个变体，n_round = 2轮）；"单步"（single-step），从任意位点开始的迭代过程，后续变体基于已发现的最佳变体构建（n_total = 19 × n_site + 1个变体，n_round = n_site轮）；"top96重组"（top96 recomb），在每个位置进行SSM，基于加性重组计算所有替换组合，选择前96个变体（n_total = 19 × n_site + 97个变体，n_round = 2轮，其中96为常用筛选板孔数）。
MLDE和ftMLDE训练模型集成，采用随机或ZS预测器引导的训练样本选择。训练后的模型用于预测所有变体的适应度，前96个预测变体用于评估。ALDE将总样本量平均分为多轮（n_round = 2、3或4），每轮采样由采集函数引导。类似ftMLDE对MLDE的改进，ftALDE使用ZS预测器选择信息量更大的初始训练集，而非ALDE中使用的随机采样。
考虑到实验筛选的通量和费用差异，探索了筛选的唯一变体总数（n_total）范围，从120到2,016个样本。跨景观平均而言，MLDE需要48个训练样本（n_total = 144）才能超越重组DE，需要96个（n_total = 192）才能在两个指标上超越单步DE。MLDE需要96个训练样本（n_total = 192）才能匹配平均最大适应度，需要384个（n_total = 480）才能达到与最具竞争力的DE策略（top96重组）相当的全局最优比例。通过整合各种ZS预测器，ftMLDE始终优于随机采样的MLDE（在最多960个训练样本（n_total = 1,056）时，平均最大适应度提高4%–12%；在所有训练样本量下，全局最优比例提高9%–77%），且ftMLDE以更少的训练样本达到与MLDE相同的平均最大适应度和全局最优比例水平。这些结果表明，MLDE能比DE更有效地鉴定高适应度变体，且其性能随训练数据增加而提高。
使用ZS预测器的聚焦训练可进一步改善性能。接下来，将ALDE（多轮训练和测试，由采集函数引导）与相同总筛选变体数的MLDE（单轮训练和测试，相当于两轮）进行比较。两轮时，ALDE在总样本480个后开始在平均最大适应度上超越MLDE，在288个总样本后开始在全局最优比例上超越MLDE，但直到1,056个样本才在两个指标上超越ftMLDE。四轮时，ALDE匹配或超过ftMLDE性能。使用聚焦训练时，ftALDE在相同轮数（n_round = 2）下匹配或超过ftMLDE，并随轮数增加进一步改善（ftALDE × 3和ftALDE × 4）。然而，对于活性变体少于1%的文库，即使四轮ALDE（无聚焦训练）也始终不如ftMLDE。这些观察强调了使用ZS预测器进行聚焦训练的效用，使MLDE能匹配多轮ALDE的性能，并进一步改善ALDE。
鉴于不同方法跨景观性能的标准差较大，首先使用Elo评分系统验证比较结果——一种广泛用于模型基准测试和竞技游戏的方法。然后检查每种方法在单个景观上的表现，发现某些景观比其他景观表现出更显著的改进。
【数据】景观数量：16个；n_total范围：120–2,016；MLDE训练样本：48（n_total=144）超越重组DE，96（n_total=192）超越单步DE；ftMLDE改进：平均最大适应度提高4%–12%（最多960训练样本），全局最优比例提高9%–77%；ALDE超越MLDE阈值：480总样本（平均最大适应度）、288总样本（全局最优比例）；ALDE超越ftMLDE阈值：1,056样本
讨论与解读
本研究证实，所有测试的MLDE策略在16个蛋白质适应度景观上均超过或至少匹配DE性能，且当景观属性对DE构成更大障碍时（如活性变体更少、局部最优更多），优势更加明显。使用利用各种先验知识的ZS预测器，聚焦训练中富集的训练集可带来进一步性能提升。总体而言，这项计算研究表明MLDE策略具有高度泛化性，能比DE实现更好的蛋白质工程结果，并提出了有效部署这些方法的关键考量。
需要指出的是，公平比较传统DE方法和ML策略具有挑战性。基于迭代SSM的DE方法固有地搜索空间有限。评估的重组方法仅组合了相对起始序列含单个氨基酸替换的变体。单步DE随着目标位点数量增加需要更多轮工程，拖慢进程。此外，虽然使用唯一变体评估DE，但实践中确保足够的变体空间覆盖率（如95%）通常需要接近3倍过采样。相比之下，ML辅助方法可避免这些限制。
编译者解读
这项研究最值得称道之处在于用16个真实实验景观而非模拟数据来检验MLDE策略，结论对湿实验有直接参考价值。聚焦训练与零样本预测器的结合，本质上是把进化信息、结构信息和稳定性知识"前置"到模型训练中，在数据稀缺时尤其有效。不过，研究依赖的计算模拟评估与真实多轮湿实验仍有距离——实际筛选中的实验噪声、漏检和成本约束可能改变策略排序。对合成生物学实践者而言，一个可操作的启示是：当景观崎岖度未知时，优先选择ftMLDE或ftALDE策略，并用活性变体比例作为快速判断指标。
参考来源
Li F-Z. Evaluation of machine learning-assisted directed evolution across diverse combinatorial landscapes. Cell Systems, 2025. DOI: 10.1016/j.cels.2025.101387
DOI: 10.1016/j.cels.2025.101387