从实验室进化实验推断蛋白质适应性景观

萨米尔·德科斯塔 ，艾米丽·海因兹 ，蔡斯·弗雷施林，宋惠彬 ，菲利普·罗梅罗

发布时间：1 年 2023 月 <> 日

抽象

定向实验室进化应用迭代的突变和选择轮次来探索蛋白质适应度景观，并提供有关蛋白质序列、结构和功能之间潜在关系的丰富信息。实验室进化数据由多代进化人群中采样的蛋白质序列组成，这种数据类型不适合既定的有监督和无监督机器学习方法。我们开发了一个统计学习框架，可以模拟进化过程，并可以从进化轨迹上的多个快照推断蛋白质适应度景观。我们将建模方法应用于二氢叶酸还原酶（DHFR）实验室进化数据，由此产生的景观参数捕获了DHFR结构和功能的重要方面。我们使用由此产生的模型来了解健身景观的结构，并发现了许多上位的例子，但总体上是一个全局峰值，可以从大多数起始序列进化上访问。最后，我们使用该模型对DHFR实验室进化轨迹进行计算机外推，并计算设计未来进化轮次的蛋白质。

作者摘要

实验室进化彻底改变了我们对蛋白质结构、功能和进化的理解，并产生了无数有用的蛋白质，在医学、生物催化和生物技术中得到了广泛的应用。这些实验通过迭代的突变和选择轮次探索蛋白质序列空间，并且可以提供穿越健身景观的人群的丰富数据。在本文中，我们提出了一个统计学习框架，该框架模拟了进化过程，并可以从实验室进化轨迹上的多个快照推断潜在蛋白质适应度景观的结构。我们生成了二氢叶酸还原酶（DHFR）实验室进化数据集，并应用我们的建模方法来推断景观参数。估计的参数可精确定位决定DHFR结构和功能的关键残基。我们使用所得模型来了解健身景观的局部和全局结构，并在蛋白质工程的计算机定向进化中执行。

数字

Fig 4图1图2图3Fig 4图1图2图3

引文： D'Costa S，Hinds EC，Freschlin CR，Song H，Romero PA （2023） 从实验室进化实验推断蛋白质适应度景观。公共科学图书馆计算生物学19（3）： e1010956. https://doi.org/10.1371/journal.pcbi.1010956

编辑 器： 李金燕， 澳大利亚悉尼科技大学

收到： 5年2022月16日;接受： 2023月 1， 2023;发表： <>月 <>， <>

版权所有： ? 2023 德科斯塔等人。这是一篇根据知识共享署名许可条款分发的开放获取文章，该许可允许在任何媒体上不受限制地使用、分发和复制，前提是注明原作者和来源。

数据可用性： 该实验生成的原始测序数据已上传到NCBI序列读取档案（SRA加入号PRJNA923701），标题为二氢叶酸还原酶（DHFR）的实验室进化。此外，用于重现这种计算方法的软件代码可在 Github： https://github.com/RomeroLab/dhfr_neutral_evolution 的开源许可证下获得，并在 https://doi.org/10.5281/zenodo.7622051 存档。

资金： 这项工作得到了授予PAR的美国国立卫生研究院（5R35GM119854）的支持。资助机构在研究设计、数据收集和分析、发表决定或手稿准备方面没有任何作用。

竞争利益： 作者声明没有相互竞争的经济利益。

介绍

从蛋白质序列到功能的映射形成了高维蛋白质适应度景观。了解这种景观对于理解和模拟自然进化、诊断遗传疾病以及设计用于生物技术、人类健康和化学的新蛋白质非常重要。这种景观是由高度复杂的蛋白质构象、动力学和生物物理/生化机制塑造的，并且是由天文数字般的大量可能的蛋白质序列定义的。虽然从物理角度建模序列函数映射具有挑战性，但可以利用统计学和机器学习的方法从稀疏采样的实验和进化数据中推断出底层景观[1-4]。

对蛋白质适应度格局进行建模的统计方法是围绕提供标记或未标记数据的两种常见数据类型构建的。标记的蛋白质数据由一组氨基酸序列以及每个序列如何映射到感兴趣的特定蛋白质特性（例如热稳定性、酶活性或结合亲和力）组成。这些序列功能数据通常使用蛋白质诱变文库和中通量或高通量测定法生成，以分配功能标记[5，6]。监督学习方法（如线性回归或更复杂的非线性模型）可以从标记的序列函数数据中学习，以推断序列与函数的映射[7-10]。未标记的蛋白质数据由取自基因组和宏基因组测序数据库的天然蛋白质序列组成。无监督学习方法可以从这些未标记的蛋白质数据中学习，以推断健身景观[11-13]。直接偶联分析（DCA）是一类重要的无监督学习方法，它从相关序列的多个序列比对中学习残基协同进化模式[14，15]。DCA方法已被用于预测蛋白质的三维结构[16-18]，模拟突变的影响[19]，以及设计新的蛋白质[20]。

定向实验室进化应用迭代轮次突变和选择来探索蛋白质适应性景观[7]。随着种群的进化，它会对序列空间的不同区域进行采样，并生成可用于理解适应度景观结构的进化轨迹。实验室进化数据包括从多个连续世代进化的种群中采样的蛋白质序列。这些数据自然不适合既定的有监督或无监督学习范式。以前的工作已经处理了类似于自然进化数据的实验室进化数据，并执行了无监督的DCA方法来推断景观参数[21，22]。虽然这些方法在确定三维结构中的接触残留物方面是有效的，但它们忽略了实验室进化数据的顺序性质，而是将来自多代的序列视为独立样品。[23] 演示了一种使用来自多轮深度突变扫描实验的数据推断上位健身景观的方法。此外，[24]使用遗传时间序列数据来推断上位健身景观。这些方法对随时间变化的数据如何相互连接进行建模。观察进化过程如何随着时间的推移展开，提供了有关健身景观结构的宝贵信息。

在这项工作中，我们开发了一个统计学习框架，从实验室进化数据中推断蛋白质适应性景观。我们使用群体遗传学原理来开发潜在进化过程的模型，并构建一个似然函数来估计多轮进化的景观参数。我们对二氢叶酸还原酶（DHFR）进行了15轮实验室进化，以生成一个庞大而多样化的数据集，该数据集由从多个连续世代中采样的序列组成。我们应用我们的学习方法来推断DHFR景观，发现模型参数捕获了DHFR功能的重要方面，并揭示了残基之间相互作用引起的景观上位。我们使用学习模型通过运行数千次进化模拟来了解健身景观的全局结构，并发现所有轨迹都收敛到相同的序列，这表明尽管有许多局部上位的例子，但总体最优值。最后，我们应用我们的模型从我们的实验DHFR进化停止的地方开始，并在计算机中继续进化过程。该程序用于推断进化轨迹并设计超出训练数据的新功能DHFR。

结果

实验室进化探索二氢叶酸还原酶的适应性前景

实验室进化应用迭代轮次的突变和选择来探索蛋白质适应性景观。我们对小鼠二氢叶酸还原酶（mDHFR）进行了实验室进化实验，以搜索编码DHFR活性的不同序列的适应性景观（图1）。DHFR将二氢叶酸还原为四氢叶酸，在嘌呤生物合成和细胞生长中起着至关重要的作用。

