蛋白结构的“精雕细琢”——在Hermite平台上使用RiD精修蛋白

DP Technology 2021-07-23
蛋白结构预测

蛋白质在生命活动中起着至关重要的作用,其结构与生理功能密切相关。为了获得蛋白质的三维结构,实验生物学家发展了不同的解析结构的方法,目前主要的三大实验手段为X射线晶体学,核磁共振,冷冻电子显微镜。然而蛋白质序列的数量飞速增长,UniProt数据库中机器注释的序列目前已经有2.1亿多,然而PDB数据库中解析出的三维结构的蛋白质数量不到18万,实验方法成本高,周期长。因此给定序列的蛋白质结构预测具有重要意义,同时这也是计算生物学中最有挑战性的问题之一。

计算生物学家发展了不同的蛋白质结构预测(Protein Structure Prediction, PSP)方法。其中一类方法基于模板进行蛋白质结构预测,称为同源建模。尽管精度稳定提升,基于模板的结构预测仍难以在同源序列较少时起作用,且模型精度在蛋白质环区部分显著降低。环区具有高可变性的序列,在基于模板的预测中常处于序列未对齐的片段;在三维结构中处于溶剂可及表面,构象具有巨大柔性;在X射线晶体学中,也由于电子密度低更难解析。环区的建模和结构优化对基于模板的结构预测和蛋白质结构优化都具有重要意义。文献报道的建模方法可分为从头算,基于数据,和二者联合方法。

它们共同面临的挑战是对较长环区的处理,从头算法采样的状态空间指数增长且难以摆放侧链,基于数据的方法数据不足。另一类结构预测方法不基于模板。过去数年间,通过引入多序列比对得到的共进化统计信息,和最新的深度学习方法,人们对无模板蛋白质结构预测取得了重大突破,标志性事件是2020年AlphaFold2在结构预测精度上的巨大提升。从AlphaFold2的预测结构来看,在人类蛋白质组中,有35.7%的残基落在最高精度带内,58.0%的残基被有把握地预测(pLDDT>70)。AlphaFold2 对CASP14中有47% 的蛋白、37% 的100 个残基以下的蛋白,结构预测的结果 GDT-HA打分小于 75。对这些蛋白,仍然需要进一步的结构精修。

深势科技
图1 残基预测精度分布
蛋白结构精修

尽管AlphaFold2带给了我们前所未有的视野,但是还会存在部分区域不准的情况;而蛋白结构精修(Protein Structure Refinement)可以在此基础上帮助我们获得更趋近于生理状态的蛋白结构,即得到更加精确的局域结构信息,这对于酶反应活性位点机理研究,分子对接前蛋白构象准备等工作有重要帮助。

目前,主要的精修(Refinement)手段有两类:一类是通过图神经网络(Graph Neural Network, GNN)1对蛋白编码后生成距离约束,进一步通过蛋白折叠程序(如:Rosetta)折叠得到;另一类是从初始结构出发,借助分子动力学模拟(Molecular Dynamics, MD)对局域相空间采样2,从所得构象系综中找出自由能面上的能量极小点作为生理结构。其中,基于MD方法的蛋白精修方式表现最好3,4

然而,基于MD的蛋白精修方法存在较大的局限性。其主要问题在于:难以对相空间充分采样找到自由能全局极小点。如果蛋白的两种构象之间存在较大的势垒,那么普通的MD过程需要经过很长时间才能越过这一势垒,这导致了该方法虽理论上可行,但实际上对蛋白的优化有限。一般而言,为了提高采样效率,我们可以通过基于温度或基于反应坐标的增强抽样(Enhanced Sampling)方法实现。这两种方法也存在弊端:首先,蛋白体系的动力学模拟的温度不能过高,否则蛋白会“变性”;其次,传统的基于反应坐标的增强抽样方法(如:Metadynamics)只能选择不到三个反应坐标(Collective Variable, CV),而如果我们一次性选取多个反应坐标,那么这些抽样的方法将消耗大量的计算资源。为了能够处理高维CV的增强抽样问题,我们提出了基于神经网络的增强动力学方法(Reinforced Dynamics, RiD)5

增强动力学

RiD能够选取任意个数的CV作为反应坐标,并对所选CV增强,这一优势得益于神经网络强有力的表达能力。RiD将训练一个神经网络,网络的输入为所选的反应坐标,网络的输出为对应的自由能, 即得到函数:

深势科技
功能组件,

RiD主要由三个组件组成,并依次迭代进行,分别是:增强抽样模块、限制性动力学&自由能计算模块、网络训练模块。

增强抽样模块

这一模块将神经网络给出的自由能作为偏置势(Biased Potential)加入采样过程,进行增强抽样,并得到增强采样后的动力学轨迹。这一部分可以并行多条轨迹同时采样。

限制性动力学&自由能计算模块

在这一模块中,RiD将分析上一个模块中采样得到的轨迹中的构象,从中选出没有被训练过的构象和反应坐标,计算这些构象对应的CV和平均力(Mean Force),作为训练数据进一步训练网络。平均力的计算过程通过限制性动力学(Restrained MD)实现。

