基于机器学习的模拟：通过深度突变扫描实验解析SARS-CoV-2的适应性特征

《Journal of Chemical Information and Modeling》：Machine Learning-Driven Simulations of the SARS-CoV-2 Fitness Landscape from Deep Mutational Scanning Experiments

　　高分辨率图像
下载MS PowerPoint幻灯片

预测蛋白质变异效应是在准备病原病毒株、理解与突变相关的疾病以及设计新蛋白质方面面临的关键挑战。由于复杂的变构效应和上位效应，蛋白质序列-结构-功能关系难以建模。为了研究有效的建模策略，我们使用标记了血管紧张素转换酶2（ACE2）结合亲和力的SARS-CoV-2受体结合域（RBD）序列深度突变扫描（DMS）库训练了监督机器学习（ML）模型。与仅仅添加或平均点突变效应相比，这些模型在预测组合突变效应方面表现出更优异的性能，并且在仅使用接近野生型（WT）变异体进行训练时，它们在排版奥密克戎变异体方面也具有很强的外推性能。我们通过结合ML和马尔可夫链蒙特卡洛模拟来描述RBD适应性景观，以从WT序列预测进化模式。这些模拟生成的序列谱与DMS数据中的高适应性序列相当，并能预测未见的奥密克戎变异体中的突变。这些模型为RBD序列元素之间的关系提供了见解，并为使用DMS预测新兴病毒株提供了新的视角，我们预计这将适用于其他进化预测任务。为了便于应用和未来开发这一策略，我们介绍了Mavenets：https://github.com/SztainLab/mavenets

引言
导致COVID-19大流行的SARS-CoV-2病毒仍在不断进化，积累能够逃避免疫反应和治疗的突变。(1?4) 由于可能引发毒力菌株的突变数量众多，为未来的变异株做准备仍然具有挑战性。大多数问题突变发生在刺突蛋白的受体结合域（RBD）中，该区域作为病毒感染的第一个接触点与血管紧张素转换酶2（ACE2）结合。(5) Starr等人在2020年对刺突RBD进行了深度突变扫描（DMS)，(6) 提供了105个RBD序列变异体对ACE2结合影响的序列-功能信息，这是病毒适应性的一个强指标。(7?9) 随后在2021年从B1.351变异体开始进行了额外实验(10)，以及在2022年从奥密克戎BA.1和BA.2开始进行了实验(11)，这些实验使用了优化过的库生成方法，在每个位置产生饱和点突变，每个实验仅描述了较少数量（104个）的额外序列。然而，这些研究本质上无法提供关于所有可能突变的全面信息（例如，对于一个包含200个残基的蛋白质，约有20^260种可能的突变）。此外，任何单个序列的报告适应性估计都存在噪声，这阻碍了重要突变体的分离，并妨碍了对单个突变的理解。

能够计算预测任何突变对病毒适应性的影响将使科学家能够迅速识别出未来可能成为问题的突变。然而，这种依赖性很难建模。由于蛋白质序列-功能关系的复杂性，仅通过添加或平均单个突变的影响无法准确估计组合突变的效果。(12?15) 机器学习（ML）在建模复杂相互依赖性方面取得了显著成功；(16) 例如，ML在分子生物学中具有革命性影响，例如Alphafold方法可以预测序列-结构关系。(17?21) 同样，基于蛋白质语言模型（PLMs）的方法，如ESM(18)和ProtT5(22)通过利用大规模的蛋白质序列数据集学习了通用的序列-功能关系。此外，包括EVE(23)和GEMME(24)在内的共同进化模型使用了多个序列对齐的信息来推断突变效应。最近，结合PLM、共同进化或显式结构表示的混合方法在专门的蛋白质预测任务中显示出了改进。(25?29) 虽然这些方法旨在具有广泛的适用性，但它们对高度专业化或超出分布范围的情况的泛化能力仍然存在变化，(30,31) 通常需要通过微调来进行领域特定的调整。(32?35) 此外，这些大型模型所需的庞大计算资源在某些应用中可能带来实际限制。

对于定义明确、特定于任务的应用，使用专注于特定数据训练的专用模型可能具有优势。利用多重效应分析（MAVE）(36)生成的高通量数据（如深度突变扫描（DMS）(37)大大促进了专用蛋白质适应性预测器的发展。(37) 例如，Gelman等人证明了神经网络可以从DMS数据中学习序列-功能关系，捕捉非线性效应，并指导蛋白质设计。(38) Tareen等人引入了MAVE-NN，它结合了统计和神经网络方法来从MAVE数据中建模基因型-表型关系。(39) Faure等人开发了MoCHI，这是一个灵活的框架，它将神经网络与可解释的统计指标相结合，独特地支持高阶上位性建模和多模态表型数据的整合。(40) 这些方法已被应用于学习SARS-CoV-2 RBD适应性。Taft等人描述了一种“深度突变学习”方法，该方法整合了DMS数据来合理设计针对特定感兴趣区域的突变库。(15) 在这些组合库上训练的机器学习分类器能够准确描述RBD变异体作为ACE2和各种治疗抗体的结合物或非结合物的特性。重要的是，他们发现，当编辑距离远离野生型（WT）时，点突变并不是additive的，并且组合库对于建模性能至关重要。正如作者指出的，这项研究仅关注已知在关注变异体（VOC）中突变的选定位点的组合突变，限制了预测范围。

许多其他人也开发了SARS-CoV-2刺突适应性的预测模型；(41?51) 然而，很少有模型关注在早期菌株中序列和实验数据很少的情况下进行早期预测的能力。(52) 在这项工作中，我们使用DMS库训练了RBD适应性预测器，并评估了它们预测未见序列适应性的能力。准确性是通过使用用于训练的DMS实验中的保留数据、完全保留的DMS实验的数据，以及通过使用适应性预测器的输出作为能量函数并通过马尔可夫链蒙特卡洛（MCMC）采样的适应性景观来表征的。

在从WT序列开始的大型组合DMS库（105）上训练的模型被发现能够准确描述用于训练的实验中保留数据的趋势，并能够准确预测保留的DMS实验中存在的模式，包括那些编辑距离远大于训练中发现的序列，如奥密克戎菌株。在专注于实现饱和点突变的较小库上训练的模型显示出较差的预测能力。通过评估它们预测种群中观察到的突变进化的能力来评估适应性景观模拟。所生成的序列谱优先考虑了VOCs中的热点残基位置和氨基酸替换，这与通过对组合DMS库应用适应性阈值所识别的谱相当。尽管模拟生成了多样的组合序列，但如何将这些景观转化为流行病学预测仍是一个未解决的问题。总的来说，这些发现表明，在多突变数据上训练预测器的重要性，以及ML比从点突变推断简单模型更准确地描述高阶突变的能力，但也突显了在序列空间中评估预测器性能的困难。这里介绍的工作流程在Mavenets中实现，涉及在DMS数据上训练和评估定制的机器学习模型以生成蛋白质适应性景观。我们预计这个框架可以容易地扩展到广泛的突变数据集的进化预测。

