作者研究了从蛋白质骨干原子坐标预测蛋白质序列的问题。迄今为止,机器学习解决此问题的方法一直受限于可用的实验测定蛋白质结构的数量。作者使用AlphaFold2为1200万个蛋白质序列预测的结构,从而将训练数据扩充了近三个数量级。相比现有方法,整体提升了近 10 个百分点。该模型还能推广应用于多种更复杂的任务,包括蛋白质复合物设计、部分mask结构预测、结合界面设计以及多状态预测。
来自:Learning inverse folding from millions of predicted structures, ICML2022
目录
- 背景概述
- 从预测结构中学习逆折叠
- 数据
- 模型架构
- 训练
- 结果
- 固定骨架的蛋白质设计
- 零样本预测
背景概述
设计能编码具有所需特性蛋白质的新型氨基酸序列,即从头蛋白质设计,是生物工程中的一项核心挑战。解决这一问题最成熟的方法是使用能量函数,该函数直接对蛋白质折叠状态的物理基础进行建模。
最近,一类基于深度学习的新方法被提出,这些方法利用生成模型来为特定结构预测序列、生成蛋白质骨干结构、联合生成结构和序列,或者直接对序列进行建模。深度学习生成模型具有直接从数据中学习蛋白质设计规则的潜力,这使得它们有望成为当前基于物理原理的能量函数的替代方法。
然而,实验测定的蛋白质结构数量相对较少,这限制了深度学习方法的应用。实验测定的结构覆盖的蛋白质序列空间不到已知序列空间的0.1%。
在此,作者探究预测结构是否可用于克服实验数据的局限性。随着蛋白质结构预测领域的进展,如今已能够考虑大规模地从预测结构中学习。为大型数据库中的序列预测结构,可将蛋白质序列的结构覆盖范围扩大几个数量级。为了训练用于蛋白质设计的逆向模型,作者使用AlphaFold2为UniRef50中的1200万个序列预测结构。
作者专注于从蛋白质骨干结构预测序列,这一问题也被称为逆折叠(inverse folding)或固定骨干设计(fixed backbone design)。作者将逆折叠视为一个序列到序列的问题,采用自回归编码器-解码器架构。在这个架构中,模型的任务是根据蛋白质骨干原子(backbone atoms)的坐标恢复其天然序列。
- 图1:利用预测结构增强逆折叠研究。为了评估使用预测结构训练蛋白质设计模型的潜力,作者使用AlphaFold2为1200万个UniRef50蛋白质序列预测结构。训练一个自回归逆折叠模型来进行固定骨干蛋白质序列设计。训练集和测试集在拓扑结构层面进行划分,以便在结构独立的骨干上对模型进行评估。作者将具有不变几何输入处理层的Transformer模型与先前工作中使用的完全几何模型进行比较。对输入坐标应用span masking(随机选择连续的一段区域进行mask)和噪声处理。
作者利用大量结构未知的序列,将它们作为额外的训练数据。在没有实验测定结构时,以预测结构为条件来训练模型(见图1)。这种方法类似于机器翻译中的反向翻译,在机器翻译中,一个方向的预测译文被用于改进相反方向的模型。研究发现,即使预测的输入(即结构)质量较低,反向翻译也能有效地从额外的目标数据(即序列)中学习。
作者发现,现有方法受到数据的限制。目前最先进的逆折叠模型在使用预测结构扩充训练数据时性能会下降,然而更大规模的模型和不同的模型架构能够有效地从这些额外数据中学习,这使得在结构上独立的天然骨干序列恢复率提高了近10个百分点。依据先前研究中的固定骨干设计基准对模型进行评估,并在一系列任务中评估其泛化能力,这些任务包括蛋白质复合物和结合位点设计、部分mask骨干结构预测以及多构象分析。作者还进一步探讨了这些模型在零样本预测方面的应用,即预测突变对蛋白质功能和稳定性、复合物稳定性以及结合亲和力的影响。
从预测结构中学习逆折叠
逆折叠的目标是设计出能够折叠成所需结构的序列。在这项研究中,我们专注于蛋白质的骨干结构,暂不考虑侧链。虽然20种氨基酸各自都有特定的侧链,但它们都有一组共同的原子构成氨基酸骨干。在这些骨干原子中,选取氮(N)、α-碳原子(Cα)和碳原子(C)的坐标来代表骨干结构。
- 图2:蛋白质设计任务示意图。
利用天然存在的蛋白质结构,我们可以通过监督模型从其骨干原子在三维空间中的坐标预测蛋白质的天然序列。形式上,这个问题表示为学习条件分布 p ( Y ∣ X ) p(Y | X) p(Y∣X) 。对于长度为 n n n的蛋白质,给定结构中每个骨干原子(N、Cα、C)的空间坐标序列 X = ( x 1 , . . . , x i , . . . , x 3 n ) X = (x_{1}, ..., x_{i}, ..., x_{3n}) X=(x1,...,xi,...,x3n),目标是预测氨基酸的天然序列 Y = ( y 1 , . . . , y i , . . . , y n ) Y = (y_{1}, ..., y_{i}, ..., y_{n}) Y=(y1,...,yi,...,yn) 。这种概率密度通过序列到序列的编码器-解码器进行自回归建模: p ( Y ∣ X ) = ∏ i = 1 n p ( y i ∣ y i − 1 , . . . , y 1 ; X ) p(Y | X)=\prod_{i = 1}^{n} p(y_{i} | y_{i - 1}, ..., y_{1} ; X) p(Y∣X)=i=1∏np(yi∣yi−1,...,y1;X)通过最小化负对数似然来训练模型。
数据
预测结构
作者为UniRef50中的序列生成了1200万个结构,以探究预测结构如何改进逆折叠模型。为了选择用于结构预测的序列,首先使用MSA Transformer预测所有UniRef50序列的多序列比对(MSA)距离图。我们根据距离图的LDDT分数对序列进行排序,以此作为预测质量的指标。选取了长度不超过500个氨基酸的前1200万个序列,并使用AlphaFold2模型对它们进行正向折叠,最后采用Amber力场进行优化松弛。最终得到的预测数据集规模约为实验结构训练集的750倍。
训练与评估数据
作者在CATH数据库中一个结构独立的子集上对模型进行评估。在拓扑结构层面按80/10/10的比例划分CATH数据库,结果是16153个结构被分配到训练集,1457个到验证集,1797个到测试集。需要格外注意防止测试集中的信息通过预测结构泄露。使用Gene3D拓扑分类来筛选训练中用于监督的序列,以及作为AlphaFold2预测输入的多序列比对(MSA)。还在CATH测试集的一个较小子集上进行评估,该子集使用Foldseek通过TMscore进一步筛选,以排除与训练集中任何相似的结构。
- 图3:一个UniRef50序列(UniRef50: P07260;PDB编号:1AP8)的AlphaFold预测结果与实验结构对比示例。实验结构以带透明度的粉色显示。预测结果根据pLDDT置信度分数上色,高置信度区域为蓝色。
模型架构
作者研究使用几何向量感知器(GVP)层的模型架构(Learning from protein structure with geometric vector perceptrons),该架构能够学习向量特征的旋转等变变换以及标量特征的旋转不变变换。
作者展示了三种模型架构的结果:(1) 来自Jing等人的GVP-GNN,它目前是逆折叠领域的最先进模型;(2) 宽度和深度增加的GVP-GNN(GVP-GNN-large);(3) 一种混合模型,由一个GVP-GNN结构编码器和其后的通用Transformer组成(GVP-Transformer)。评估中使用的所有模型都训练至收敛。
在逆折叠过程中,预测出的序列应与结构坐标的参考系无关。对于输入坐标的任何旋转与平移变换 T T T,我们希望模型的输出在这些变换下保持不变,即 p ( Y ∣ X ) = p ( Y ∣ T X ) p(Y|X)=p(Y|TX) p(Y∣X)=p(Y∣TX) 。本文研究的GVP-GNN和GVP-Transformer逆折叠模型均具有这种不变性。
GVP-GNN
从Jing等人所述的具有3个编码器层和3个解码器层的GVP-GNN架构开始,采用Jing等人描述的向量门。作为GVP-GNN的输入,蛋白质结构被表示为邻近图,其中每个氨基酸对应图中的一个节点。节点特征由源自二面角的标量节点特征和源自骨干原子相对位置的向量节点特征组合而成,而边特征则反映了相邻氨基酸的相对位置。
在利用预测结构进行训练时,作者发现一个更深且更宽的GVP-GNN表现更佳,它有8个编码器层和8个解码器层(GVP-GNN-large)。
GVP-Transformer
作者使用GVP-GNN编码器层来提取几何特征,随后接入一个通用的自回归编码器-解码器Transformer。在GVP-GNN中,输入特征具有平移不变性,且每一层都具有旋转等变性。作者对来自GVP-GNN的向量特征进行基变换,转换到为每个氨基酸定义的局部参考系中,以得到旋转不变特征。在对比研究中,增加GVP-GNN编码器层的数量可提升整体模型性能,这表明GVP-GNN中的几何推理能力与Transformer层具有互补性。模型规模扩大至拥有4个GVP-GNN编码器层、8个通用Transformer编码器层和8个通用Transformer解码器层,参数达1.42亿的GVP-Transformer模型时,性能得到提升。
训练
结合实验数据和预测数据
在训练过程中,作者在每个训练轮次将实验得到的结构训练集(约1.6万个结构)与AlphaFold2预测训练集的10%随机样本(1200万个结构中的10%)混合,使得实验数据与预测数据的比例达到1:80。对于更大的模型而言,训练时使用较高比例的预测数据有助于避免在规模较小的实验训练集上出现过拟合现象。虽然添加预测数据能够提升模型性能,但仅使用预测数据进行训练会导致性能大幅下降。
对于目标序列中的每个氨基酸,损失的权重是相等的。ESM-IF1会屏蔽掉AlphaFold2置信度分数(pLDDT)低于90的预测输入坐标,这大约占预测坐标的25%。pLDDT置信度分数的可视化可参见图3。这些低置信度区域通常出现在序列的起始和末尾,可能对应于无序区域。作者在每个序列开头添加一个标记,以表明该结构是实验测定的还是预测得到的。对于每个残基,将AlphaFold2给出的pLDDT置信度分数作为一个特征,通过高斯径向基函数进行编码。
在训练过程中,向预测结构添加尺度为0.1埃的高斯噪声,模型性能会稍有提升。这一发现与埃杜诺夫等人的研究结果一致,他们观察到,与最大后验(MAP)预测相比,使用带噪声的合成数据进行反向翻译能提供更强的训练信号。
跨度掩码
为了实现对部分掩码骨干结构的序列设计,作者在训练过程中引入骨干掩码。对独立随机掩码和跨度掩码都进行了实验。在自然语言处理中,跨度掩码相较于随机掩码能提升性能。作者随机选择长度不超过30个氨基酸的连续跨度,直至15%的输入骨干坐标被掩码。
结果
对于固定骨架设计,作者首先在给定所有骨架坐标的序列设计标准设定下进行评估。然后,作者从三个维度使序列设计任务更具挑战性:(1)对坐标引入mask;(2)推广到蛋白质复合物;(3)基于多种构象进行条件设计。此外,作者表明逆折叠模型对于蛋白质复合物稳定性-protein complex stability、结合亲和力-binding affinity和插入效应-insertion effects是有效的零样本预测器。
固定骨架的蛋白质设计
作者从给定蛋白质骨干原子(N、Cα、C)坐标来预测天然蛋白质序列这一任务入手。针对此任务,困惑度以及在留出的天然序列上的序列恢复情况,是两个常用指标。困惑度衡量的是天然序列在预测序列分布中的逆似然性(似然性越高,困惑度越低)。序列恢复(准确率)衡量的是采样得到的序列在每个位置与天然序列匹配的频率。表1使用困惑度和序列恢复指标,对在结构上留出的骨干上的模型进行了比较。
- 表1:固定骨架的序列设计。基于每个残基的困惑度(越低越好;最低困惑度加粗显示)和序列恢复率(越高越好;最高序列恢复率加粗显示)对模型进行比较。Short代表短序列,Single-chain代表单链,All代表所有类型的序列。
注意到,当前最先进的逆折叠模型受到CATH训练集的限制。在CATH数据集上,将当前100万参数的模型(GVP-GNN)扩展到2100万参数(GVP-GNN-large)会导致过拟合,序列恢复率从42.2%降至39.2%。另一方面,当前100万参数规模的模型无法有效利用预测结构:用预测结构训练GVP-GNN会使序列恢复率降至38.6%,且随着训练中预测结构数量的增加,性能会进一步恶化(图6a)。
- 图6:训练数据的对比研究。a.增加预测结构数量的影响。b.训练期间增加预测结构与实验结构混合比例的影响。对于GVP-GNN-large和GVP-Transformer而言,预测结构比例越高,性能提升越明显。c.GVP-GNN和GVP-Transformer的模型规模。
更大规模的模型得益于利用AlphaFold2预测的UniRef50结构进行训练。相较于仅使用实验测定结构进行训练,使用预测结构训练,GVP-GNN-large的序列恢复率从39.2% 提升至50.8%,GVP-Transformer的序列恢复率从38.3% 提升至51.6% 。这种提升也体现在困惑度上。
部分掩码骨架
作者在部分masking骨架上对模型进行评估。虽然训练期间的掩码操作不会显著改变模型在未掩码骨架上的测试性能,但掩码操作确实能让模型对掩码区域的序列进行有实际意义的预测。尽管GVP-GNN-large在短长度掩码上的困惑度较低,但在掩码长度超过5个氨基酸时,其性能会迅速下降至背景分布的困惑度(图4)。相比之下,GVP-Transformer模型即使在较长的掩码区域也能保持一定性能,如果使用跨度掩码而非独立随机掩码进行训练,其性能下降幅度更小(图4)。
- 图4:不同长度掩码坐标区域的困惑度。对于有几个以上token的掩码区域,GVP-GNN架构的困惑度降至背景分布的困惑度,而GVP-Transformer在长掩码跨度上保持一定的准确率。
蛋白质复合物
尽管训练数据仅包含单链蛋白质(前面的实验也都是单链下的子集),但作者发现这些模型能够推广应用于多链蛋白质复合物。通过在各条链之间添加10个mask token来连接它们,以此表示复合物,并在连接时将用于序列设计的目标链置于开头。作者纳入了CATH测试集中长度不超过1000个氨基酸的所有复合物。对于作为蛋白质复合物组成部分的链而言,与仅输入单链坐标相比,当把完整的复合物坐标作为输入时,两种模型的困惑度都有显著改善(表2)。这表明,GVP-GNN和GVP-Transformer都能够利用来自三维结构上相近但序列上相距较远的氨基酸之间的链间信息。
- 表2:在CATH测试集中,当仅给出一条链的骨干坐标(“链”列)以及给出复合物所有骨干坐标(“复合物”列)时,复合物序列设计的性能。两列中的困惑度均针对复合物中的同一条链进行评估。
多种构象
multi-state设计在工程酶和生物传感器领域备受关注。一些蛋白质以多种不同的折叠形式处于平衡状态,而另一些蛋白质在与伴侣分子结合时可能呈现出不同的构象。
对于一个骨架 X X X,逆折叠模型会针对该骨架预测出可能序列 Y Y Y 的条件分布 p ( Y ∣ X ) p(Y | X) p(Y∣X)。为了设计出与状态 A A A 和状态 B B B 都兼容的蛋白质序列,我们希望找到在每个状态的条件分布 p ( Y ∣ A ) p(Y | A) p(Y∣A) 和 p ( Y ∣ B ) p(Y | B) p(Y∣B) 中具有高似然性的序列。我们将这两个条件似然性的几何平均值,作为以与两种状态都兼容的序列为条件的期望分布 p ( Y ∣ A , B ) p(Y | A, B) p(Y∣A,B) 的近似值。
作者在PDBFlex数据集中87个具有多种构象的测试蛋白质上,比较了单态和多态序列设计的性能。在局部柔性残基上,多态设计比单态设计产生更低的序列困惑度(图7)。
-图7:双态设计。对于PDBFlex中结构独立的蛋白质,以两种构象为条件的GVP-Transformer,在局部柔性残基处产生的序列困惑度,比以单一构象为条件时更低。
零样本预测
接下来的结果表明,逆折叠模型在实际设计应用中是有效的突变效应零样本预测器,包括对复合物稳定性、结合亲和力以及插入效应的预测。为了评估某一突变对特定序列的影响,在已知通过实验测定的野生型结构的情况下,依据逆折叠模型,采用突变序列与野生型序列的似然比来进行衡量。由于GVP-GNN和GVP-Transformer均基于自回归解码器,所以二者都能够进行精确的似然估计。然后,将这些似然比得分与在同一组序列上通过实验测定的适应度值进行比较。
似然比计算
对于给定的蛋白质骨架构象 X X X,逆折叠模型会预测出可能序列 Y Y Y 的条件分布 p ( Y ∣ X ) p(Y|X) p(Y∣X)。这里的 X X X 可以是通过实验测定的野生型蛋白质结构,而 Y Y Y 可以是野生型序列 Y w t Y_{wt} Ywt 或者突变序列 Y m u t Y_{mut} Ymut。
根据逆折叠模型,计算野生型序列的似然值 p ( Y w t ∣ X ) p(Y_{wt}|X) p(Ywt∣X),以及突变序列 Y m u t Y_{mut} Ymut 在相同结构 X X X 下的似然值 p ( Y m u t ∣ X ) p(Y_{mut}|X) p(Ymut∣X)。
似然比(Likelihood Ratio,LR)的计算公式为: p ( Y m u t ∣ X ) p ( Y w t ∣ X ) \frac{p(Y_{mut}|X)}{p(Y_{wt}|X)} p(Ywt∣X)p(Ymut∣X)
复合物稳定性
作者利用Atom3D基准测试对模型进行蛋白质复合物界面突变效应的零样本预测评估。该基准测试将SKEMPI数据库中的结合自由能变化作为一个二元分类任务。作者发现,即使不以预测结构作为训练数据,GVP-GNN的序列对数似然值也是蛋白质复合物稳定性变化的有效零样本预测指标(表C.5),其性能与使用迁移学习的最佳监督方法相当。虽然观察到在训练中加入预测结构后困惑度有显著改善(表2),但这并未进一步提升对SKEMPI中单点突变的复合物稳定性预测(表C.5),这表明仅基于单点突变评估模型可能存在局限性。
- 表C.5:SKEMPI测试集上的蛋白质复合物稳定性(单点突变稳定性增加的二元分类)。尽管逆折叠模型仅在单链上进行训练,但它们能够推广应用于蛋白质复合物。与仅将单链作为输入(chain)相比,将完整复合物作为输入(complex)可提升性能。此处将零样本预测与在SKEMPI训练集上训练的全监督及迁移学习方法进行比较。
结合亲和力
虽然SKEMPI数据集每个蛋白质只有一个突变条目,但我们还想评估逆折叠模型是否能够对同一蛋白质的不同突变进行排序,这可能有助于实现结合亲和力的优化,而结合亲和力优化是治疗性设计中的一项重要任务。借助斯塔尔等人(2020年)生成的一个数据集,评估逆折叠模型是否能够预测突变对结合的影响。在该数据集中,针对新冠病毒受体结合域(RBD)的所有单氨基酸替换,都通过实验测量了其与人血管紧张素转化酶2(ACE2)的结合亲和力。鉴于其在界面优化或设计方面的潜在应用,作者重点关注受体结合基序(RBM)内的突变,RBM是RBD中与ACE2直接接触的部分。当给定所有RBD和ACE2的坐标时,表现最佳的逆折叠模型所产生的RBD序列对数似然值与实验测得的结合亲和力的斯皮尔曼相关性为0.69。当不向模型提供ACE2坐标时,观察到相关性较弱,这表明逆折叠模型会利用结合伙伴的结构信息。当对RBM坐标进行掩码处理时(195个残基中有69个被掩码,掩码范围比模型训练期间更长),不再观察到RBD对数似然值与结合亲和力之间存在相关性,这表明模型依赖界面处的结构信息来识别能够保留结合能力的界面设计。通过逆折叠进行的零样本预测优于基于序列的变体效应预测方法,后者利用每个位置上突变氨基酸与野生型氨基酸之间的似然比来预测突变对结合亲和力的影响。
序列插入
通过在插入区域使用掩码坐标标记,逆折叠模型也能够预测插入效应。在腺相关病毒(AAV)衣壳变体上,发现序列对数似然值的相对差异与布莱恩特等人(2021年)通过实验测量得到的插入效应相关。