缩略图 下载：

.PPT幻灯片

.PNG大图

.TIFF原始图像

图1.

实验性DHFR实验室进化 （a）实验室进化结合了随机诱变和功能选择的迭代轮次，以进化分子群。我们在DHFR上进行了15轮进化，并对第1-5代和第15代的种群进行了测序，以获得沿进化轨迹的快照。（b）实验进化种群的序列统计显示，与野生型DHFR的距离增加，种群内序列之间的距离增加，表明扩散样扩散到蛋白质序列空间。（c） 在自然DHFR序列的背景下可视化第15轮种群。实验室进化实验探索了DHFR序列空间的一小部分。使用在自然DHFR序列上训练的变分自动编码器（VAE）的二维潜在空间对序列进行可视化（VAE的详细信息可在S2文本中找到）。（d）来自第15轮群体的序列的突变统计显示，活性位点残基的突变率较低。突变N20D在第15轮种群中超过了野生型天冬酰胺残基。

https://doi.org/10.1371/journal.pcbi.1010956.g001

我们采用了一种常用的选择策略，该策略应用抗生素甲氧苄啶来抑制大肠杆菌的天然DHFR，并使细胞依赖于异源表达的mDHFR，即耐药的甲氧苄啶。我们使用容易出错的PCR诱变mDHFR，每个基因的目标突变率为四个核苷酸替代。然后，我们将该文库转化为大肠杆菌，并进行了一项选择，以鉴定能够在甲氧苄啶存在下支持大肠杆菌生长的DHFR变体。在此之后，我们从所有幸存的变异中提取质粒DNA，用容易出错的PCR对这些变异进行再变向，重复选择过程，总共进行了15轮实验室进化。对于每一轮，我们跟踪转化体的总数和功能性DHFR变体的比例;这两个数字的乘积给出了人口规模的估计值（S1表）。在15轮进化中，该种群的平均规模为300，000个DHFR变体，并且从未少于40，000个。与像[25]这样的实验在每一轮进化中只选择最适者克隆并尽可能增加选择压力相反，我们的低严格选择和大种群规模创造了一个中立的进化过程，产生不同的序列，保持野生型DHFR活性。

我们在定向进化第1-5轮和第15轮进行了Illumina DNA测序，以获得进化轨迹的样本。从这些测序数据中，我们看到进化过程产生了与野生型mDHFR具有不同汉明距离的序列分布，并且进化的种群越来越多地偏离起始序列（图1b和S1图）。

第 15 轮种群平均有 11.9 个氨基酸取代，相当于每轮累积 0.79 个氨基酸取代。我们还观察到种群的成对汉明距离在进化过程中线性增加。与野生型的平均距离和平均成对距离都随着每一轮线性增加的事实表明，大多数进化轨迹都在探索景观上的独立方向。我们还可视化了定向进化实验中产生的DHFR序列如何适应自然进化产生的更大蛋白质家族（图1c）。我们的定向进化实验开始捕捉到与一些自然序列类似的变异，但只探索了自然界跨越的序列空间的一小部分。

我们进一步分析了第15轮序列，以了解突变如何在一级序列中分布。突变通常分布在整个蛋白质序列中，但在蛋白质核心和活性位点残基中观察到的突变率较低。我们还观察到N20D突变在第15轮人群中超过了野生型天冬酰胺残基，可能表明该突变的正选择。残基20位于mDHFR的活性位点区域，与核苷酸磷酸盐结合。

从连续几轮实验进化中学习的统计框架

定向进化提供了在蛋白质适应性景观上进化的种群样本。我们开发了一个统计框架，从这些实验进化轨迹中推断潜在的景观结构。我们假设来自连续进化轮的数据比来自单轮的数据或来自独立考虑的多轮数据提供了固有更多的信息。观察进化过程如何在多轮中展开，使我们能够做出更有力的推断和推断行为。对于简化的说明性示例，如果我们观察氨基酸的频率在每一轮中如何变化，我们可以推断这些轨迹来估计进化过程将在哪里收敛。我们建立了实验室进化过程的生成模型，使用广义Potts模型对健身景观进行参数化，并从多轮进化的测序数据中推断景观参数。

我们将实验室进化的动力学建模为马尔可夫链过程，其中序列根据其突变可及性和相对适应性过渡到其他序列。我们对几轮进化之间的序列转换机制做出了一些假设。首先，我们假设突变过程在每个DNA位置独立发生，并且每个位置的突变概率从实验中已知。其次，我们假设转移机制是时间均匀的，即适应度值和突变概率在轮次之间都没有变化。第三，我们假设实验转化体的数量足够大，因此给定轮次的序列分布仅取决于它们的相对适应度水平。最后，我们假设简化的马尔可夫链动力学，它假设直系后代之间的局部竞争，非常接近真实动力学（参见S1文本中无限种群动力学的马尔可夫链近似）。

我们使用（广义）Potts模型参数化健身景观，该模型描述了所有氨基酸残基及其成对相互作用如何促进健身。Potts模型已被广泛用于恢复蛋白质残基之间的相互作用图，并且强相互作用已被证明对应于蛋白质3D结构中的远程接触[15，17]。定向进化数据由从多轮进化中采样的序列组成，每轮观测到的序列是马尔可夫模型概率分布的随机实现。我们可以通过最大化统计似然函数来估计 Potts 模型参数，以获得给定观测到的进化数据的最可能的模型。我们的进化模型和参数估计方法的详细信息在方法部分给出。

学习的景观参数捕获蛋白质结构和功能

我们应用上面开发的统计学习框架，从我们的实验实验室进化数据中推断出DHFR的健身景观。我们估计了进化第1-5轮和第15轮的Potts模型的典型参数，并使用这些参数来获取有关蛋白质结构和功能的信息。

学习的Potts模型揭示了单个氨基酸取代如何影响野生型DHFR的活性（图2a）。这种学习的突变图清楚地突出了酶的关键催化残基的重要性，并显示了具有相似理化特性的预期突变模式残基。我们将该突变图谱与使用DCA [19]和贝叶斯VAE使用EVE [26]从天然DHFR序列中学习的突变图进行比较（S2图）。与DCA相比，EVE方法旨在捕获更高阶的相互作用。我们推断的突变图与DCA的突变图谱之间的基于Spearman秩的相关性为0.54，与EVE的相关性为0.62，表明从进化实验中学到的健身景观与从自然序列中学到的相似但不相同。G18A氨基酸取代在推断模型中具有最大的有益作用，并且位于与NADPH相互作用的环中。我们将平均突变效应幅度映射到三维DHFR结构上，以了解学习的参数与结构的关系（图2b）。我们观察到蛋白质核心的突变往往对活性有更大的影响，可能是因为这些位点的突变破坏了酶的三维结构。

缩略图 下载：

.PPT幻灯片

.PNG大图

.TIFF原始图像

图2. 模型参数与DHFR的结构和功能有关。