结果和讨论
机器学习架构
预测RBD结合亲和力
首先通过在同一DMS实验中提取的数据训练和评估模型来量化预测性能。通过最小二乘回归训练了多个预测器，以将RBD序列与ACE2结合亲和力联系起来。用于训练的DMS数据集(6)在实验上以野生型（WT）RBD序列初始化，包含105,525个变异体及其相应的相对KD值，以Δlog KD相对于WT报告，即高于零的值比WT结合得更好。

预测器包括线性模型、全连接神经网络（MLP）带有一hot嵌入和protT5蛋白质语言模型（PLM）嵌入(16,53)、消息传递图神经网络（MPGNN）(54,55)、变压器(53,56)以及带有dropout和不带dropout的梯度提升决策树(DART)(57,58)；(59) （图1）中显示了部分表现优异的结果。在所有搜索的模型和超参数（表S1–S3）中，使用一hot序列嵌入的浅层MLP实现了最低的平均绝对误差（MAE）0.27 Δlog KD，其性能优于线性模型、DART、T5嵌入、变压器和基于结构的MPGNN，后者的测试MAE分别为0.93、0.46、0.34、0.29和0.38 Δlog KD。尽管可能存在更优的超参数和嵌入，但简单的一hot嵌入和MLP架构在这个数据集上的强大性能对于使用这些模型进行下游任务来说是令人鼓舞的，因为它们具有复杂的表示。这些指标是使用包含总数据10%的随机选定的测试集计算的。为了确保测试集在训练过程中完全不受影响，还使用了相同大小的单独验证集进行超参数优化。对于带有一hot嵌入的MLP，除了完全重新划分数据外，还测试了两种额外的随机训练和验证分割，每种方法都显示出相似的性能（图S1）。

图1
图1. 在保留的测试集上各架构的性能总结。预测的Δlog KD值与实验测量的Δlog KD值一起绘制，并在每个图中标记了MAE和Spearman的秩系数。

图1. 如果通过合适的图表定义，表达力足够的MPGNN应该能够胜过MLP。测试的MPGNN是使用在WT结构中发现的最近邻关系创建的。预测蛋白质变异体的结构是一个突出的挑战，因为关于蛋白质变异体的结构真实信息非常少。结构预测器通常无法预测单个突变的结构影响，因为它们依赖于已解决的结构和观察到的进化数据。(60,61) 此外，许多突变会阻止成功折叠，这对结构预测构成了挑战。尽管设计模型已被证明可以成功识别稳定突变，但它们是通过识别符合蛋白质数据库（PDB）的结构分布的序列来实现的，而不是明确预测单个突变的结构。(62) 如Gelman等人所示，(38) 单个野生型结构可以高精度地用于图边的定义。然而，与Wang和Gamazon(44)类似，我们没有在这个数据集上看到图表示的性能提升。此外，由于结合亲和力是集合属性，静态结构信息可能无助于预测突变Δlog KD。

带有一hot嵌入的MLP预测高阶突变的适应性，其准确性远高于简单地添加或平均点突变的影响。虽然这并不令人惊讶，但它强调了在分析DMS数据时考虑非添加效应的重要性，以及生成组合突变库的价值。(12?15) 不管突变数量如何（最多为9，即训练集中存在的最大数量，图2A），MLP的MAE保持较低且相对一致。此外，我们测试了训练期间存在的组合突变数量对模型性能的影响，并发现了直接的相关性（图S2）。然而，随着突变计数的增加，数据点数量减少，尽管保持了较低的MAE，但Spearman的秩相关系数在超过六个WT突变时显著下降，这可能是因为这些突变倾向于都落在低结合、高噪声区域（图2B,C）。正如之前所指出的，(63,64) 验证变异效应预测器的准确性是一项复杂任务。预期的使用目的不同可能会改变哪些方面的性能相关；因此，考虑多个指标是至关重要的。

图2
图2. 按WT中的突变数量分层的预测。（A）盒线图显示了使用MLP带有一hot嵌入、线性模型或通过对序列中的每个点突变求和或平均值得到的预测序列与真实Δlog KD之间的绝对误差。（B）A中每个预测的Spearman秩相关性。注意，求和和平均值保持相同的排名，因此显示为单条青线。（C）库中包含每种WT突变数量的序列数量。

高分辨率图像
下载MS PowerPoint幻灯片

几个VOC在疫情早期就超过了WT菌株。(65) 因此，随后的RBD–ACE2结合研究使用了不同的起始序列。2021年，Bloom及其同事发表了从N501Y、E484K和B1.351（N501Y、E484K和K417N）衍生出的额外库（我们将其称为“beta研究”），随后在2022年又发表了从奥密克戎BA.1和奥密克戎BA.2数据集衍生出的库（我们将其称为“奥密克戎研究”）。这些结合测量遵循了原始WT研究的设计；然而，新的库构建技术旨在高效饱和点突变，导致每个库的数据量减少了大约一个数量级（约10^4）。这些新研究都包含一个从WT序列构建的对照库，因此在原始研究的基础上新增了七个新库。每个库都使用上述相同的方法训练了具有独热嵌入的MLP。性能是在每个实验特定的随机保留测试集上评估的。在这些较小数据集上训练的模型的实验内性能相对较低，MAE范围为0.57至0.64 Δlog KD，而使用原始数据集训练的模型的MAE为0.27 Δlog KD；Spearman系数范围为0.730至0.825，而原始研究的Spearman系数为0.896。这些结果来自单独的超参数扫描，以获得每个数据集的最优网络架构（图3A）。

图3：初始模型对其他7个DMS库的外推性能总结。每个数据集都根据用于生成库的起始序列进行了标记。这些研究由Starr等人独立进行，每个研究都有一个武汉（WT）参考序列对照。研究包括原始（基础）、（6）B.1.351（武汉beta、beta、N501Y、E484K）、（10）和奥密克戎（武汉奥密克戎、BA.1、BA.2）。（11）（A）独立训练模型的预测Δlog KD值与实际Δlog KD值。（B）原始模型在每个单独库的测试集上的预测Δlog KD值与实际Δlog KD值。

