DRUGAI
今天为大家介绍的是来自澳大利亚国立大学C.J.Jackson团队的一篇论文。蛋白质语言模型(PLMs)将氨基酸序列转换为训练机器学习模型所需的数值表征。许多PLMs规模庞大(超过6亿个参数)且在广泛的蛋白质序列空间上进行训练。然而，这些模型在预测准确性和计算成本方面存在局限性。在此，作者使用多重祖先序列重建方法来生成小型但功能性聚焦的蛋白质序列数据集用于PLM训练。与大型PLMs相比，这种局部祖先序列嵌入产生了具有更高预测准确性的表征。作者表明，由于祖先序列重建数据的进化特性，局部祖先序列嵌入产生了更平滑的适应度景观，在这种景观中，适应度值相近的蛋白质变体在表征空间中的数值距离也更接近。这项工作有助于在数据稀疏且计算资源有限的现实环境中实施基于机器学习的蛋白质设计。
为了在没有功能标签的情况下构建能够捕获蛋白质生物物理和进化特征的表征，PLMs通常通过无监督掩码语言建模(MLM)进行训练，在这种方法中，模型的任务是预测已从周围序列上下文中被掩码的残基的身份。经过训练的PLM的隐藏状态被提取并合并为表征，隐式地捕获了蛋白质序列的相关物理和生物学特性。MLM任务本身在蛋白质工程和进化推断中也显示出应用价值。
蛋白质表征也可以通过物理氨基酸描述符的正交主成分经验性推导获得：疏水性、溶剂可及表面积或电荷。然而，与深度表征模型不同，这种线性表征无法捕获上下文依赖性。因此，像PLMs这样的深度表征模型具有更强的上下文感知能力，这大大提高了作者捕获复杂生物现象的能力。
近期研究表明，经过显式潜在空间正则化训练的深度神经网络能够找到蛋白质空间的平滑表征，而那些考虑上位效应系数（epistasis coefficients）的网络改善了表征在下游监督任务中的使用。这意味着平滑的嵌入空间(相对于蛋白质适应度而言)可能比非平滑的空间更具信息量。对于没有显式潜在空间正则化的无监督PLMs，表征空间中的任何隐式平滑都是通过从MLM目标学习上下文序列依赖关系而实现的。因此，作者提出，在能最大化上位效应预测理解的序列数据上训练的模型可能会产生高度信息丰富的蛋白质表征。
祖先序列重建(Ancestral sequence reconstruction，ASR)是一种统计方法，用于从系统发育树的内部节点推断已灭绝的分子序列。由ASR产生的序列通常具有功能性，并经常表现出适合蛋白质工程的新表型或特性。因此，ASR被广泛用于生成热稳定蛋白质、探索功能序列空间以及设计用于工程改造的蛋白质支架。事实上，最近的研究表明，在功能性蛋白质生成方面，ASR的表现优于包括大型PLMs在内的最先进的深度神经网络。ASR还为分子进化过程和理解序列-功能关系提供了见解。研究已灭绝祖先序列的生物物理、化学和生物学特性可以揭示蛋白质获得新表型和功能的机制，以及那些混淆或促进这些机制的特征，如上位效应。在本研究中，作者描述了使用ASR生成的序列来训练家族特异性蛋白质表征模型。
多重ASR
训练深度表征模型需要覆盖功能序列空间不同区域的大型序列数据集。对于具有n个末端的完全解析和有根系统发育树，如果仅从每个内部节点的后验概率分布中采样最可能(最大后验概率，MAP)的序列，最多可以重建n-1个祖先序列。这是因为ASR依赖于先验系统发育假设，这种假设本质上将ASR生成的序列多样性限制在潜在的树结构(和序列进化模型)内。
在此，作者开发了一种称为"多重ASR"(mASR)的方法，该方法采样统计等效的拓扑结构作为ASR的先验，以生成用于表征学习的大型且多样化的序列数据集。由于真实的系统发育关系永远无法被确定地获知，系统发育拓扑结构是通过启发式树搜索算法重建的，这些算法使特定树的负对数似然值最小化。通过并行执行多个独立的树搜索，并通过近似无偏检验过滤掉统计上不等效的树，作者生成了一组同样有效但不同的系统发育树，用于重建祖先序列。采样序列空间的规模意味着随机独立树搜索很少返回相同的树拓扑结构，这有效地将ASR生成的序列数量增加了一个因子，该因子等于被近似无偏检验接受的树搜索重复次数。
传统ASR方法在生成功能序列时的另一个限制是插入和缺失(indels)事件被最大似然序列进化模型错误处理或忽略。由于indel事件经常驱动功能多样化，作者还开发了一个自动化流程，用于对mASR生成的大型祖先序列数据库中的indels进行最大似然重建。对于位点n，所有末端都被分配一个二进制标签（0表征无indel，1表征有indel），这取决于是否在特定位点观察到缺口字符。作者假设序列中的插入率大约等于缺失率，并使用等速率模型重建所有祖先节点位点n处存在缺口的概率。作者对在现存序列中≥1%观察到缺口的所有位点重复这个过程（假设位点独立），并移除所有重建标签<0.5的祖先节点中的残基。作者的最大似然indel处理与mASR相结合，使作者能够从进化信息生成大型的真实序列数据集。
在此，作者研究了上述indel处理流程的mASR在三个不同蛋白质系统上的应用：细菌磷酸三酯酶(PTE)、细菌聚对苯二甲酸乙二醇酯(PET)水解酶和真核咪唑甘油磷酸脱氢酶。细菌PTE家族，是一组水解磷酸三酯类农药的同源蛋白，已被广泛用作功能适应和分子进化的模型。PET水解酶(PETase)是一组最近发现的细菌和环境宏基因组酶，它们可以水解PET塑料，并已经过广泛的工程改造，包括最近描述使用mASR探索PETase序列空间新区域的工作。最后，咪唑甘油磷酸脱氢酶(酿酒酵母中的His3p；在此称为His3p)是真核酶，对内源性组氨酸生物合成至关重要。
从293个非冗余现存PTE同源蛋白的初始比对开始，尽管大多数重复返回高度相似的拓扑结构，作者通过100次树搜索重复生成了10,148个独特的PTE类序列数据集（图1a）。作者将祖先和现存PTE序列嵌入到PLM ESM-1b中，该模型在UniRef50-S数据库中约2.5亿个非冗余蛋白质序列上通过MLM预训练，以研究mASR生成的序列多样性。然后作者使用t-分布随机邻域嵌入(t-SNE)将现存和祖先PTEs的ESM-1b表征投影到二维空间进行可视化（图1b）。通过mASR生成的大多数祖先PTE序列属于未被用于重建它们的现存PTE同源物采样的序列空间区域。在His3p上的mASR结果类似，通过100次树搜索重复从331个现存序列生成了7,864个不同的序列。PETase上mASR的使用最近已在其他地方描述，简而言之，与His3p和PTE的结果一致。
图 1
因此，考虑到重建的祖先蛋白质通常表现出与其模板序列相当的催化活性和增加的热稳定性，mASR似乎是序列新颖性的重要来源。由于祖先序列在大型预训练PLMs中没有表征，因此专门设计用于蛋白质工程中新序列生成的生成式PLMs不太可能产生类祖先序列。这些结果表明，mASR（包括最大似然indel重建）是一种独特且强大的统计方法，用于生成大型祖先序列数据集。
LASE模型
作者使用mASR生成的数据集来训练用于表征学习的小型transformer编码器（图1c）。Transformer由一个位置嵌入层、六个编码器块（每个块包含一个四头多头注意力层和一个前馈全连接层）以及一个时间分布式全连接输出层组成。嵌入和前馈维度分别为128和512，所有模型都使用MLM目标进行训练。Transformer的编码器块从MLM目标的训练数据中学习长程序列依赖关系；这些层的隐藏状态捕获了被嵌入的蛋白质序列的物理、化学和生物学特性，并被合并成密集的固定长度蛋白质表征，这与之前基于transformer的PLMs相似。这些表征模型的小规模（约230万参数）使其训练在计算上成本较低，在常规硬件（单个Nvidia A10）上每个epoch需要不到1小时。这与大型PLM微调形成鲜明对比，后者需要在家族或任务特定数据上重新训练所有参数，或采用尖峰数据调度，其中微调数据集大小被不相关的蛋白质序列稀释100倍以避免遗忘通用蛋白质特征。作者将这种方法称为"局部祖先序列嵌入"(LASE)，因为它使用合成的祖先序列数据为通过mASR采样的局部序列空间构建信息丰富的表征空间。
LASE模型预测性能
图 2
接着,作者测试了各个系统特异性LASE模型的预测性能。PTE变体对2-萘基己酸酯(2NH)的芳基酯酶活性此前已在四个独立实验中得到描述。其中三项研究是通过定向进化获得的:两条"前向"进化轨迹,即将PTE进化为增强的芳基酯酶活性(R轨迹和S轨迹),以及一条"反向"轨迹,即将R轨迹末端获得的最进化的芳基酯酶"反向"进化以重获PTE活性。第四个数据集是通过对PTE活性位点附近的六个氨基酸残基进行组合突变来研究上位效应而生成的。类似的方法最近也被开发并用于研究二氢叶酸还原酶的适应度景观。为了生成标准训练数据集,通过对所有实验中共有的数据点进行活性测量值的归一化,将每个实验的数据进行了整合。为了减少由细胞生长差异引起的蛋白质表达噪声,所有澄清裂解液的活性数据都通过细胞密度(在600 nm处测量的光密度或OD600)进行了归一化。PTE表达不太可能显著影响相对活性测量值(超出催化效率的真实变化),因为尽管积累了多个突变,PTE变体仍保持着相当的(大于两倍的)溶解度水平,而相对活性则相差几个数量级。最终PTE数据集的结构被描绘为汉明图(图2b),其中可以看到两种不同的数据结构:轨迹数据,即突变沿着线性定向进化轨迹积累,以及组合数据,即探索六个活性位点替换的所有组合。PETase和His3p数据集各自来自单个实验,不需要跨实验归一化。PETase数据集包含了在表皮蛋白酶序列空间中稀疏采样的196个祖先PETase变体(图2c),而His3p数据集包含了2,028个野生型S. cerevisiae His3p的组合突变体(图2d)。因此,作者的分析对生物学上不同背景(原核生物异源物降解酶和必需的真核生物代谢酶)采样了实验上多样化的场景(定向进化、随机序列采样和完整的组合突变)。
与其他模型进行性能比较
LASE与以下几种方法进行了基准比较:使用随机权重初始化的LASE架构(LASE_random)、在等量现存序列上训练的LASE架构(LASE_extant)、独热编码、在大型数据库上训练的蛋白质语言模型所学习的表征(UniRep、ProtTrans、ESM-1b和ESM-2),以及从氨基酸理化性质衍生的经验表征(Georgiev、Z-scale、ProtFP)。每个表征模型的性能都通过皮尔逊R2相关系数来衡量,该系数用于计算每个酶系统测试集中观察值与预测适应度值之间的相关性。作者使用皮尔逊R2(这是蛋白质工程中的标准度量)来评分预测值和观察到的适应度值之间的线性关系,因为蛋白质变体之间的等级和绝对差异都很重要。对于PETase和His3p系统,模型性能被报告为在五重交叉验证训练-测试数据分割上确定的平均R2值,每次分割都在80%的监督数据上训练模型,并在剩余20%未参与训练的数据上测试模型。为了确保这种方法的预测性能能够反映现实世界中训练和测试数据来自多个实验的蛋白质工程案例,作者合成并测定了26个在定向进化轨迹中未被采样的PTE突变体(图2a),并比较了预测和观察到的芳基酯酶催化效率。在这个PTE测试集中的每个序列都携带着定向进化世代n和n+1之间的中间单一突变,这适用于所有相隔超过一个突变的世代。作者对每个蛋白质系统的监督回归任务都尝试了一系列不同的模型架构(表1)。
表 1
当使用0.5%的遮蔽目标进行训练时,LASE在预测PTE对芳基酯的催化效率方面表现出最具信息量的表征,在测试数据上获得了0.76的皮尔逊R2值。目前最先进的ESM表征模型ESM-1b和ESM-2在催化效率预测方面的表现明显不如LASE,分别获得了0.69和0.62的皮尔逊R2值。作者在PETase的LASE和ESM表征的预测性能之间观察到了类似的趋势。事实上,尽管ESM-1b比LASE表征有多300倍以上的参数(ESM-1b有6.7亿参数并在2.5亿个序列上训练,而LASE有230万参数并在1万个序列上训练;图3c),在相同的训练-测试数据分割上,使用30%遮蔽训练的LASE模型在祖先表皮酯酶序列上获得了平均0.50的皮尔逊R2值,而ESM-1b表征仅为0.09(Tukey honestly significant difference（HSD）P值=0.09)。
图 3
相比之下,ESM模型和LASE在His3p上的性能差距很小。在这里,ESM-1b略微优于His3p的LASE表征(在五个测试-训练分割上的平均皮尔逊R2分别为0.86和0.89(Tukey HSD P=1))。值得注意的是,对于His3p,所有嵌入方案都表现相对较好,包括预期表现较差的对照组LASErandom。这表明His3p数据集作为一个大型组合数据集,对嵌入选择并不敏感。相反,回归模型能够利用详尽数据集中存在的强序列到适应度信号。
基于此,作者提出LASE在PETase和PTE系统上相对于His3p的性能差异可能源于用于训练ESM-1b的序列语料库中相关序列信息的表征,以及监督数据和蛋白质工程任务的结构。确实,PTE和PETase都是罕见的异物降解酶进化案例,可能是在蛋白质宇宙中的一次性场景中进化而来。因此,这些蛋白质类别在ESM-1b/2训练数据中相对代表性不足(因此在UniRef数据库上训练的所有其他大型PLM中也是如此)。这强调了使用mASR生成合成训练数据的价值,这些数据在查询序列(或查询序列组)周围的序列空间中局部密集,而这些序列在自然蛋白质进化过程中可能被稀疏采样。此外,His3p的监督任务在野生型序列周围的少量突变上是密集且组合完整的。这与cutinase和PTE监督数据集形成对比,后者对序列空间的更大区域进行采样。LASE捕获的进化信息很可能在与祖先序列MLM训练数据相当的尺度上转化为更具信息量的表征,而不是密集的单一突变尺度。这使得LASE特别适用于那些实验采样序列稀疏且多样,或者属于在其他PLM训练数据或现存序列数据库中代表性不足的家族的蛋白质工程任务。
LASE表征模型的性能和MLM训练期间使用的遮蔽百分比存在显著差异。例如,使用15%遮蔽目标训练的PTE LASE模型表现明显较差(皮尔逊R2=0.54),不如使用0.5%遮蔽目标训练的等效表征模型(皮尔逊R2=0.76),这表明遮蔽百分比是应该在PLM训练过程中优化的超参数。作者发现,0.5%的遮蔽目标在PTE和His3p数据集上优于使用5%、15%和25%遮蔽目标训练的架构相同的模型。遮蔽0.5%的训练序列(大约两个位点)在直观上可能提供能够捕获比传统15%遮蔽(每个训练序列~50个位点)目标更丰富的测试数据集规模信息的表征。事实上,当LASE模型使用30%遮蔽目标训练时,PETase表征最具信息量,这与其测试数据集比PTE和His3p更加多样化相吻合(图2)。此外,作者观察到使用较低遮蔽百分比训练的模型性能具有更大的敏感性，这表明虽然较低的遮蔽百分比可能会产生更具信息量的表征,但要达到能提供这种表征的参数空间可能是一个更具挑战性的任务。
平滑的表征空间促进预测性能
在确定LASE的表征在PTE、PETase和His3p的监督任务中超越或可与基于大量现存序列训练的模型的性能相媲美后，作者进一步探讨了各模型表征空间的结构与预测性能之间的关系。作者提出，模型嵌入空间中观察到的适应性越平滑，该模型表征在下游监督任务（通过回归顶层模型）的预测性能可能越好。该假设基于以下前提：能够将序列嵌入为平滑适应性函数的表征，比起那些生成复杂表征空间（表征向量与适应性之间存在复杂关系）的表征，需要更简单的变换即可映射到适应性域。实际上，已有研究表明，明确地对transformer的表征空间进行适应性平滑处理，可以在下游回归任务中生成更有信息量的表征。
作者采用了光谱图方法，评估PTE数据集中每个模型表征嵌入下观察到的适应性是崎岖还是平滑的。具体而言，使用归一化的Dirichlet能量作为嵌入PTE数据集的k近邻（KNN）图上平滑度的度量。KNN图中的节点表征嵌入的序列坐标，边连接每个嵌入序列之间的k个最近邻。通常，Dirichlet能量是衡量函数崎岖（非平滑）程度的有效指标。在此，局部Dirichlet能量表征KNN图中相邻节点之间观察到的适应性差值的平方；如果在KNN图中，邻近节点映射到相似的数值适应性，则Dirichlet能量较低，可以认为嵌入的适应性景观是平滑的。
在PTE的所有测试表征中，作者发现表征嵌入的KNN图的归一化Dirichlet能量与随机森林回归模型预测性能之间存在负相关关系（皮尔逊R2相关系数；R2 = 0.67，P = 0.0016）（图4a,b）。这种关系在PETase系统中也很明显。这表明表征空间的平滑度（相对于适应性）与PTE数据集中更具信息量的表征相关。LASE生成了最平滑的嵌入KNN图，这与PTE测试集中芳基酯酶活性预测的性能一致（表1）。
图 4
当这些表征投影到二维t-SNE基座以进行可视化时，ESM-1b的投影与底层的Hamming图结构相似：具有高芳基酯酶活性的PTE变体从S轨迹中被单独投影，与属于R轨迹的变体分开（图4c）。相比之下，LASE嵌入空间的投影（图4d）根据适应性将这些功能变体聚集为邻近的群体，尽管它们具有不同的进化历史。
编译|黄海涛
审稿|王梓旭
参考资料
Matthews, D. S., Spence, M. A., Mater, A. C., Nichols, J., Pulsford, S. B., Sandhu, M., ... & Jackson, C. J. (2024). Leveraging ancestral sequence reconstruction for protein representation learning. Nature Machine Intelligence, 6(12), 1542-1555.