（a） 模型在186个DHFR序列位置的预测突变效应的热图。野生型氨基酸呈白色。（b）映射到DHFR结构上的每个位点的平均突变效应（PDB ID：3K47）。最大的突变往往发生在蛋白质核心和底物结合位点。（c） DHFR中所有站点对之间的相互作用强度。相互作用强度计算为Potts模型相互作用系数上每对位点所有氨基酸组合的Frobenius范数。（d） 在DHFR结构上绘制的前十个长程相互作用（PDB ID：3K47）。其中许多相互作用是通过基质的相互作用发生的。（e） 在接触图上绘制的残基之间的顶部 L/2 （93） 相互作用，显示重原子接近 5 ? 和 8 ? 的残基。 （f） 比较在 R15 数据上训练的伪似然 DCA 模型的接触预测，重新加权以解释进化偏差的 R15 数据，或使用我们在进化轨迹上的方法。水平线表示随机机会，垂直线在常用的 L/2 阈值处绘制。

https://doi.org/10.1371/journal.pcbi.1010956.g002

Potts模型也可用于了解景观上位和残基之间的相互作用。我们计算了野生型所有单突变体和双突变体的适应度评分，并在大约6万个可能的氨基酸和位置对中确定了景观中倒易征上位的47个例子。在推断的景观中，一个这样的例子发生在K99A和R24695E之间，其中单个突变单独是有益的，但是当组合在一起时会导致适应性下降。我们在景观中观察到了<>个标志上位的例子，其中一个突变在存在另一个突变时具有相反的效果。我们通过计算弗罗贝尼乌斯范数 F 来计算残基-残基相互作用评分ij一对残基之间的所有相互作用参数之间（图2c）。推断的残基相互作用从酶的三维结构以及酶活性位点中的功能相互作用中捕获了许多接触。顶部相互作用在残基R71和G117之间，它们在3D结构中不直接相互作用，而是形成核苷酸结合口袋的相对末端（图2d）。十大残基-残基相互作用中有七种涉及这些关键位点之一。前 20 个交互在 S2 表中提供。

Potts模型相互作用评分可用于识别在三维蛋白质结构中接触的残基对。具有最高 L/2 （93） Frobenius 分数的残基对对应于距离小于 10 ? 的 5 个结构接触和距离小于 21 ? 的 8 个接触（图 2e）。我们将模型的接触预测性能与从进化数据推断接触的既定方法进行了比较，包括标准的DCA建模程序[19]以及改进的DCA模型，该模型根据序列与野生型的距离对序列进行加权，以解释进化过程[22]（图2f和S3表).当在第 15 轮数据上进行训练时，在预测距离小于 2 ? 的顶级 L/79 （2） 远程接触时，这两种方法都无法超过随机机会预期 （93.5%）。相比之下，我们的模型考虑了进化轨迹，能够正确识别出前 L/10 接触中的 2 个（10.7%），这远高于随机机会。我们的模型还通过恢复与重原子5?8 ?和小于8 ?的更多结构接触来优于其他两种方法（图2f和S3表）。

在计算机进化模拟中，绘制健身景观的全局结构并推断进化轨迹

我们的统计方法从实验实验室的进化轨迹推断出潜在的健身景观。该模型可用于在计算机模拟中运行，以了解景观、进化过程和设计新蛋白质。

我们使用我们的模型来了解健身景观的全球结构和适应性进化步行的收敛。我们对所有 1600 种氨基酸的 186 个随机氨基酸序列（长度 20）进行了均匀采样，以获得广泛的景观采样。对于这些序列中的每一个，我们通过评估所有单个突变体，选择最合适的变体并重复此过程直到达到局部适应峰值来进行适应性步行。我们发现每一个适应性进化轨迹都收敛到同一个适应度峰值，这个适应度峰值的序列是来自野生型DHFR的79个氨基酸替代。从景观的不同区域开始的适应性步行汇聚到同一山峰，这一事实意味着富士山式的健身景观几乎没有局部最佳。

我们还使用学习的模型继续实验DHFR进化过程，并将进化轨迹外推到后代（图3b）。我们从实验室进化实验最后一轮中观察到的最常见的氨基酸序列开始模拟。该起始序列是来自野生型DHFR的10个氨基酸替换。我们通过评估所有单突变体和双突变体，选择最合适的变体，并重复此过程，直到没有进一步的上坡步骤来进行适应性行走。模拟的进化轨迹继续远离野生型DHFR和第15轮起始序列，并收敛到上述全球景观搜索中发现的相同适应度峰值。

缩略图 下载：

.PPT幻灯片

.PNG大图

.TIFF原始图像

图3. 进化模拟以了解景观结构。

（a）我们均匀采样长度为186的随机氨基酸序列，并进行适应性行走，直到达到适应度峰值。所有适应性步行都汇聚到同一个健身高峰。显示了11个具有代表性的进化轨迹。（b）我们通过继续计算机的进化过程来继续我们的实验室DHFR进化。我们沿着进化轨迹对几个DHFR序列进行了实验测试，发现许多是无活性酶。

https://doi.org/10.1371/journal.pcbi.1010956.g003

我们想测试定向进化轨迹是否可以在计算机中外推，作为设计新蛋白质的方法。我们沿从第 15 轮序列到自适应适应度峰值的轨迹实验表征了十个 DHFR 序列（图 3b）。前五个序列沿着这个轨迹连续拾取，接下来的五个序列通过跳过每隔一个序列来拾取。我们没有对整个轨迹进行采样。我们发现第15轮序列是一种活性DHFR酶，第15轮序列的双突变体也是活跃的，但除此之外的所有序列都是无活性的，在甲氧苄啶存在下无法补充大肠杆菌的生长（图3b和S3图）。

讨论

定向实验室进化是探索蛋白质适应性前景的有力方法。定向进化生成由多代采样序列组成的数据，并提供有关潜在健身景观结构的宝贵线索。定向进化数据自然不适合既定的有监督或无监督学习方法，因为它们不考虑进化数据生成过程。在这项工作中，我们开发了一个统计学习框架，从实验室进化数据中推断蛋白质适应性景观。我们建立了定向进化过程的生成模型，使用广义Potts模型对健身景观进行了参数化，并从多轮进化的测序数据中推断了景观参数。我们将这种学习方法应用于大型且多样化的二氢叶酸还原酶（DHFR）定向进化数据集。推断的景观模型揭示了许多由残基之间的相互作用引起的上位的例子，但总体全球适应度峰值可以从大多数起始序列进化上获得。最后，我们探索了景观模型推断蛋白质工程进化轨迹的潜力。

为了计算的可追溯性，我们的实验室进化模型做了两个重要的近似。了解这些近似值何时有效或无效有助于确定我们方法的范围和局限性。第一个近似假设进化中的种群具有无限的种群规模，而在实践中，所有实验室进化都具有由实验约束决定的有限种群规模。我们的实际DHFR实验种群规模约为105?106而这个数字很大程度上取决于文库转化为大肠杆菌的效率。人口规模越大，这种近似值会变得更加准确。