由于我们在原始WT库上训练的模型使用独热嵌入时，在距离WT序列最多九个突变的情况下显示出相对较低的MAE，并且在距离WT序列最多五个突变的情况下显示出较高的Spearman系数（图2A、B），我们试图确定仅使用原始WT库训练的模型预测独立实验生成的数据的能力。Beta实验包含从WT衍生出最多七个突变的库，而奥密克戎实验包含从WT衍生出14到19个突变的库。令人惊讶的是，与使用同一实验数据训练的模型相比，使用原始WT训练的模型得到的预测Spearman相关系数更高。我们注意到这对于独热嵌入是成立的，并且这些额外数据集可能存在更优的嵌入和超参数。例如，使用T5嵌入训练的BA.1模型表现出更好的性能，MAE为0.52 Δlog KD，Spearman系数为0.828（图S3）。尽管如此，我们得出结论，我们的模型在原始数据集上的外推性能高于预期，能够准确地对独立实验中的突变进行排序，包括那些突变总数超过训练集的突变。

为了进一步研究独热嵌入和T5嵌入的外推性能，我们比较了在原始数据集或BA.1数据集上训练的模型使用这两种嵌入方式的预测能力（图S3）。在原始数据集的情况下，使用独热嵌入训练的模型的Spearman相关系数更高；然而，在仅使用BA.1数据集训练的模型中，T5嵌入在七个额外数据集中的四个数据集中表现出更好的外推性能。我们假设独热嵌入的惊人高外推性能可能是由于DMS数据集包含与WT具有高度序列一致性的单一蛋白质，而用于训练PLM的大量蛋白质序列可能在区分这些细微差异方面信息量较少。尚未确定微调、选择不同的PLM或修改数据准备策略是否可以进一步提高性能。观察到的变异性表明，不应仅基于单次实验的性能观察来断定哪种架构对所有DMS数据集都是最优的。结合速度和准确性的独热模型特别适用于下游应用，包括下一节中描述的马尔可夫链蒙特卡洛采样。因此，除非另有说明，以下实验都是使用具有独热嵌入的模型进行的。

我们关于超出分布性能的结论基于Spearman相关系数，该系数表明模型在排序突变方面的能力更强，而不是基于MAE。绝对测量值（如KD）对实验条件非常敏感，并且实验之间存在显著差异。因此，在分析多个DMS实验时，通常会进行数据转换以考虑实验响应曲线中的差异，这通常使排名成为更稳健的适应性指标。视觉检查表明，较大的MAE是由于实验特定的偏移造成的，其中WT训练的模型在所有库中都高估了亲和力，而仅在奥密克戎库中低估了亲和力（图3B）。通过实验特定的线性校准可以将MAE降低到0.26至0.56 Δlog KD的范围，而不影响排名（图S4）。然而，这种校准需要预先知道实验偏移，而这往往是不易获得的。

自然会想知道，通过在多个实验中汇集训练数据是否可以改善WT训练的网络在其他较小数据集上的强大外推性能。首先，尝试使用从小实验中提取的数据对预训练的WT网络进行微调；这并未显著提高低数据实验的性能。相反，从头开始使用独热嵌入在所有实验的组合数据上训练MLP被证明更有效，从而将MAE范围降低到0.38–0.46 Δlog KD，并将Spearman相关系数范围提高到0.893–0.907（对于所有数据集），除了从WT衍生的库，其性能下降到0.47–1.02 Δlog KD MAE和0.831–0.864 Spearman相关系数（图4A）。从WT衍生的各种库在其序列上有显著重叠，但根据每个实验设置有不同的标签。为了考虑实验响应曲线可能的不同，我们在MLP架构中添加了一个调整头，根据记录的实验修改序列依赖的预测。（66,67）这使WT库的MAE降低到0.28–0.41 Δlog KD，并将Spearman相关系数范围提高到0.892–0.919。非WT库在这个调整后没有显著变化（图4B）。

我们注意到，对于某些分析来说，实验特定的调整可能很有用，但使用该模型来预测保留实验的结果并不直接，因为相应的响应曲线是未知的。正如没有实验调整的模型所示，当这些实验在序列空间中表征不同位置时，仅凭MLP本身就具有捕捉多个实验结果的能力。理想情况下，这些模型的预测准确性应该通过与不受实验响应曲线影响的“真实”结合亲和力进行验证；然而，作者没有这样的测试集。相反，准确性是基于预测病毒适应性的流行病学趋势的能力来评估的，如下一节所述。

为了将我们训练的预测器与观察到的群体突变趋势联系起来，使用具有独热嵌入的训练好的MLP准确预测了RBD结合亲和力，从而利用MCMC探索了未来VOCs的估计适应度景观。（68）使用结合亲和力的预测创建了适应度表面。首先，模拟旨在仅从WT序列的信息评估预测VOC突变趋势的能力，因此使用了原始的WT训练预测器。由于在用于训练的原始DMS中没有VOC的知识，这代表了对外推性能的严格评估。

首先，我们提出可以通过假设系统具有类似玻尔兹曼的分布??(??)∝??(??)=e???Δlog??D(??)来定义一个信息量丰富的景观，其中p(x)是特定序列的概率，f(x)与这个概率成正比，β是一个未知的标量，Δlog KD是给定突变的真实结合亲和力变化。（1）在这个公式中，Δlog KD作为相对结合自由能变化（ΔΔG）的代理。通过定义类似玻尔兹曼的分布，我们将高亲和力变体（Δlog KD > 0）视为在进化景观中占据较低能量、更可能的状态。这个定义意味着使用训练模型的预测输出代替Δlog KD进行MCMC，其中不同的β值会集中采样更稳定的突变体。由于MCMC允许系统探索与训练数据中存在的序列高度不同的序列，这种类型的模拟代表了对外推能力的复杂评估。为了从方程1的分布中采样，通过选择1到201之间的随机位置来随机突变一个20个标准氨基酸中的一个。提出的突变体根据Metropolis标准被模型评估并接受或拒绝。尽管在实验外排名性能很有希望，但使用方程1进行的所有模拟结果都产生了完全突变的序列（每个位置都包含与WT不同的氨基酸），表明在高突变水平下模型的准确性从根本上是有限的（图S5）。我们还注意到，组合序列空间随着序列长度的增加而呈指数级增长；因此，随机突变提议实际上偏向于高度突变的变体，而不是少数位点的突变。

为了产生更具信息量的适应度模拟，进行了两项修改。首先，为了避免过度突变导致采样超出训练数据中的典型突变水平，禁止采样序列超过WT序列的9个突变。然而，所得到的模拟仍然累积了最大设定数量的突变，很少有采样低于7个突变（图S5），这表明对于WT的9个突变范围内的特定突变，可能存在适应度预测的错误。因此，为了产生物理上合理的突变统计并进一步集中采样，我们进一步修改了采样分布，使其与方程2成正比：????(??)=1[??(??)>??]e?????(??)?????(??)fλ(x)=1[N(x)>ω]e?βg(x)?λN(x)（2），其中g对应于预测的亲和力，??N是与WT序列中的残基匹配的残基数，λ控制向WT序列的偏置，1表示0–1的指示函数，ω限制了模拟期间允许的最大突变数量。λ可以被视为控制模拟期间发现的平均突变数量的调整参数。方程2通过增加在每个位置提出WT氨基酸的概率来采样，同时保持提出非WT氨基酸的相等概率，类似于Rosetta中的“偏好天然残基”功能（图S5）；这种修改的提议机制控制了马尔可夫链的平稳分布。