网络训练模块

我们将根据上一个模块中得到的训练数据(同时整合之前迭代的训练数据)训练神经网络。这些网络为全连接的结构,在当前的Web版本中,我们将同时使用ResNet的结构。训练好神经网络后,RiD完成本次迭代,并使用训练好的网络进行下一次迭代,直到完成指定的迭代次数。

深势科技
图2 RiD工作流

在网络训练步骤中,RiD会训练多个(一般为4个)神经网络模型,这些网络各自独立初始化,但使用相同的训练集训练。RiD将采用不确定性指示函数(Uncertainty Indicator),即4个模型的输出结果的标准差(也叫模型偏差,Model Deviation),来指导如何添加偏置势和如何筛选训练集。RiD内部将会对不确定性的阈值进行自适应的调整,以确保每次采样都能探索到新的相空间区域。

计算方法
使用RiD作为增强抽样方法对蛋白的局域相空间抽样,克服了传统增强抽样方法无法选取高维反应坐标的问题。我们将这一功能集成在Hermite中(Uni-ProteinRefine)。在Hermite Uni-ProteinRefine中,我们将使用RiD方法辅助蛋白结构精修。经过测试,我们将CV选取限定在二面角上,用户可以自由选取主链上任意个数的标准氨基酸的二面角作为反应坐标。采样结束后,我们将从所有的构象中用RWplus6和EvoEF27两种打分函数对构象进行筛选,从中选出一批分数较高的构象系综,计算其平均结构后,再进行能量优化作为最终输出结构。

RiD的采样过程可以并行进行,因此用户可以设置不同个数的轨迹分别进行采样,每一段轨迹我们称为一个walker。同时,用户可以设置RiD的总迭代(Iteration)次数。

在结构精修过程中,为了防止蛋白结构过分偏离初始结构,我们需要依照初始结构在采样过程中为蛋白的Alpha C 原子施加限制势。为了尽量避免限制势对采样过程的影响,采用平底谐振子势(Flat Bottom Harmonic Potential)8对蛋白进行约束往往能够给出较好的结果。为此,我们额外提供了一个选项为bottom_width,其中bottom_width即为平底势的宽度。需要注意的是:势底越宽,采样空间越大,采样越自由,但结果也容易偏离初始构象;势底越窄,采样空间受限,但结果会与初始结构贴近。
结果分析
我们使用了CASP13竞赛中三个例子(R0974s1, R0986s1 和 R1002-D2)测试了RiD在蛋白结构精修的作用,并与目前精修结果排名最好的Feig课题组的结果进行对比9

对于三个测试用例,我们分别选取了136,182,116个二面角作为CV,采用了8条轨迹并行采样,使用了含有4层隐藏层的DNN网络,每层含有1200个节点,经过16个迭代后,得到精修后的结果。三个结构的初始GDT-HA分数为65.6,59.2和72.9,经过RiD精修后,它们的分数均达到78以上,其中最好分数分别为92.4,82.6和79.7,结构优化提升显著。从图3a中可以看出,对比目前表现最好的FEIGLAB的工作,RiD在三个测试例子中达到了80.5的平均GDT-HA分数,相较于初始结构的平均65.9分表现极佳,且显著高于Feig课题组得到的平均75.9分。这一分数表明:相对于经典的分子动力学模拟,RiD的确探索到了更大的构象空间。另外,通过计算三个例子的测试标准差得到的Error Bar可以看出(图3b),RiD的模拟结果方差小,方法具有较好的鲁棒性。
深势科技
图3 RiD在蛋白结构精修中的结果(以CASP13中的R0974s1, R0986s1 和 R1002-D2为例)

进一步分析MD过程的轨迹,我们将增强抽样得到的轨迹做tICA分析,将其投影至两个主成分构成的二维子空间中(图4d、4f、4h),可以清楚的看出经典轨迹探索的相空间与RiD方法所探索的相空间对比。其中浅蓝色叉号代表天然结构(Native Structure);黑色叉号代表初始结构;虚线为经典MD所探索的相空间区域;填充色为RiD方法所探索的相空间区域。可以看出,RiD给出了更大的相空间区域,确实可以从一个初始的状态,采样到天然的状态。
深势科技
图4 RiD在蛋白结构精修中GDT-HA分数随时间的变化以及采样空间

尤其是对于R1002-D2一例(图4h,图5),经典MD轨迹明显陷入了局域势阱中而根本无法触及天然结构,RiD方法则成功越过势垒并采集到了天然结构。这充分体现了RiD在大体系和高维空间中的采样能力。
深势科技
图5 RiD在蛋白R1002-D2中的结果 白色为天然结构,蓝色为预测的结构,红色是通过RiD精修过的结构

了解更多RiD的细节,请阅读参考材料5和9。
使用范例

我们选用CASP13中的一个例子R1002-D2蛋白作为本教程的用例蛋白 https://predictioncenter.org/casp13/target.cgi?id=206&view=all 。 登录Hermite,进入项目管理界面,创建一个名为"ProteinRefine"的新项目。单击项目名进入。