第二个近似假设种群根据简化的马尔可夫链动力学（S1 Text）进化，其中竞争只发生在共享直接祖先的序列之间。在我们的实际实验中，所有序列都直接参与生长选择，无论它们的谱系如何（S4图）。当竞争序列在实验的每个阶段具有相似的适应度时，马尔可夫链动力学近似是有效的，如果初始野生型序列附近的适应度景观是平坦的（即，大多数选定的突变是中性的），并且进化实验以较低的突变率进行。在低突变率下，实验室进化实验仅探索野生型序列的局部邻域，并且由于中性景观，所有选定的突变序列将具有相似的适应性。此外，马尔可夫链动力学在前几轮进化中是一个更好的近似，而所有演化序列仍然接近初始野生型序列。

一旦我们指定了进化动力学，我们就通过近似时刻匹配来推断健身景观参数。轮 r 和 r + 1 之间的一阶和二阶矩的近似递归公式是基于感兴趣的变量与其余变量之间的弱依赖性假设得出的，这允许我们仅根据蛋白质残基图中的局部信息近似 r + 1 轮的边际。这种近似只有在残差图分解时才是精确的，但我们希望这为大多数节点和/或节点对提供了合理的近似，其中它们与其他变量之间的交互强度是有限的。在我们的推理过程中，我们有三个主要的正则化超参数。前两个影响主效应的大小和成对效应参数。这些设置方式与DCA模型相同[17]。第三个超参数是正则化参数，它从平均参数推断规范参数。当它设置为比估计协方差矩阵中最小特征值的负数大几倍时，我们观察到最佳性能。在这三个参数中，结果对最后一个参数最敏感。

我们推断的DHFR景观参数能够识别三维蛋白质结构中的接触残基以及通过酶活性位点间接相互作用的众多功能偶联残基。我们的轨迹学习方法在常用的L/2阈值下，对于重原子接近5 ?的残基，显示出比既定的接触预测方法有显着改进。这种改进是由于对进化轨迹进行建模以及包括前几轮进化的测序数据。推断的几个顶级相互作用（S2表）位于3D结构中不接近的残基中，而是位于与底物或辅因子相互作用的功能区域中。这些相互作用在接触预测中被标记为假阳性，但它们可能是由配体相互作用或辅因子重新定位引起的上位相互作用[27]。S17表中前20个相互作用的残基对中有2个涉及色氨酸残基。虽然这代表了仅与三个残基（25W，58W和114W）的重复相互作用，但这种过度表示可能是由于模型假设中关于容易出错的PCR转移概率的偏差。然而，已知色氨酸残基在大肠杆菌DHFR的结构和稳定性中起重要作用[28]，特别是残基25在天然DHFR序列中高度保守，并且在蛋白质的功能中起作用[29]。所以另一种解释是，这些色氨酸残基出现在S1表中是因为它们的重要性。

以前的工作已经证明了通过将DCA应用于最后一代实验室进化实验来识别残留接触的能力[21，22]。Stiffler等人对DCA应用了重新称重程序以解释进化偏差，并且能够使顶级L / 50接触的接触恢复率超过2%。以前使用TEM1 β-内酰胺酶[30]的实验室实验无法检测局部上位相互作用，但是更大和最近的实验[21]能够识别上位相互作用并使用DCA检测一些接触者。值得注意的是，这些DCA方法无法从我们的第15轮进化数据（S5图）中可靠地检测DHFR接触。我们将标准DCA和[22]的加权DCA方法应用于我们的第15轮DHFR数据，当预测2D结构中接近5 ?的顶部L / 3接触时，正确标记的接触的百分比低于随机机会预期。有许多因素可能导致这些数据集之间接触恢复的差异，包括群体序列多样性、测序的变异数量或进化突变率。Stiffler等人提供了氨基糖苷类乙酰转移酶（AAC6）的下采样数据集分析[22]，因此我们可以直接与我们的mDHFR数据进行比较。在 105序列，它们的平均成对距离为10.9%，并识别近40%的接触者。相比之下，我们的mDHFR数据为105序列中，我们的平均成对距离为11.1%，但只能恢复2.15%的接触。基于这些发现，序列多样性和数据集大小似乎没有导致接触恢复的差异。这两个实验在突变率方面确实存在显着差异，其中AAC6每轮0.8%，而DHFR每轮进化0.4%。我们假设这种不同的突变率导致了不同的种群结构，尽管具有相似的大小和多样性水平。即使对于固定的中性景观，遵循具有不同突变率的准物种动态的有限种群也会进化出不同程度的稳健性[31]，这种稳健性与上位相互作用有关[32]。使用简化的马尔可夫链动力学对这些类型的实验室进化实验的模拟由[33]提供，这些模拟没有将突变率视为重要参数。然而，由于实验室进化实验可能遵循更复杂的准物种动力学，因此突变率可能起着更重要的作用。

我们使用推断的模型参数来探索上位和健身景观的全局结构。我们发现了许多加性突变的例子，大约0.4%的突变通过符号上位相互作用，只有34个互惠符号上位的例子。这种上位相互作用的频率与其他研究一致[35，35]。尽管存在景观上位，但我们发现估计的景观结构具有整体全球适应度峰值，可以通过大多数起始序列的适应性步行来获得。在其他蛋白质中也观察到了类似的健身景观特征[36，<>]，但是，由于模型在其训练数据之外进行泛化的能力存在各种近似和限制，我们警告不要将模型的估计景观特征解释为准确存在于实际蛋白质健身景观中。

我们使用推断的Potts模型进行计算机定向进化实验，以设计新的，以前未观察到的DHFR变体。我们发现该模型可以设计接近训练数据制度的功能性DHFR，但进一步的进化外推导致了非功能性酶。这一结果与其他基于机器学习的蛋白质工程研究一致，这些研究表明模型准确性降低，同时从训练数据推断[10，37]。这一结果也表明，在79个突变处预测一个适应度峰值是不可靠的。模型不准确可能是由于数据不足或质量低下、我们在进化模型中所做的近似或参数估计的计算挑战的结果。更准确的进化模型将考虑有限的种群规模和一代中所有序列之间的竞争，这被称为有限准物种模型[31，38]。提高模型准确性的另一种可能性是使用没有交互作用的更简单的一阶模型，以便我们需要估计更少的参数。提高蛋白质工程可靠性的另一种方法是运行多个独立的进化模拟并测试一组不同的设计。

我们的统计景观推断方法自然地补充了连续定向进化的最新进展[39-41]。这些实验方法在连续喂养的生物反应器中结合了种群水平的诱变和选择，以进化种群而无需离散的突变/选择步骤。可以通过下一代DNA测序对种群进行采样和分析，以观察种群如何随时间变化并穿越健身景观。我们的学习方法可以从这些连续的进化数据中推断出景观，以了解蛋白质的结构，功能和进化。

蛋白质序列、结构和功能之间的关系涉及数千种精致的分子相互作用，这些相互作用在空间和时间上动态耦合。机器学习正在彻底改变我们对这些关系的理解，它以超出人类理解的规模和分辨率剖析蛋白质的复杂内部运作。数据驱动的蛋白质科学的未来进展将提高我们理解自然进化过程，预测遗传疾病以及设计新蛋白质以在生物技术中广泛应用的能力。

材料和方法

统计方法概述