我们注意到，实验特定的调整可能对某些分析有用，但使用该模型来预测保留实验的结果并不直接，因为相应的响应曲线是未知的。正如没有实验调整的模型所示，当这些实验在序列空间中表征不同位置时，单独的MLP具有捕捉多个实验结果的能力。理想情况下，这些不同模型的预测准确性应通过与不受实验响应曲线影响的“真实”结合亲和力进行验证；然而，作者并没有这样的测试集。相反，准确性是基于预测病毒适应性流行病学趋势的能力来评估的，如下一节所述。

为了将我们训练的预测器与观察到的群体突变趋势联系起来，使用具有独热嵌入的训练好的MLP来探索使用MCMC预测的未来VOCs的估计适应度景观。（68）使用多种定义根据结合亲和力的预测来创建适应度表面。首先，模拟旨在评估仅从WT序列的信息预测VOC突变趋势的能力，因此使用了原始的WT训练预测器。由于在用于训练的原始DMS中没有VOC的知识，这代表了对外推性能的严格评估。首先，我们提出可以通过假设系统具有类似玻尔兹曼的分布??(??)∝??(??)=e???Δlog??D(??)p(x)∝f(x)=e?βΔlogKD(x)来定义一个信息量丰富的景观，其中p(x)是特定序列的概率，f(x)与这个概率成正比，β是一个未知的标量，Δlog KD是给定突变的真实结合亲和力变化。在这种表述中，Δlog KD作为结合自由能变化的代理。通过类似玻尔兹曼的分布定义景观，我们将高亲和力变体（Δlog KD > 0）视为占据进化景观中较低能量、更可能的状态。这个定义意味着使用训练模型的预测输出代替Δlog KD进行MCMC，其中不同的β值会集中采样更稳定的突变体。由于MCMC允许系统探索与训练数据中存在的序列高度不同的序列，这种类型的模拟代表了对适应度预测器外推能力的复杂评估。为了从方程1的分布中采样，通过在1到201之间选择一个随机位置来提出随机突变到20个标准氨基酸中的一个。提出的突变体由模型评估并根据Metropolis标准接受或拒绝。尽管在实验外的排名性能很有希望，但使用方程1进行的所有模拟都产生了饱和突变的序列（每个位置都包含与WT不同的氨基酸），这表明在高突变水平下模型的准确性基本上是有限的（图S5）。我们还注意到，组合序列空间随着序列长度的增加而呈指数级增长；因此，随机突变提议实际上偏向于高度突变的变体，而不是少数位点的突变。