深势科技
图6 Hermite项目管理界面

将鼠标悬停在页面左上角Function处会出现Hermite的功能,我们选取Uni-ProteinRefine功能。在弹出的对话窗中单击 From file,然后从本地上传R1002-D2蛋白的PDF文件,单击Next。

深势科技
图7 Hermite Uni-ProteinRefine功能

接下来进入Uni-ProteinRefine的参数设置界面。Uni-ProteinRefine支持用户选取蛋白质主链上任意氨基酸的二面角作为CV,我们可以点击Residue左侧的详情来查看这个蛋白中包含的残基序号和种类。我们提供了两种选取蛋白质主链上任意氨基酸残基的方式:通过输入残基在PDB文件中的序号选取或者通过选定残基范围选取。这里我们采用第二种方式,单击Select Residues by Range 旁边的开关,使之成为ON状态。由于R1002-D2蛋白有59个残基,均为标准氨基酸,我们选取主链上所有的二面角作为CV,在对话框中分别输入60和118,单击Add。下面将显示你添加的残基的信息,拖动滑条可以查看这些氨基酸。

深势科技
图8 Uni-ProteinRefine中选取氨基酸残基

确定无误后,向下滑动滑条,我们将看到剩余的三个配制窗口:

我们可以设置 Number of walkers 为 4,这表明我们将有4条独立的轨迹进行采样(参看RiD原理模块一)。Walker个数越多,采样的空间越大,但相应的计算耗费成倍増长。我们建议使用4~8个walker,过多walker的结果可能趋同。

我们可以设置Number of Iteration 为 4,RiD将会运行4个迭代。我们建议进行4~8次迭代,迭代较少时神经网络训练不充分,采样过程不能得到有效增强(迭代数必须大于1,否则神经网络将不会被训练);迭代次数过多则会消耗过多计算资源且结果不会线性变优。

设置势底宽度为0.6 nm。我们推荐0.4~0.6 nm尺度下的限制,过强的限制会显著影响采样过程;过弱的限制会使结果显著偏离初始结构。

深势科技
图9 Uni-ProteinRefine中设置强化动力学参数参数

设置完成后,给该任务取名为"ProteinRefine",单击Submit,单击右上角的Jobs可以查看已提交的任务。本用例蛋白存在59个氨基酸残基,按上述参数设定,计算将在约四小时后完成。任务完成后,可以点击Jobs看到状态显示为:Finished,随后可以点击show按钮展示结果,也可以点击download icon进行结果下载。

深势科技
图10 Hermite任务界面

点击show后,右侧结果栏中展示了5个可能的结果,它们分别是平均了不同个数的构象得到的平均结构,并做局部结构优化得到的,我们分别使用RWplus和EvoEF2对该五个结构的打分。RWplus为基于统计数据的打分函数,EvoEF2则结合了力场参数,理论上,两者均是数值越低结构越好。

点击结果栏中的proteinID,我们可以把结果加载到前端叠合分析。下图中,玫红色是原始结构,灰色是优化后得到的结构。用户可以选择合适的结构作为最终结果。

深势科技
图11 Uni-ProteinRefine结果展示
总结起来,RiD的优势在于:能够借助神经网络处理高维CV的增强抽样问题,同时对多个CV添加偏置势,克服了传统增强抽样方法的不足,大大提高了大体系的采样效率;通过大量选取CV,降低了CV选择个数较少时动力学模拟对CV选取的敏感性;能够有效适用于没有先验知识或物理直觉的大分子体系。

请通过链接填写问卷:https://wj.qq.com/s2/8730755/dae3 我们会分批开放名额,并以微信和邮件的方式通知您进入内测
References:
1. Xiaoyang Jing, Jinbo Xu, Fast and Effective Protein Model Refinement Using Deep Graph Neural Networks, Nature Computational Science, 2021, ASAP. DOI: 10.1038/s43588-021-00098-9.
2. Lim Heo, Collin F Arbour, and Michael Feig. Driven to near-experimental accuracy by refinement via molecular dynamics simulations. Proteins: Structure, Function, and Bioinformatics,87(12):1263–1275, 2019.
3. https://predictioncenter.org/casp13/zscores_final_refine.cgi
4. https://predictioncenter.org/casp14/zscores_final_refine.cgi
5. Zhang L , Wang H , Weinan E , . Reinforced dynamics of large atomic and molecular systems[J]. 2017.
6. Zhang J, Zhang Y (2010) A Novel Side-Chain Orientation Dependent Potential Derived from Random-Walk Reference State for Protein Fold Selection and Structure Prediction. PLoS ONE 5(10): e15386.
7. Xiaoqiang Huang, Robin Pearce, Yang Zhang. EvoEF2: accurate and fast energy function for computational protein design. Bioinformatics (2020) 36:1135-1142
8.https://manual.gromacs.org/documentation/2019/reference-manual/functions/restraints.html
9. Wang, D., Zhang, L., Wang, H., E, W.. Efficient sampling of high-dimensional free energy landscapes using adaptive reinforced dynamics[J]. arXiv preprint arXiv:2104.01620, 2021.