我们将实验室进化的动力学建模为马尔可夫链过程，其中序列根据其突变可及性和相对适应性过渡到其他序列。我们考虑长度为 L Ω的所有可能密码子序列的集合，其中每个序列用 x ? （x ） 表示1， ...， xL).每个 x我对应于编码第 i个残基位置的密码子，每个密码子来自排除终止密码子的 61 个密码子集合。我们让 Π（x） 表示每个序列 x ∈ Ω 的适应度，即序列在单位时间内对自身产生的拷贝数，π（x） 表示基于相对适应度的相应患病率，定义为 π（x） = Π（x）/∑你∈Ωπ（u）。

在（S1文本）的假设下，我们有一个转移概率，由下式给出

(1)

其中 g（x → y） 表示在没有选择的情况下突变为序列 y 的概率序列 x。每个序列 x(n， r)从舍入 R 是边际概率 p 的随机样本(r)在第 R步转换之后。

我们使用（广义）Potts模型π参数化健身景观，该模型描述了所有氨基酸残基及其成对相互作用如何促进适应性。我们将适应度水平参数化π（x） = πθ(x） 序列 x） 的序列 x，具有 Potts 模型，在氨基酸上设置了规范参数，其中

(2)

AC（?）是从密码子集合到集合到氨基酸的映射，Z（θ）是归一化常数，使得概率之和为1。此模型由 q一个l 主效应参数（h我），其中 和（耦合）交互作用效应参数 （eij).

定向进化数据由从多轮进化中采样的序列组成。这些观察到的序列 x(n， r)， n = 1， 2， ...， nr在每一轮 r 是概率分布 p 的随机实现(r)(?;我们可以通过最大化以下对数似然来估计 Potts 模型参数 θ。

(3)

其中表示存在测序数据的实验轮数集，p（x → y;（1） 中的 θ） 由波茨模型下的马尔可夫链动力学给出。

确定最大化该对数似然函数 （3） 的模型参数 θ 具有挑战性。首先，即使πθ形成马尔可夫随机场（MRF），分布为p(r)不再对与适应度分布π关联的图形进行分解θ.特别是，有条件的独立关系，它πθ不要持有 p(r)， r = 1， 2， ...，这禁止在马尔可夫随机场设置中应用参数推理技术。此外，状态空间的大维数排除了对马尔可夫链过程的任何精确跟踪。例如，转移矩阵的维数为 |Ω|× |Ω|，即使计算给定 θ 的转移矩阵中的单个元素在计算上也是不可行的，因为 （1） 中分母中的和难以处理（如 |Ω| = 61186).

我们使用近似矩匹配方法来克服这些计算挑战。特别是，我们首先推导出了一阶和二阶边际之间的近似关系，在每个回合r和μ下边际πθ（第二节 边缘人之间的近似关系）。然后，我们寻找哪些与经验和预期的一阶和二阶边际相匹配，并且也是局部兼容的（推断平均参数）。最后，我们使用平均场DCA [1]，使用估计的边际作为输入来获得Potts模型（2）的参数θ的估计值。（请参阅 S15 文本中的规范参数估计）。

