DRUGAI
今天为大家介绍的是来自Ben Lehner团队的一篇论文。序列长度为100的蛋白质可能的氨基酸序列组合超过个，而已知宇宙中的原子数远少于这个数量。在如此庞大的序列空间中，无论是实验还是计算，都只能探索其中极小的一部分。近年来，深度神经网络越来越多地用于探索高维序列空间。然而，这些模型极其复杂。作者通过从大于的序列空间中进行实验采样，证明至少某些蛋白质的遗传结构非常简单，可以通过完全可解释的能量模型在高维序列空间中准确预测遗传特性。这些模型能够捕捉自由能与表型之间的非线性关系，但总体上由自由能的加性（additive）变化组成，只有少量的两两能量耦合作用。这些能量耦合作用稀疏，且与结构接触、主链的接近性相关。结果表明，蛋白质的遗传实际上既简单又容易理解。
大规模并行实验可以全面量化蛋白质中单个氨基酸（amino acid，aa）变化的影响。同样，双突变体的实验分析也是可行的，至少对于小型蛋白质来说是如此。然而，高阶突变体的分析由于基因型可能组合的指数爆炸性增长，变得不可行。例如，在蛋白质的34个不同位点上进行单个突变的组合方式为，大约等于。目前的技术下，实验探索如此大量的基因型极具挑战性，迄今为止，实验分析的序列空间约为。
图 1
此外，即使是适量的随机突变组合，总会产生非功能性蛋白质。例如，在一个小型蛋白质结构域中，假设能量组合是加性的，只有2%到8%的5个氨基酸变体和不到0.2%的10个氨基酸变体预计会正确折叠（n = 2个结构域；图1a）。因此，即使在大多数蛋白质中采样数千万的随机组合基因型，也几乎无法提供有关遗传结构的信息——即突变如何组合以决定表型的规则。这对训练和评估预测模型并没有太大帮助，除了验证大多数基因型是未折叠的这一简单预测外。
探索高维序列空间的一种策略是使用深度学习。具有数百万个拟合参数的深度神经网络在各种预测和蛋白质设计任务中取得了成功，包括预测组合突变体的效果。然而，这些模型的结构极其复杂，难以解释。可能是蛋白质基因型-表型的关系图非常复杂，突变之间需要有大量的相互作用才能进行准确预测。或者，正如能量测量和推断所表明的，这些关系图可能要简单得多，并且可以使用统计模型进行解释。对作者来说，简单模型是指具有少量参数的模型（因此能够实现大规模的数据压缩）以及具有可解释参数的模型（因此能够提供理解）。
作者使用了一种实验设计，通过富集功能性蛋白质序列，来探索维度超过30、基因型数量超过的高维蛋白质序列空间的遗传结构。研究发现，蛋白质的结构非常简单，基于加性能量模型能够提供非常好的预测性能。量化突变之间的两两能量耦合（energetci coupling）进一步提高了预测能力，使其在高维基因型空间中表现出色。这些耦合作用较为稀疏，并且与蛋白质的三维结构相关。因此，至少某些蛋白质的遗传结构非常简单，表现为加性能量作用以及来自稀疏两两结构耦合的小部分贡献。
在序列空间中采样
作者之前展示了，通过合成、筛选和测序实验的组合，可以大规模测量数千个单个突变对蛋白质稳定性的能量效应。在这些实验中，测量了每个突变对野生型蛋白质以及少量具有不同折叠稳定性变体的细胞丰度的影响。例如，利用一个浅层双突变体库，作者能够推断适配蛋白GRB2的C末端SH3结构域中几乎所有突变（1,056/1,064 = 99%）的折叠吉布斯自由能变化（ΔΔGf）。类似的大规模并行测量现已应用于其他信号结构域的单突变折叠稳定性测定，包括癌蛋白KRAS，以及在体外测定的许多主要为原核生物的小型结构域。
在GRB2-SH3结构域中随机组合突变会很快导致蛋白质不能折叠，预计当有五个和十个突变时，分别约有98%和超过99.9%的基因型无法折叠（基于加性能量，见图1a）。这种随着突变组合导致的稳定性快速下降，跟在其他蛋白质上的实验测量结果一致。为了在高维序列空间中探索能折叠的基因型，作者使用了一种启发式技术，来富集组合变体中保守的折叠和功能。对于每一个可能的初始单一氨基酸替换，作者逐步选择进一步的替换——每个残基位置一个——以同时最大化1、组合突变体的预测丰度；2、相互作用伙伴的结合。对于GRB2-SH3，预测能够保留两种分子表型的最大突变集包含34个单一氨基酸替换：25个位于表面残基（相对溶剂可及表面积（RSASA）≥ 0.25），3个位于蛋白质核心（RSASA < 0.25），6个位于GAB2配体结合界面（配体距离<5Å；见图1b，右）。
作者合成了一个包含这34个突变组合的库（“库1”），并使用高度验证的组合选择和丰度蛋白质片段互补测定（AbundancePCA）对约种基因型中的部分样本进行了细胞丰度的定量测量。总计获得了129,320个变体的三重复丰度测量数据，这占整个序列空间的0.0007%。这些测量具有高度的可重复性（Pearson相关系数r > 0.91；见图1d）。
基因型频率分布的对称豆荚状形态，在实验采样的文库中得到了再现，基因型数量在与野生型（零阶）和第34阶突变等距的中间汉明距离17处达到峰值（图1e）。随着氨基酸替换数量的增加，丰度的中位数测量值有所下降，但仍有数千个带有大量突变的基因型，其丰度得分与野生型蛋白无显著差异（n = 2,706，带有超过20个突变，双侧z检验，名义P值 > 0.05；见图1f）。
利用能量模型做基因预测
图 2
量化大量多重突变体的效应使作者能够在遗传景观的更广泛区域中测试基因型-表型模型的预测性能，而不仅限于训练时使用的局部邻域。为了构建和评估模型，作者将所有分析限制在所有三个生物重复中都有定量测量的变体上（n = 71,233）。值得注意的是，作者最初的能量模型（图2a）是基于丰度和配体结合选择（doubledeepPCA，ddPCA）训练的，只量化了单一和双氨基酸突变的效应，这一模型解释了多重组合突变体适应性方差的约一半（R² = 0.5；见图2b右下角），其中大多数（94%）包括了野生型序列中的至少13个氨基酸替换。
在这个简单模型中，唯一训练的参数是野生型的吉布斯自由能（ΔGf）和单一氨基酸替换的自由能变化（ΔΔGf），以及一个将折叠分子比例与AbundancePCA得分（适应性）联系起来的两参数（仿射）变换（见图2a）。如此大比例的表型方差由基于加性能量模型解释（无特定的表观效应或遗传相互作用），该模型仅在包含一到两个基因改变的基因型上训练，表明突变对蛋白质的能量效应在很大程度上与背景无关。
另一方面，一个线性模型——即隐含假设在多突变体中突变效应在表型水平是加性组合的——基于相同的ddPCA数据训练，其表现明显较差（R² = 0.32）。该线性模型系统性地低估了突变组合的实际表型效应（图2b左下角），这是由于没有考虑蛋白质热力学引起的突变效应缩放（即全局表观效应）的结果。例如，在已经未折叠的蛋白质中引入一个去稳定化突变，对折叠分子的比例没有影响（见图2c中模型的下平台），而这一点无法通过线性模型捕捉。这些结果表明了拟合能量模型的一个关键优势：考虑全局表观效应可以提升预测在训练数据局部邻域之外的泛化能力。
将线性模型和能量模型拟合到组合数据中，分别提高了30%和13%的解释方差（图2b，上方图），这可能是由于训练数据量的增加以及更相关的遗传背景，在这些背景中量化了每个单一氨基酸的效应：库中大约50%的变体（n ≈ 30,000）包含并报告了34个单一氨基酸替换中的任意一个效应，这意味着每个单突变的测量数据比相对较浅的ddPCA文库多了近三个数量级。虽然一阶线性模型和能量模型解释的方差比例相当（R² = 0.62 和 0.63；见图2b，上方图），但线性模型的回归残差偏差表明该模型不太合适。能量模型很好地拟合了数据，准确捕捉了观察到的AbundancePCA适应性与推断的折叠自由能变化（ΔGf）之间的全局非线性关系（全局表观效应）（图2c）。此外，从组合数据集和ddPCA数据集中推断的自由能变化（模型参数）之间的相关性也非常高（Pearson相关系数r = 0.87），但前者的估计值往往更极端，再次证明了在更多遗传背景中测量突变效应的价值，从而使其能量效应的估计更加准确（图2d）。
耦合作用提升了遗传预测
接下来，作者测试了量化突变之间的非加性能量耦合是否能提升预测性能。在组合数据集中，每对突变平均出现在17,923种基因型中，从而能够稳健地测量二阶遗传相互作用项（能量耦合，ΔΔΔGf）。纳入所有二阶能量耦合作用后，模型性能提升了额外的9%（R² = 0.72），这与预期一致，即两两效应是蛋白质中特定表观效应的重要来源（见图2e）。尽管一阶项效应强度更大，且偏向于去稳定化效应，二阶能量耦合的效应通常较温和，且以零为中心（见图2f）。
物理接触和主链邻近性
图 3
以表型水平测量，蛋白质中的遗传相互作用此前已被证明至少部分反映了蛋白质结构。结合组合深度突变扫描和热力学建模，作者推断出总共561对能量耦合，为探究其机制起源及其与蛋白质结构的关系提供了机会。将耦合能量大小（绝对折叠ΔΔΔGf）与折叠结构中突变对的三维距离（最小侧链重原子距离）进行比较，结果显示出一个L形分布，最强的能量耦合发生在结构上邻近的残基之间（图3a）。能量耦合最强的前五对残基都位于5.5Å以内，而前20对中有15对（75%）的残基距离小于8Å。尽管接触距离与耦合能量强度之间存在弱负相关关系（Spearman相关系数ρ = -0.12；见图3a），但这种趋势在主链序列上不相邻的残基对中不明显（Spearman相关系数ρ = -0.02，主链距离>5个残基）。
另一方面，将耦合强度与主链序列中残基之间的分离距离（沿肽链主链）进行比较，显示出明显的负相关关系，并且这种关系在较大距离上依然存在（Spearman相关系数ρ = -0.28），即使排除残基间的直接物理接触（<5Å，Spearman相关系数ρ = -0.27；见图3b）。图3c中的相互作用矩阵总结了这些观察结果：最强的能量耦合与直接物理接触一致（黑色圆圈；见图3d），并且能量耦合强度沿蛋白质主链逐渐衰减（图3c，靠近对角线与远离对角线的单元格对比）。矩阵还突出了二级结构元件之间的物理相互作用，作为强能量耦合的热点。
为了理清这些不同潜在结构因素对能量耦合强度的相对重要性，作者收集了一组定量特征，描述了残基对原子之间存在的化学键或相互作用的数量和类型，以及它们在折叠结构中的相对位置（图3e）。基于这12个结构特征的线性回归模型可以预测耦合强度（图3f）。值得注意的是，当该模型在独立的组合突变实验中推断的能量耦合（“库3”，）上进行评估时，表现也同样良好（Pearson相关系数r = 0.46，R² = 0.21）。这表明尽管模型较为简单，但它能够捕捉到能量耦合强度的结构决定因素，并且能量耦合作用是由结构相互作用引起的。
高阶突变体的折叠和功能
图 4
作者的实验鉴定了大量GRB2-SH3基因型，尽管带有多次突变，但仍然具有较高的细胞丰度（例如，包含超过五个突变的25,564种基因型；见图1f）。为了进一步确认这些高丰度的多突变体是正确折叠并具有功能的，作者进行了第三次组合突变实验，并使用蛋白质-蛋白质相互作用测定（BindingPCA）测试了GRB2-SH3变体与肽配体结合的能力（见图4a）。只有当蛋白质采用其天然构象时，才能识别肽配体（见图1b、c）。作者设计了一个包含15个单一氨基酸替换的文库（库3），这些替换发生在一个22个氨基酸的窗口内，避免了原始文库中与配体结合界面相关的突变（侧链重原子距离配体< 5Å）。该文库包含32,768个变体（），并且与原始的个库共享6个单一氨基酸替换。总共获得了25,967个变体的结合测量和31,936个变体的丰度测量，分别占整个序列空间的79%和97%。这些测量具有高度的可重复性（结合和丰度的Pearson相关系数分别为r > 0.85和r > 0.94；见图4c）。
将结合能力的变化与第三阶、第六阶和第九阶变体的丰度变化绘制在一起，可以看出，大多数影响结合的突变也影响了孤立结构域的浓度，这与之前的研究结果一致，并符合蛋白质稳定性变化是突变对结合影响的主要原因这一预期（见图4d）。然而，值得注意的是，大多数高丰度的高阶突变体也能结合GAB2配体，表明它们是正确折叠的（见图4d）。例如，4%（204/4,805）的包含九个突变的变体的丰度与野生型蛋白无显著差异（名义P > 0.05），其中96%（177/184）也能结合配体（预测结合分子比例>0.5）。因此，大多数高丰度的高阶GRB2-SH3突变体是正确折叠的。
多表型遗传预测
在大量遗传背景中，测量了这两个相关分子表型的单氨基酸和双氨基酸突变效应，这为热力学建模提供了丰富的数据来源。首先，仅考虑丰度表型时，作者观察到一个具有未折叠和折叠能量状态的加性双态热力学模型，在验证集中表现优于线性模型（R² = 0.93 对比 R² = 0.87）。要实现与一阶能量模型相似的预测性能，线性模型需要纳入二阶和三阶遗传相互作用项，这大大增加了模型的复杂性（715个参数对比16个参数，增加了40倍以上）。这种使用大量特定的两两和高阶遗传相互作用项来捕捉数据中的全局非线性（全局表观效应）的模型复杂性，被称为“幻影表观效应”。
接下来，作者在之前的研究基础上，使用了一个神经网络实现的三态平衡模型，其中包含未折叠、折叠和结合的能量状态（图4e），以同时推断折叠和结合的自由能变化（ΔΔGf 和 ΔΔGb），以及折叠和结合的能量耦合项（ΔΔΔGf 和 ΔΔΔGb）（图4f）。该模型与数据非常吻合（图4g），几乎解释了所有的适应性方差（图4h），并且推断出的折叠和结合自由能变化（一阶项）与之前使用独立的浅层双突变文库获得的结果高度相关（Pearson相关系数分别为r = 0.9 和 r = 0.7）（图4i，双突变）。据作者所知，这是首次测量到大量蛋白质的折叠（ΔΔΔGf）和结合（ΔΔΔGb）能量耦合项。
编译 | 黄海涛
审稿 | 王梓旭
参考资料
Faure, A. J., Martí-Aranda, A., Hidalgo-Carcedo, C., Beltran, A., Schmiedel, J. M., & Lehner, B. (2024). The genetic architecture of protein stability. Nature, 1-9.