为了产生更具信息量的适应度模拟，进行了两项修改。首先，为了避免过度突变导致采样超出训练数据中的典型突变水平，禁止采样序列超过WT序列的9个突变。然而，所得到的模拟仍然累积了最大设定数量的突变，很少有采样低于7个突变（图S5），这表明对于WT的9个突变范围内的特定突变，可能存在适应度预测的错误。因此，为了产生物理上合理的突变统计并进一步集中采样，我们进一步修改了采样分布，使其与方程2成正比：????(??)=1[??(??)>??]e?????(??)?????(??)fλ(x)=1[N(x)>ω]e?βg(x)?λN(x）（2），其中g对应于预测的亲和力，??N是与WT序列中的残基匹配的残基数，λ控制向WT序列的偏置，1表示0–1的指示函数，ω限制了模拟期间允许的最大突变数量。λ可以被视为控制模拟期间发现的平均突变数量的调整参数。方程2通过增加在每个位置提出WT氨基酸的概率来采样，同时保持提出非WT氨基酸的相等概率，类似于Rosetta中的“偏好天然残基”功能（图S5）；这种修改的提议机制控制了马尔可夫链的平稳分布（见方法部分）。模拟序列与自然进化序列的比较

从fλ中抽取的序列与截至2025年1月31日存入Genbank的人类群体中的真实序列分布进行了比较。这些Genbank序列代表了病毒的总体适应性，而模拟分布仅依赖于预测的ACE2结合亲和力和与野生型的序列相似性。免疫逃逸是SARS-CoV-2进化的主要驱动力。由于中和抗体与ACE2竞争结合RBD，高适应性变异体通常会进化以降低抗体亲和力，同时保持或增加对ACE2的亲和力。(1?4)因此，ACE2结合亲和力适应性为潜在传染性提供了一个合理的代理指标。然而，这没有考虑到那些虽然能够强烈结合ACE2但因为维持或增强了抗体结合而适应性不高的序列。此外，具有高ACE2亲和力的预测突变可能对病毒生命周期的其他阶段产生负面影响，例如成熟病毒颗粒的复制和释放，而这在我们的模型中并未得到体现。尽管这些情况可能导致对病毒总体适应性的误判，我们仍然尝试将我们的ACE2结合适应性景观与当前在人群中检测到的序列进行比较。

用于训练的DMS数据中的突变频率也与Genbank数据进行了比较。由于这些DMS数据包含的突变是随机的，并不依赖于ACE2结合亲和力，而是反映了所使用的实验程序，因此根据不同的Δlog KD分数阈值选择了样本（图S6）。实验中的Δlog KD阈值为零，表示与野生型相比具有更好或相等的结合亲和力，这一阈值对于预测变异体（VOC）残基位置和突变最为有效。由于序列空间非常庞大，我们没有专注于精确的序列匹配，而是分析了每个位置的突变倾向。我们还将这种方法与仅从点突变信息分析最大结合亲和力的方法进行了比较，排除了DMS数据集中包含的组合序列信息。

为了确定仅基于WT训练的ML预测器是否能够预测VOC中的突变，我们抽取了基线模型给出的分布。我们发现25个M步骤足以在三个独立的MCMC复制中得到一致的突变谱（图S7）。我们将模拟预测的前20个突变与截至2025年1月31日在Genbank中发现的突变进行了比较。超过70%的模拟预测突变发生在VOC残基上，包括位置452、493、498和501。相比之下，WT DMS数据集中超过30%的分数阈值识别的突变与VOC残基重叠，包括位置339和452（图5A）。DMS点突变中得分最高的突变仅包含VOC位置498、501和505。许多残基位置被预测会突变为多种氨基酸。为了考虑到这一点，我们检查了突变频率最高的20个残基，而不考虑它们突变成了哪种氨基酸，这揭示了额外的热点位置。模拟还识别出了残基373，而DMS数据突出了残基346、493和501，点突变数据仅识别出了373、477和484。在前20个位置中，模拟预测的中有4个与VOC残基重叠，而DMS数据中有1个，点突变中有2个，尽管考虑到前20个位置，两种方法预测的热点总数相同（图S8）。总体而言，这表明我们的模拟比单独分析DMS数据更能有效地发现热点位置。

图5

图5. 序列谱分析。(A) 饼图显示了来自Genbank的序列、基于原始WT数据训练的MLP的模拟结果，以及训练数据序列中分数大于或等于0 Δlog KD（优于WT）的前20个突变。百分比对应于给定集合中具有该突变的总序列数量。楔形的颜色根据残基编号进行区分，对应于Genbank突变前20名中的残基。(B) 对于28个VOC残基的模拟序列中的突变进行Logo图表示。粉红色星号表示在给定位置上最可能的氨基酸对应于VOC突变。灰色星号表示最可能的氨基酸具有与VOC突变相同的理化性质，例如极性、碱性和芳香性。(C) A中显示的一个热编码序列的UMAP，根据Δlog KD进行着色。独特的Genbank序列用分数为零的颜色标记，分数小于零（即比WT的结合能力差）的用灰色标记。(表1)中的VOC以黑点显示，并标出了代表性的变异类别。

接下来，分析了VOC菌株中发现的28个突变在每个位置上最常见的替换情况（表1）。我们的模拟始终预测VOC氨基酸是九个位置上最常见的突变：G339D、R346T、R403K、N440K、V445P、N460K、S477N、T478K和F486P，而K356T和L368I在某些复制中被预测为最常见的突变，但在其他复制中是第二常见的（图5B和S7）。除了S477N和T478K之外，所有正确预测的突变都特定于奥密克戎菌株。此外，在其他几个位置上，模拟并未预测出VOC中发现的突变，但预测出了具有相似理化性质的氨基酸。例如，N501经常被预测为突变成F而不是Y，L452被预测为突变成K而不是R。将这一点与DMS或点突变数据相比，也预测了九个突变，然而它们的身份不同。模拟预测了R346T、S477N和F486P，但DMS阈值检测没有预测到这些突变；而N481K和N501Y则通过DMS阈值检测到但模拟没有预测到。对得分最高的点突变进行分析后，发现只有VOC位置498、501和505。许多残基位置被预测会突变为多种氨基酸。为了解释这一点，我们检查了突变频率最高的20个残基，不论它们突变成了哪种氨基酸，这揭示了额外的热点位置。模拟还识别出了残基373，而DMS数据突出了残基346、493和501，点突变数据仅识别出了373、477和484。在前20个位置中，模拟预测的有4个与VOC残基重叠，而DMS数据只有1个，点突变有2个，尽管考虑前20个位置时，两种方法预测的热点总数相同（图S8）。总体而言，这表明我们的模拟比单独分析DMS数据能更有效地发现热点位置。

图5. 序列谱分析。(A) 饼图显示了Genbank序列、基于原始WT数据训练的MLP的模拟结果，以及训练数据序列中分数大于或等于0 Δlog KD（优于WT）的前20个突变。百分比对应于给定集合中具有该突变的总序列数量。楔形的颜色根据残基编号进行区分，对应于Genbank突变前20名中的残基。(B) 对于模拟序列中28个VOC残基的突变进行Logo图表示。粉红色星号表示在给定位置上最可能的氨基酸对应于VOC突变。灰色星号表示最可能的氨基酸具有与VOC突变相同的理化性质，如极性、碱性和芳香性。(C) 图A中显示的一个热编码序列的UMAP，根据Δlog KD进行着色。独特的Genbank序列用分数为零的颜色标记，分数小于零（即结合能力比WT差）的用灰色标记。(表1)中的VOC以黑点显示，并标出了代表性的变异类别。

高分辨率图像
下载MS PowerPoint幻灯片

接下来，分析了VOC菌株中发现的28个突变在每个位置上最常见的替换情况（表1）。我们的模拟始终预测VOC氨基酸是九个位置上最常见的突变：G339D、R346T、R403K、N440K、V445P、N460K、S477N、T478K和F486P，而K356T和L368I在某些复制中被预测为最常见的突变，但在其他复制中是第二常见的（图5B和S7）。除了S477N和T478K之外，所有正确预测的突变都特定于奥密克戎菌株。此外，在其他几个位置上，模拟没有预测出VOC中发现的突变，但预测出了具有相似理化性质的氨基酸。例如，N501经常被预测为突变成F而不是Y，L452被预测为突变成K而不是R。将这一点与DMS或点突变数据相比，也预测了九个突变，然而它们的身份不同。模拟预测了R346T、S477N和F486P，但DMS阈值检测没有预测到这些突变；而N481K和N501Y则通过DMS阈值检测到但模拟没有预测到。分析得分最高的点突变后，发现F486P在模拟中被预测到，但在阈值检测中没有预测到；而T478K则在两种方法中都被预测到。与模拟类似，N501位置上的点突变到F的亲和力比到Y的亲和力更高。尽管MCMC参数β、C和ω最初是在不了解VOC的情况下选择的，但我们分析了它们对VOC热点残基和氨基酸替换预测的影响。没有发现任何参数能够改善预测（表S4）。这项分析表明，使用模拟和单独分析DMS数据来预测VOC残基上最常见的替换情况的性能相当，尽管有些突变仅由模拟预测到，而有些突变则由DMS识别但模拟没有预测到。

表1. 关注的变异体总结
名称 菌株 RBD突变
Alpha B.1.1.7 E484K, N501Y
Beta B.1.3 1K417N, E484K, N501Y
Gamma P1 K417T, E484K, N501Y
Delta B.1.6 17.2L452R, T478K
Epsilon B.1.4 2(7/9)L452R
Reta B.1.5 25E484KIota
B.1.5 26E484KIota’B
B.1.5 26S477N
Kappa B.1.6 17.1L452R, E484Q
Zeta P2 E484K
Mu B.1.6 21(0.1)R346K, E484K, N501Y
21K Omicron BA.1G339D, S371L, S373P, S375F, K417N, N440K, G446S, S477N, T478K, E484A, Q493R, G496S, Q498R, N501Y, Y505H
21L Omicron BA.2G339D, S371F, S373P, S375F, T376A, D405N, R408S, K417N, N440K, S477N, T478K, E484A, Q493R, Q498R, N501Y, Y505H
22A/B Omicron BA.4/5G339D, S371F, S373P, S375F, T376A, D405N, R408S, K417N, N440K, L452R, S477N, T478K, E484A, F486V, Q498R, N501Y, Y505H
23A Omicron XBB.1.5G339H, R346T, L368I, S371F, S373P, S375F, T376A, D405N, R408S, K417N, N440K, V445P, G446S, N460K, S477N, T478K, E484A, F486P, F490S, R493Q, Q498R, N501Y, Y505H
24A Omicron JN.1G339H, K356T, S371F, S373P, S375F, T376A, R403K, D405N, R408S, K417N, N440K, V445H, G446S, N450D, L452R, L455S, N460K, S477N, T478K, N481K, V483-, E484K, F486P, Q498R, N501Y, Y505H

用粗体显示了模拟预测的突变。奥密克戎变体的名称也用Nextstrain分支方案进行了标记。(80)虽然这里的分析侧重于基于序列的预测，但需要注意的是，不同RBD序列的生物学效应是由RBD和ACE2之间相互作用的物理差异驱动的。VOC的结构、能量和理化因素已在其他地方进行了广泛研究。(74?79)然而，我们强调了仅使用WT序列知识进行VOC预测的基于序列的工作流程具有显著的优势（例如，在研究结构尚不清楚的新出现病毒时）。

由于本文中的模型是训练用于最小化保留数据上的误差，因此它们的预测可能会消除DMS测量中固有的噪声。这引出了一个问题：用预测值替换实验观察结果能否改善DMS数据集的阈值划分？为了评估这是否能够改善VOC突变的预测，我们在训练集上进行了推理，并使用得到的Δlog KD值重复了分析以进行分数阈值划分。这导致了与原始阈值划分相似的热点预测（图S9），表明模拟在热点预测方面仍然表现更好。然而，当考虑选定VOC位置上最常见的氨基酸时，“去噪”后的分数阈值得到的一种序列谱显示，13个最常见的替换与VOC突变正确重叠，包括模拟和实验分数阈值识别的突变，除了G339D和T478K。此外，使用去噪分数阈值还识别出了S371L、R408S、G446S、Q498R等突变，但模拟、实验分数阈值或点突变分析都没有识别到这些突变（图S10）。尽管模拟与带阈值的DMS数据的表现相似，但在回顾性地预测给定位置上最可能的突变方面，ML模型去噪后的DMS数据表现最佳。

接下来，使用基于所有DMS库训练的模型从WT开始进行了模拟，这些库包含了关于beta和omicron BA.1及BA.2菌株的信息。通过关注在BA.2之后出现的28个突变中的12个（包括BA.4/5、XBB.1.5和JN.1，表1），评估了预测结果。然而，多实验模型的模拟或相应DMS数据的阈值划分都没有识别出原始模型未预测到的新的热点或替换。原始和聚合训练数据都将BA.2之后的残基452识别为最常突变的残基，而原始数据还识别出了残基490。考虑位置而不是频繁突变时，两者都额外识别出了残基346。无论使用哪种模型训练，残基452和460都在前20名中，其中460是模拟预测的其他热点，但DMS单独训练的模拟没有预测到（图S11）。每个位置上最常见的替换与原始WT训练数据的模型相同，模拟的结果也几乎相同，不过模拟没有预测到R346T。去噪对于这个模型也没有帮助，没有识别出额外的替换，L368I、V445P、N481K、F486P的预测不正确，而L452被预测为K而不是R（图S12）。这表明包含额外的DMS库并没有改善最常见的替换预测，反而恶化了热点预测。这可能是由于库构建技术的不同导致围绕起始序列的分布较小，或者实验响应偏移的不同；然而，确切原因尚不清楚。将模拟偏差集中在BA.1序列周围也没有带来额外的热点残基或VOC突变预测。使用原始数据集训练的模型进行的模拟的预测能力不如使用WT序列作为起始数据的模拟。结合在DMS数据集上的高预测准确性和无偏模拟快速累积饱和突变的趋势表明，尽管这些预测器包含了关于序列空间的大量信息，但训练数据中未表征的某些区域给模拟带来了挑战。虽然引入额外的DMS训练数据集增加了覆盖范围，但并没有提高VOC的预测能力。这种数据量的增加被序列空间的组合增长所掩盖；例如，奥密克戎数据集提供了104个额外的数据点；然而，在其特征编辑距离内，存在超过10^18种可能的突变。此外，如使用多个数据集训练的预测器进行DMS去噪所示，在数据集上的MAE改进可能不会转化为对相关突变的更好下游预测。尽管面临这些挑战，通过将采样偏向于序列空间中置信度较高的区域，仍然可以生成包含信息性突变的模拟结果。为了更好地理解模拟结果、训练数据集以及实验观察到的序列之间的重叠情况，我们使用了主成分分析（PCA）随后进行均匀流形逼近和投影（UMAP）（81），将one-hot编码的序列投影到二维空间中。当仅使用模拟序列和VOC序列进行降维训练时，使用仅在WT数据上训练的模型，所有VOC都出现了聚集（图S14）。然而，当使用来自Genbank、DMS训练集和模拟序列的组合数据进行降维时，VOC明显聚类为三个不同的区域，包括Genbank序列中的preomicron、delta和omicron变体。模拟序列主要探索了preomicron和delta序列，同时还探索了一些新的未探索区域。DMS分数阈值化的序列似乎仅代表了preomicron序列（图5C）。VOC的覆盖范围依赖于用于降维的参考集，这凸显了诸如UMAP之类的可视化方法对输入数据选择的敏感性。尽管存在这种变异性，但用于训练的数据与模拟产生的数据之间的覆盖范围差异仍然表明了外推对模拟结果的显著影响。由于模拟准确性所需的偏差，模拟序列与omicron序列之间缺乏重叠并不令人惊讶。准确地描述序列景观仍是一个开放的研究领域，现有多种策略可供选择（82?84），为未来的研究提供了有希望的方向。

结论

我们证明了结合监督学习（ML）、DMS和MCMC的方法是模拟适应性景观并提取可能突变的有意义信息的一种可行途径。通过使用多层感知器（MLP），我们实现了对SARS-CoV-2 RBD与ACE2结合的高度准确预测。值得注意的是，这些预测在捕捉组合突变的效应方面优于加性或平均估计，表明了对上位性关系的理解。通过将预测结果与人群数据进行回顾性比较，该预测器在与最大熵偏差结合使用时，能够近似地描绘出一个信息丰富的适应性景观，生成了大量独特的组合RBD序列，识别出热点残基，并从仅在接近WT序列的变体上训练的模型中预测出VOC（包括omicron分支）中存在的突变。尽管识别出了独特的VOC突变，但总体序列分布并未比单独分析DMS数据中的高适应性序列发现更多的VOC突变，因此如何将高预测适应性的序列生成转化为具体的流行病学预测仍是一个未解决的问题。当前进化模拟实现的局限性在于将病毒适应性简化为单一的ACE2结合指标。此外，由于我们的模拟来源于在受控实验室中进行的DMS实验，因此它们本质上没有考虑额外的生物变量，如种群规模、选择压力缩放或环境稳定性，这些都代表了未来改进的途径。我们推出了开源软件包Mavenets：https://github.com/SztainLab/mavenets，用于重现和扩展这些发现，为开发包括多种适应性指标、替代采样算法和前瞻性变异测试的改进进化预测器奠定了基础。

方法

数据准备

WT ACE2–RBD结合DMS数据集从https://media.github.com/media/jbloomlab/SARS-CoV-2-RBD_DMS/master/results/binding_Kds/binding_Kds.csv获取。有关数据集生成的更多细节，请参见补充信息（SI）。我们过滤了包含195,081个条形码序列的原始数据集，去除了“delta_log10Ka”列中包含NaN值的序列，最终得到146,437个序列。请注意，“delta_log10Ka”指的是Δlog KD的表观值，不要与结合常数KA混淆。“表观”这一命名是为了反映研究中使用了二聚体ACE2。

每个库中有11,672个重复序列。由于Δlog KD值的标准差较低，我们取这些重复序列的平均值，最终得到108,588个序列。接下来，两个库中剩余的重复序列数量为3,063个，其标准差同样较低。对这些序列进行平均处理后，得到了本研究中使用的最终105,525个独特的序列，这些序列的Δlog KD值被用于研究。此外，Bloom及其同事针对beta（10）和omicron（11）变体进行的实验还获取了额外的数据，通过购买基因来进行位点饱和突变。合成和PCR错误不可避免地会导致序列中出现多个突变，从而使库的大小约为10^4。所有额外数据集的通用资源位置都在表S5中提供。这些数据都经过了与原始数据相同的过滤和处理流程，但Δlog KD是通过从原始库中测量得到的WT log KD值来计算的。

数据被随机划分为80%、10%和10%用于训练、验证和测试；但在DART模型的训练中，是在这些训练和验证集的组合上进行了5折交叉验证。

预测模型

多层感知器（MLP）

使用PyTorch（85）实现的多层感知器（MLP）被训练用于预测Δlog KD。网络在所有层上使用了ReLU（86）激活函数，除了输出层不需要激活函数。训练过程中使用的损失函数是预测Δlog KD值与真实Δlog KD值之间的均方误差（MSE）。批量大小为32，使用AdamW优化器（87,88）进行训练，直到通过验证损失的上升或达到最大训练周期数来观察收敛；有关训练的更多细节，请参见补充信息（SI）。优化器没有使用特定的学习计划；然而，在早期超参数调整期间尝试了无需计划的优化器（89），但没有获得更好的结果。在one-hot编码上训练的MLP具有以下不同大小的隐藏层：8、16、32、64、128和256。对于基于T5嵌入训练的MLP，还考虑了隐藏层大小为128、256、512和1024的次要网络；虽然这些较大层的MLP有所改进，但对基于one-hot的MLP并没有显著提升。通过独立改变每层的大小，创建了隐藏层深度为0–3的网络（深度为0对应于线性模型）。此外，还扫描了不同的学习率（3 × 10^–4和1 × 10^–4），并扫描了不同的权重衰减值（5 × 10^–3、1 × 10^–3、5 × 10^–4和1 × 10^–4），以及不同的dropout水平（0、0.1、0.3和0.5）。这些不同的架构和超参数通过网格搜索进行了探索（表S1和S2）。我们注意到许多MLP超参数的选择接近于最优选择的结构；因此，我们不尝试解读任何实验中发现的最优超参数。一旦确定了最优的架构和超参数设置，就用这些设置训练了3个具有不同随机初始化的网络；然后在测试集上评估了表现最好的模型实例。

序列嵌入

在早期架构探索期间，使用了多种不同的序列嵌入（即特征化方法）作为MLP的输入。基于20种典型氨基酸使用了一维编码嵌入。另一种方法是定义了一种糖胺酸one-hot嵌入，通过在one-hot嵌入中添加一个额外的标记来表示根据Nx[S/T]序列基序的糖基化天冬酰胺残基。还使用了UniRep（91）嵌入与TAPEtokenizer（92），采用了1900维的隐藏状态层。还测试了通过Bio Embeddings生成的SeqVec嵌入（93），以及ProtTrans-Bert-BFD、ProtT5-XL-UniRef50（53）和ESM1b（94）嵌入，默认设置。考虑的超参数在表S2中列出。有关输入表示的更多细节，包括用于MPN和Transformer的细节，请参见补充信息（SI）。

其他建模方法

仅使用one-hot序列嵌入进行测试，DART模型除外，它使用了基于残基类型的等效分类变量；Transformer则使用了自定义的嵌入方案。

消息传递图神经网络

使用PyG（95,96）创建了一个消息传递图神经网络，该网络基于Barros等人（97）通过模拟获得的WT RBD的单一结构进行操作。除非另有说明，否则训练细节与MLP的训练细节相同。首先使用cpptraj（98）通过分层方法并配以10的随机筛选从轨迹中提取结构，共得到了20个簇。然后使用mdtraj（99）从最庞大的簇中提取平均结构来生成残基图。每个序列使用one-hot残基表示法来描述节点，边的定义基于结构中每个残基的α碳原子之间的距离。在最优架构中，边的连接基于序列空间中N端和C端方向上10个残基的窗口，以及3D空间中15 ?的距离截止。距离通过高斯展开特征化为10维特征向量，并进行了75倍的比例缩放。消息传递了7次，每次传递都包含一个具有128个特征隐藏层的MLP。为了确定最优架构，扫描了以下超参数：窗口大小：5、10、15和20 ?；距离截止：10、15、20、30 ?；距离特征向量维度：0、6、8、10、12；缩放因子：1、10、25、50、75、100。

Transformer

在Pytorch中实现了一个BERT风格的（仅编码器）Transformer（100），使用了SiLU激活函数（101?103）。除非另有说明，否则训练细节与MLP的训练细节相同。学习率设置为3 × 10^–4。嵌入是通过根据残基类型从查找表中提取的嵌入向量与特定位置的嵌入向量相加来初始化的。这些嵌入通过多头注意力、层规范和MLP进行处理，如预LN架构中描述的那样（104）。在基于多头注意力的更新之后，每个结果嵌入都通过一个共享的MLP进行处理以产生一个标量值；这个值在整个序列上求和以得到模型预测。超参数优化是通过扫描注意力块的数量、注意力头的数量、嵌入的维度、块MLP的dropout率、注意力模块的dropout率、最终MLP的层数以及最终MLP的dropout率来进行的。扫描的值列表在表S3中提供。

多个带Dropout的加法回归树

使用Light Gradient Boosted Machine库（105）创建了多个带Dropout的加法回归树（DART）和梯度提升决策树模型（57,58）（59,58）。为了简洁起见，将没有Dropout的梯度提升决策树和带Dropout的多个加法回归树统称为DART。最优模型使用了511个叶子节点、dropout、1399棵树、0.3的dropout率、最大50的叶子数据大小、0的分类l2和5的分类平滑。通过多次局部网格搜索扫描了超参数，包括Dropout的存在（即DART或梯度提升决策树）、叶子节点的数量、树棠除掉的比例、树棠除掉的比例、分类平滑的存在以及l2。所有超参数组合都进行了5000次的训练更新，以确定最优的树木数量。我们注意到，尽管训练点较少，相对较多的叶子节点数量对回归有利，这可能暗示了高度的上位性。与其他方法不同，超参数是通过对组合的训练和验证数据集进行5折交叉验证来优化的。所有尝试都将输入视为分类变量。

调整头部

某些神经网络进行了每次实验的单独调整（图4B）。这一方法的实施首先是将序列通过内部网络转换为一个标量值，然后使用由实验身份参数化的函数对这个标量值进行精细化处理；这一过程在方程3中有描述，其中g表示完整的预测模型，h指的是内部网络，m表示实验调整。数学表达式为：??(??seq,??exp)=??(?(??seq),??exp)。??exp被定义为每个实验独有的整数，而??seq对应于蛋白质序列。对于图4B中显示的结果，内部网络被设置为多层感知器（MLP）；使用变换器架构得到的结果也类似（未展示）。实验调整被定义为一个具有一个隐藏层的MLP（输入和输出维度均为1），该隐藏层通过残差连接与h的输出相连；为了避免神经元死亡，使用了ELU（Exponential Linear Unit）作为激活函数。在非残差调整方法导致训练收敛不稳定之后，选择了残差连接方法，这阻碍了早期架构实验中的超参数搜索。选择ELU是为了避免每个实验的响应曲线出现突变。这个MLP的最终（即输出）层是特定于单个实验的，而其他参数在所有实验之间是完全共享的。在所有情况下，每个实验的调整部分的超参数都与内部网络的超参数一起进行了优化。在使用调整部分训练模型时，前20个训练周期使用了损失的线性组合：第一个损失对应于内部网络的均方误差，第二个损失对应于经过实验调整后的预测的均方误差。训练开始时这两个损失被赋予相同的权重，然后从第20个周期开始仅对实验调整后的预测进行惩罚；剩余的训练则仅使用经过实验调整后的输出结果进行。这种调度方式被观察到能够稳定训练过程。隐藏层的扫描规模分别是1、2、4、8、16、32、64、128和256。

为了比较高阶突变预测与仅使用MLP、或对点突变的影响进行加总或平均的方法，从原始野生型（WT）数据集中获得了点突变的Δlog KD值。该数据集包含了3602个独特的点突变值，覆盖了89.6%的饱和点突变情况。数据集中包含没有对应点突变的替换序列被从分析中移除。序列根据突变数量进行分组，并使用MLP、线性模型或通过对每个突变的Δlog KD进行加总或平均来进行评估。如图1所示，“线性”模型被训练为一个没有隐藏层的MLP，采用独热编码（one-hot encoding）。

为了探索训练模型输出的适应性景观，进行了马尔可夫链蒙特卡洛（Markov Chain Monte Carlo, MCMC）模拟。在MCMC模拟中，从1到201之间随机选择一个位置，将某个氨基酸突变为20种标准氨基酸中的一种。提出的突变体根据Metropolis准则被模型评估并接受或拒绝。如果使用方程1中的Δlog KD代替训练好的模型进行MCMC模拟，并采用Metropolis-Hasting（107）移动策略，可能会导致突变率异常高，这可能是由于对训练集中不存在的突变体预测准确性有限所致。因此，MCMC模拟改为针对以下基于野生型的概率函数进行：
$$
\mathbf{p}_{\lambda}(x) \propto \mathbf{f}_{\lambda}(x) = \frac{1}{\mathcal{N}(x) > \sigma} e^{-\lambda \mathbf{g}(x) - \lambda \mathcal{N}(x)}
$$
其中$x$表示氨基酸序列，$\mathcal{N}(x)$表示$x$中与野生型序列匹配的残基数，$\lambda \in (-\infty, \infty)$表示对野生型序列的偏好程度，1表示0-1指示函数，$\omega$限制了模拟过程中允许的最大突变次数。$\omega$被设置为大于9个突变的突变体不被考虑。测试了几种$\beta$值，包括-5、-10、-15和-20，最终选择了-10以使Δlog KD的输出分布与野生型的适应度值重叠或高于野生型的适应度值。

为了提出单残基突变并接受相应的序列，采用了改进的Metropolis程序：
$$
\mathbf{q}(x'|x) = \mathbf{q}(n, r|x) \text{ if } \mathcal{N}(n, r) = \mathcal{N}_{WT}(n) \text{ else 0}
$$
其中$n$决定了要突变的残基的索引，$r$决定了突变后会变成哪种氨基酸，$L$表示可突变残基的数量。需要注意的是，这是一个关于$n$和$r$的联合概率质量函数；然而，由于提式的对称性，$n$并不显式出现在公式中。提出的突变根据以下Metropolis准则被接受：
$$
\mathbf{\eta}(x'|x) = \min\left[1, \frac{e^{-\lambda \mathbf{g}(x')}{e^{-\lambda \mathbf{g}(x)}} \cdot \alpha(x')
$$
应用方程5和6后，由于提式的不对称性，最终的稳态分布由方程4给出。也测试了不同的$\mathbf{C}$值，最终选择了0.95以平衡突变景观的多样性（图S5）。选定的$\mathbf{C}$值对应于$\lambda = 5.9$。最终的接受是通过从0到1的均匀分布中抽取一个随机变量并与方程6进行比较来决定的。如果被接受，突变体就被设为MCMC模拟的当前状态；如果被拒绝，模拟状态保持不变。观察到进行2500万步的MCMC模拟在多次运行中都能得到相同的结果，表明模拟过程已经收敛。

2025年1月31日，从NCBI Genbank（108）下载了总共9,017,742条完整的SARS-CoV-2受体结合域（RBD）基因组序列，筛选条件包括SARS-CoV-2分类单元（taxid:2697049）、人类宿主以及至少29,000个核苷酸的长度。使用Wuhan-Hu-1参考基因组（访问号NC_045512.2，位置22552–23155，NCBI RefSeq）中的603 bp RBD区域作为比对目标。

比对流程使用了多进程技术来并行处理257 GB的数据集，并将序列加载到一个有界队列中以管理内存限制。比对使用了BioPython的PairwiseAligner工具，参数设置为不匹配得分：-1，间隙开启惩罚：-2。在对齐过程中首先聚焦预期的RBD区域，并在必要时逐步扩展到全序列比对。模糊的核苷酸（N）被替换为其对应的武汉参考核苷酸。603个核苷酸的RBD编码序列根据标准遗传密码表转换为201个氨基酸。在RBD区域与参考序列差异超过20%的序列被排除（少于500条序列）。经过比对和过滤后的最终数据集包含8,255,834条RBD序列（占总数的91.55%），代表了27,980种独特的RBD序列变异，这些变异由1311种不同的突变组合而成。