我们使用规范参数θ的估计值来获取有关蛋白质结构和功能的信息。每个参数 eij(a，b）表示残基I处的氨基酸A和残基J处的氨基酸B之间的相互作用。我们使用规范参数e以与DCA方法相同的方式计算残基i和j之间的相互作用分数[17]。参数集 e 过度参数化，因此我们首先转换为零和规范，然后计算弗罗贝尼乌斯范数

弗罗贝尼乌斯规范 Fij是残基对 i 和 j 之间的交互作用得分。我们只关注长程相互作用（即氨基酸序列中相差超过5个位置的残基之间）。

图4给出了一个概述我们方法的示意图。虽然我们以这种方式组合在一起来推断健身景观的框架没有出现在文献中，但有几个组成部分来自以前的工作。通过 Potts 模型参数化 （2），我们可以将 （1） 给出的简化动力学视为类似于 [33] 中给出的进化动力学。我们框架中的概括是，我们允许多个同时突变的可能性，其中[33]仅使用单个突变方案来有效地实施模拟。此外，估计典型 Potts 模型参数上的组件遵循相当标准的平均场 DCA [15] 方法。处理矩匹配（边际之间的近似关系）和（秒推断均值参数）的组件应被视为新的贡献。

thumbnail 下载：

.PPT幻灯片

.PNG大图

.TIFF原始图像

图4. 用于推断健身景观参数的算法概述。

https://doi.org/10.1371/journal.pcbi.1010956.g004

重现此方法的软件代码可在 github.com/RomeroLab/dhfr_neutral_evolution 获得，并在DOI 10.5281/zenodo.7622051存档。

边缘之间的近似关系

我们首先用一阶和二阶边际μ π的函数近似一阶和二阶边际θ.获得所有可用回合的近似矩向量后，μ π的平均向量θ然后通过将近似矩向量与经验计数匹配来估计。

我们让 X(r)是概率分布 p 之后的随机变量(r).P（X(r) = x） ? p(r)(对于任何密码子，首先我们将回合 r 处的一阶和二阶边际定义为 ，其中 δ（?， ?） 是克罗内克三角洲函数，使得 δ（a， b） = 1 如果 a = b，否则为 0。我们让矩向量作为这些项的集合

同样，我们还定义了对应于健身景观的一阶和二阶边际πθ作为 、 以及π的矩向量μθ如

如果我们在 r 轮有测序数据，我们可以通过对数据求和来计算在该轮中观察到的单个和成对密码子频率。

我们可以将这些项放在一个频率矢量f中(r)

当可以根据数据计算时，它为我们提供了矩矢量的估计μ(r)对于回合 R。

从不观察健身景观本身，而只观察健身景观上的实验动态。因此，与轮次的矩向量不同，我们有排序数据并计算 f(r)作为μ的估计(r)，我们没有对健身景观μ的当前向量的估计。为了解决这个问题，对于每个节点i或一对节点（i，j），我们根据矩向量μ近似密码子c在i处的流行率和在i， j处r + 1处的密码子对（c，d）的流行率(r)的 R 轮和健身景观的矩向量μ。特别是，假设一组感兴趣的节点之间的依赖关系X。S和其余变量 XSc对于近似，我们在 S1 文本中导出以下递归关系：

(4)

我们注意到，一旦指定，所有平均向量随后通过 （4） 中的递归关系指定为 μ 的函数。由于在实验开始时，我们只有野生型序列的副本（用w ∈ Ω表示），因此它的频率为1，所有其他序列的频率为0。因此，如果 w我 = C 和 Wj = d 在所有其他情况下为 0。

现在，我们通过最小化所有位置i，j ∈ [L]，i≠ j之间的总对数损失来最小化预期的一阶和二阶频率（作为μ的函数）与观察到的一阶和二阶频率之间的 总对数损失来估计μ。在此过程中，我们对μ进行了重新参数化，以便我们可以在氨基酸水平参数上优化目标函数。我们还强制执行局部一致性条件，并正则化参数以防止优化期间过度拟合（秒推断均值参数）。

推断均值参数

在这里，我们通过最小化预期和观察到的一阶和二阶频率之间的总对数损耗来描述近似的矩匹配。换句话说，我们想解决以下目标函数：

(5)

其中 和 是足够的π统计数据θ，和 换句话说，该集合对应于全局一致的平均向量的集合，即所有一阶和成对边际概率，可以通过 {0， 1} 上的某个分布来实现d其中 d 是充分统计量φ的维度，集合对应于平均向量集，使得平均值仅取决于输入序列的氨基酸值。

首先，为了优化氨基酸水平参数的目标，我们将 Potts 模型μ Ω的平均参数重新参数化如下：

(6)

我们还定义 ν集(b， a;γ） = νij(甲、乙;γ） 对于 j > i。

虽然该集合可以用有限数量的线性不等式来表征，但线性不等式的数量会根据维度d而快速增长，并且通常，除非维度d很小，否则即使是线性目标也很难优化[42]。我们继续考虑通过强制规范化条件来放宽优化问题 （5）

(7)

和局部一致性条件

(8)

注 （7） 通过重新参数化 （6） 满足。我们添加以下带有拉格朗日乘数ρ的惩罚项，以促进局部一致性条件（8)

最后，为了处理γ的高维性，我们添加l2-正则化项和超参数 λ主要和 λ国际

总之，我们解决了以下优化问题：

(9)

我们将优化γ参数初始化为我们拥有测序数据的最后一轮成对频率的对数。将一个小的伪计数添加到频率中，以便为缺失的频率定义此初始化过程。我们使用自动微分计算目标的一阶导数。此外，我们使用梯度下降优化器和早期停止规则。

优化后，我们得到单个和成对估计，并可以使用等式（6）计算平均参数的估计值。

统计框架参数

统计推理方法是使用PyTorch在Python中实现的[43]。优化方法使用学习率为 0.03 的 Adam 优化器，其他学习参数设置为默认值并训练 300 步。（9）中的正则化超参数设置为λ主要= 10?3, λ国际= 10?4和 λ注册= 50。确保主效应参数边缘化为成对参数的惩罚项设置为 ρ = 105.

该实验中的确切突变分布尚不清楚，因为所有轮次都是在生长选择后测序的，但是，统计方法似乎对突变分布的选择相当稳健，因为结果看起来与使用类似方案完成的其他实验的突变偏差分布相似（结果未显示）。我们从[2]，[44]中表45的Taq DNA聚合酶列中选取突变偏向分布来模拟诱变，然后缩放以匹配平均每轮4个DNA突变。统计方法中使用的最终突变分布在S4表中给出。

mDHFR实验室演变

我们的选择菌株由克隆到pET22b质粒中的鼠二氢叶酸还原酶（mDHFR）基因组成，并转化为大肠杆菌BL21（DE3）。我们使用反应的MnCl进行了容易出错的PCR。2浓度以调节Taq DNA聚合酶的突变率[46]。我们确定最终浓度为 200 μM MnCl2每个基因产生3.25±0.74个氨基酸取代。我们进行了15个容易出错的PCR循环，用DpnI处理反应过夜以去除模板，用DNA离心柱纯化PCR产物（Zymo Research），使用环状聚合酶延伸克隆（CPEC）将插入片段克隆回pET-22b[47]，使用DNA离心柱纯化CPEC反应（Zymo Research），并将CPEC反应转化为电感受态BL21（DE3）细胞（Lucigen）。对转化进行几次稀释以确定总文库大小，其范围为 105?106菌落形成单位 （CFU）。转化的其余部分用作含有 100 μg/mL 羧苄青霉素、100 μ M IPTG 和 500 μg/mL 甲氧苄啶的 5 mL LB 竞争性生长选择的输入。我们让这些选择培养物在16°C下振荡生长37小时。 使用大约20 ODU的过夜培养物通过miniprep收获质粒DNA。然后将选定的质粒DNA群用作下一轮容易出错的PCR的模板。部分选择后培养物也作为15%甘油储备液存档，并储存在-80°C。

我们通过从每个转化板中挑选菌落到含有 LB 肉汤、96 μg/mL 羧苄西林、100 μM IPTG 和 500 μg/mL 甲氧苄啶的 5 孔板的单个孔中，确定了每轮进化中功能变异的比例。我们将这些板培养物在16°C下振荡孵育37小时，测量每个孔的OD 600，如果DHFR变体的OD 600大于0.5，则将其归类为功能性，否则它们被认为是非功能性的。

进化人群的DNA测序。

我们在几轮实验室进化中进行了下一代DNA测序。我们使用Illumina测序分析了第1-5轮，使用Pacific Biosciences测序分析了第15轮。对于 Illumina 文库，我们使用 NdeI 和 SacI 限制性内切酶从质粒中去除 DHFR 插入片段，将 Illumina 衔接符序列连接到该插入片段，并将样品提交给威斯康星大学麦迪逊分校生物技术中心 DNA 测序核心，使用 2x300 v3 试剂盒在 Illumina MiSeq 仪器上运行。每个样本有2-5百万次读取。对于PacBio测序，我们用NdeI / SacI移除了DHFR插入片段，并将样品提交给威斯康星大学麦迪逊分校生物技术中心DNA测序核心，以分析其Pacific Biosciences Sequel仪器。PacBio 运行返回了 10 多个5读。原始测序数据可在NCBI序列读取档案（SRA加入ID PRJNA923701）中找到。

序列预处理。

对于第 1-5 轮，Illumina 测序数据采用类似于 [22] 的步骤进行处理。正向和反向读数首先使用FLASH程序拼接在一起[48]。在第一个过滤步骤中，仅保留每个碱基的最小长度为500且最低质量得分为15的序列。在第二步中，保留复合质量得分至少为10的序列，这意味着没有读取错误的概率为90%。对于第 15 轮，进行了质量过滤以保持最小长度为 564、最低质量得分为 0.99 的序列。然后使用领结2将其余序列与参考序列对齐[49]。

实验测试进化设计的DHFR。

我们使用推断模型设计了十个DHFR变体来模拟进化过程。编码这些变异的基因由Twist Bioscience合成并克隆到pET21（+）质粒中。我们将质粒转化为大肠杆菌BL21（DE3）并进行生长测量以评估DHFR变体的活性。我们通过首先接种含有 5 μ g/mL 羧苄青霉素的 100 ml LB 起始培养物并在 37°C 下振荡生长过夜来进行生长测定。 然后，我们将该起始培养物稀释 100 倍到含有 100 μ g/mL 羧苄青霉素、500 μ MIPTG 和 5 μg/mL 甲氧苄啶的 LB 培养物中，并通过在 600°C 下 30.16 小时的孵育期内以 5 分钟的间隔测量 OD 37 来监测生长。 这些测量一式三份进行。在这些条件下，非活性DHFR变体显示无生长，而活性变体显示标准生长曲线。

支持信息

与野生型的平均汉明距离和每轮序列之间的平均成对汉明距离。

显示 1/11： pcbi.1010956.s001.pdf

跳到无花果共享导航

补充信息：从实验室推断蛋白质适应性景观进化实验补充数字图 S1.与野生型的平均汉明距离和平均成对汉明距离每轮序列显示大致线性进展1

1 / 11

下载

无花果分享

S1 图 与野生型的平均汉明距离和每轮序列之间的平均成对汉明距离。

https://doi.org/10.1371/journal.pcbi.1010956.s001

（英文）

S2 图 对自然DHFR序列的突变效应预测。

https://doi.org/10.1371/journal.pcbi.1010956.s002

（英文）

S3 图 使用推断模型设计的序列的生长曲线。

https://doi.org/10.1371/journal.pcbi.1010956.s003

（英文）

S4 图 理想化实验与真实实验的示意图。

https://doi.org/10.1371/journal.pcbi.1010956.s004

（英文）

S5 图 DCA方法在第15轮DHFR进化数据上进行了训练。

https://doi.org/10.1371/journal.pcbi.1010956.s005

（英文）

S1 表。 菌落采矿单位和分数函数的估计。

https://doi.org/10.1371/journal.pcbi.1010956.s006

（英文）

S2 表。 残基对之间的前20个远程相互作用得分。

https://doi.org/10.1371/journal.pcbi.1010956.s007

（英文）

S3 表。 通过不同方法恢复的远距离接触者的阳性预测值（PPV）。

https://doi.org/10.1371/journal.pcbi.1010956.s008

（英文）

S4 表。 容易出错的PCR突变偏倚。

https://doi.org/10.1371/journal.pcbi.1010956.s009

（英文）

S1 文本。 补充方法和数学细节。

https://doi.org/10.1371/journal.pcbi.1010956.s010

（英文）

S2 文本。 VAE架构和超参数。

https://doi.org/10.1371/journal.pcbi.1010956.s011

（英文）

引用

1.杨国庆， 吴志， 阿诺德.用于蛋白质工程的机器学习引导的定向进化。纳特方法。2019;16(8):687–694.密码：31308553

查看文章PubMed/NCBI谷歌学术搜索

2.Ferguson AL，Ranganathan R.大分子科学100周年 观点：数据驱动的蛋白质设计。宏研究学报 2021;10（3）：327–340.密码：35549066

查看文章PubMed/NCBI谷歌学术搜索

3.Bepler T，Berger B.学习蛋白质语言：进化，结构和功能。细胞系统. 2021;12（6）：654–669.e3.密码：34139171

查看文章PubMed/NCBI谷歌学术搜索

4.弗雷施林CR，法尔伯格SA，罗梅罗PA。机器学习导航蛋白质工程的健身环境。Curr Opin Biotechnol.2022;75(102713):102713.密码：35413604

查看文章PubMed/NCBI谷歌学术搜索

5.维特曼BJ，约翰斯顿KE，阿尔姆谢尔PJ，阿诺德FH。EvSeq：对蛋白质文库中的每个变体进行经济高效的扩增子测序。合成生物学学报， 2022;11（3）：1313–1324.pmid：35172576

查看文章PubMed/NCBI谷歌学术搜索

6.Fowler DM，Fields S.深度突变扫描：蛋白质科学的新风格。纳特方法。2014;11(8):801–807.密码：25075907

查看文章PubMed/NCBI谷歌学术搜索

7.罗梅罗PA，阿诺德FH。通过定向进化探索蛋白质健身景观。分子细胞生物学杂志， 2009;10（12）：866–876.密码：19935669

查看文章PubMed/NCBI谷歌学术搜索

8.Biswas S，Khimulya G，Alley EC，Esvelt KM，Church GM。 具有数据高效深度学习的低N蛋白质工程。纳特方法。2021;18(4):389–396.密码：33828272

查看文章PubMed/NCBI谷歌学术搜索

9.Hsu C，Nisonoff H，Fannjiang C，Listgarten J.从进化和测定标记的数据中学习蛋白质适应度模型。纳特生物技术公司。2022;40(7):1114–1122.密码：35039677

查看文章PubMed/NCBI谷歌学术搜索

10.Gelman S，Fahlberg SA，Heinzelman P，Romero PA，Gitter A.从深度突变扫描数据中学习蛋白质序列 - 功能关系的神经网络。美国国家科学院院刊.2021;118（48）：e2104878118.密码：34815338

查看文章PubMed/NCBI谷歌学术搜索

11.里维斯， 迈尔 J， 塞尔库 T， 戈亚尔 S， 林 Z， 刘 J， 等.生物学结构和功能是通过将无监督学习扩展到 250.2021 亿个蛋白质序列而出现的。美国国家科学院院刊118;15（2016239118）：e33876751.密码：<>

查看文章PubMed/NCBI谷歌学术搜索

12.Trinquier J，Uguzzoni G，Pagnani A，Zamponi F，Weigt M.使用简单的自回归模型对蛋白质序列进行有效的生成建模。纳特公社。2021;12(1):5800.pmid：34608136

查看文章PubMed/NCBI谷歌学术搜索

13.雷佩卡 D， 黄尼斯基斯 V， 卡普斯 L， 伦贝扎 E， 罗凯蒂斯 I， 兹里梅克 J， 等.使用生成对抗网络扩展功能蛋白质序列空间。自然机器智能。2021;3(4):324–333.

查看文章谷歌学术搜索

14.Weigt M，White RA，Szurmant H，Hoch JA，Hwa T.通过消息传递鉴定蛋白质 - 蛋白质相互作用中的直接残基接触。美国国家科学院院刊， 2009;106（1）：67–72.pmid：19116270

查看文章PubMed/NCBI谷歌学术搜索

15.Morcos F， Pagnani A， Lunt B， Bertolino A， Marks DS， Sander C， et al.残基协同进化的直接偶联分析可捕获许多蛋白质家族的天然接触。美国国家科学院院刊.2011;108（49）：E1293–E1301.密码：22106262

查看文章PubMed/NCBI谷歌学术搜索

16.Marks DS， Colwell LJ， Sheridan R， Hopf TA， Pagnani A， Zecchina R， et al.根据进化序列变异计算的蛋白质 3D 结构。普洛斯一。2011;6（12）：e28766.pmid：22163331

查看文章PubMed/NCBI谷歌学术搜索

17.Ekeberg M，L?vkvist C，Lan Y，Weigt M，Aurell E.改进蛋白质中的接触预测：使用伪似然推断Potts模型。物理评论 E. 2013;87（1）：012707.密码：23410359

查看文章PubMed/NCBI谷歌学术搜索

18.Kamisetty H，Ovchinnikov S，Baker D.评估基于协同进化的残基 - 残基接触预测在序列和结构丰富的时代的效用。美国国家科学院院刊.2013;110(39):15674–15679.

查看文章谷歌学术搜索

19.霍普夫 TA， 英格拉姆 JB， 波尔韦克 FJ， 舍尔夫 CP， 施普林格 M， 桑德 C， 等.从序列协变预测的突变效应。自然生物技术。2017;35(2):128–135.密码：28092658

查看文章PubMed/NCBI谷歌学术搜索

20.Russ WP， Figliuzzi M， Stocker C， Barrat-Charlaix P， Socolich M， Kast P， et al.用于设计绒毛膜变位酶的基于进化的模型。科学。2020;369(6502):440–445.密码：32703877

查看文章PubMed/NCBI谷歌学术搜索

21.Fantini M，Lisi S，De Los Rios P，Cattaneo A，Pastore A.体外进化的蛋白质结构信息和进化景观。分子生物学和进化。2020;37(4):1179–1192.密码：31670785

查看文章PubMed/NCBI谷歌学术搜索

22.Stiffler MA， Poelwijk FJ， Brock KP， Stein RR， Riesselman A， Teyra J， et al.实验进化的蛋白质结构。细胞系统。2020;10(1):15–24.密码：31838147

查看文章PubMed/NCBI谷歌学术搜索

23.Fernandez-de Cossio-Diaz J，Uguzzoni G，Pagnani A.从深度突变扫描中无监督推断蛋白质适应性景观。分子生物学和进化。2021;38(1):318–328.密码：32770229

查看文章PubMed/NCBI谷歌学术搜索

24.Sohail MS， Louie RH， Hong Z， Barton JP， McKay MR. 从遗传时间序列数据推断上位。分子生物学和进化。2022;39（10）：MSAC199.密码：36130322

查看文章PubMed/NCBI谷歌学术搜索

25.Salverda ML， Dellus E， Gorter FA， Debets AJ， Van Der Oost J， Hoekstra RF， et al.初始突变指导蛋白质进化的替代途径。公共科学图书馆遗传学。2011;7（3）：e1001321.密码：21408208

查看文章PubMed/NCBI谷歌学术搜索

26.弗雷泽 J， 诺廷 P， 迪亚斯 M， 戈麦斯 A， 敏 JK， 布洛克 K， 等.使用进化数据的深度生成模型进行疾病变异预测。自然界。2021;599(7883):91–95.密码：34707284

查看文章PubMed/NCBI谷歌学术搜索

27.Miton CM，Buda K，Tokuriki N.蛋白质进化中的上位和分子内网络。结构生物学的当前观点。2021;69:160–168.密码：34077895

查看文章PubMed/NCBI谷歌学术搜索

28.Ohmae E，Sasaki Y，Gekko K.五色氨酸突变对大肠杆菌二氢叶酸还原酶结构，稳定性和功能的影响。生物化学杂志。2001;130(3):439–447.密码：11530021

查看文章PubMed/NCBI谷歌学术搜索

29.比尔德·瓦，阿普尔曼，黄S，德尔坎普TJ，弗赖斯海姆JH，布莱克利RL。诱变揭示的人二氢叶酸还原酶的保守活性位点残基色氨酸-24的作用。生物化学。1991;30(5):1432–1440.密码：1991124

查看文章PubMed/NCBI谷歌学术搜索

30.Bershtein S， Goldin K， Tawfik DS.强烈的中性漂移产生稳健且可进化的共识蛋白。分子生物学杂志.2008;379(5):1029–1044.密码：18495157

查看文章PubMed/NCBI谷歌学术搜索

31.Van Nimwegen E，Crutchfield JP，Huynen M.突变稳健性的中性进化。美国国家科学院院刊.1999;96(17):9716–9720.密码：10449760

查看文章PubMed/NCBI谷歌学术搜索

32.贝尔施泰因 S， 西格尔 M， 贝克曼 R， 托库里基 N， 陶菲克 DS.稳健性-上位链接塑造了随机漂移蛋白质的适应性景观。自然界。2006;444(7121):929–932.密码：17122770

查看文章PubMed/NCBI谷歌学术搜索

33.Bisardi M，Rodriguez-Rivas J，Zamponi F，Weigt M.蛋白质进化中序列空间探索和上位信号的出现。分子生物学和进化。2022;39（1）：MSAB321.密码：34751386

查看文章PubMed/NCBI谷歌学术搜索

34.奥尔森， 吴新， 孙炳.对整个蛋白质结构域的成对上位的全面生物物理描述。当前生物学。2014;24(22):2643–2651.密码：25455030

查看文章PubMed/NCBI谷歌学术搜索

35.Wu NC， Dai L， Olson CA， Lloyd-Smith JO， Sun R. 间接路径促进了蛋白质健身景观的适应。生活。2016;5：e16965.密码：27391790

查看文章PubMed/NCBI谷歌学术搜索

36.Weinreich DM， Delaney NF， DePristo MA， Hartl DL.达尔文的进化只能遵循很少的突变路径来获得更合适的蛋白质。科学。2006;312(5770):111–114.pmid：16601193

查看文章PubMed/NCBI谷歌学术搜索

37.布莱恩特DH，巴希尔A，西奈S，Jain NK，Ogden PJ，Riley PF等。通过机器学习实现AAV衣壳蛋白的深度多样化。自然生物技术。2021;39(6):691–696.密码：33574611

查看文章PubMed/NCBI谷歌学术搜索

38.Eigen M.物质的自组织和生物大分子的进化。Naturwissenschaften.1971;58(10):465–523.pmid：4942363

查看文章PubMed/NCBI谷歌学术搜索

39.埃斯维尔特， 卡尔森， 刘博士.用于生物分子连续定向进化的系统。自然界。2011;472(7344):499–503.pmid：21478873

查看文章PubMed/NCBI谷歌学术搜索

40.Ravikumar A，Arzumanyan GA，Obadi MK，Javanpour AA，Liu CC.基因在高于基因组错误阈值的突变率下进行可扩展的连续进化。细胞。2018;175(7):1946–1957.密码：30415839

查看文章PubMed/NCBI谷歌学术搜索

41.哈尔佩林， 图CJ， 黄EB， 莫达维C， 谢弗DV， 杜伯JE.CRISPR引导的DNA聚合酶能够在可调窗口中实现所有核苷酸的多样化。自然界。2018;560(7717):248–252.密码：30069054

查看文章PubMed/NCBI谷歌学术搜索

42.温赖特 MJ，乔丹 MI. 图形模型、指数族和变分推理。现在出版公司;2008.

43.Paszke A， Gross S， Massa F， Lerer A， Bradbury J， Chanan G， et al. PyTorch： An Imperative Style， High-Performance Deep Learning Library.在：瓦拉赫H，拉罗谢尔H，贝格尔齐默A，达阿尔奇-布克F，福克斯E，加内特R，编辑。神经信息处理系统进展 32.柯伦联合公司;2019.第8024–8035页。可用： http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.

44.Stratagene，安捷伦科技。基因形态II随机诱变试剂盒;2009. 可从： https://www.chem-agilent.com/pdf/strata/200552.pdf.

45.Shafikhani S，Siegel RA，Ferrari E，Schellenberger V.通过基于PCR的质粒多聚化在枯草芽孢杆菌中生成大型随机突变体文库。生物技术。1997;23(2):304–310.密码：9266088

查看文章PubMed/NCBI谷歌学术搜索

46.Romero PA，Tran TM，Abate AR.通过基于微流体的深度突变扫描解剖酶功能。美国国家科学院院刊.2015;112(23):7159–7164.密码：26040002

查看文章PubMed/NCBI谷歌学术搜索

47.Quan J， Tian J. 复杂基因库和途径的环状聚合酶延伸克隆.普洛斯一。2009;4（7）：e6441.pmid：19649325

查看文章PubMed/NCBI谷歌学术搜索

48.Mago? T，Salzberg SL. FLASH：快速调整短读段的长度以改善基因组组装。生物信息学。2011;27(21):2957–2963.密码：21903629

查看文章PubMed/NCBI谷歌学术搜索

49.Langmead B， Salzberg SL. 与 Bowtie 2 的快速间隙读取对齐。自然方法。2012;9(4):357.密码：22388286

查看文章PubMed/NCBI谷歌学术搜索