【殷赋云教程】蛋白与多肽的对接计算教程
1. 前言
本教程讲述如何使用CABS-dock网站、殷赋云平台分子动力学方案及相关小工具,实现蛋白与多肽的对
接计算。
本次示例中,受体蛋白分别为3BP1(PDB ID : 1ABO),配体多肽序列为APTMPPPLPP。该体系并不
清楚具体的氨基酸残基位点,所以采用盲对接方式确定结合位点。首先采用CABS-dock进行刚性对接,
挑选若干复合物进行分子动力学模拟,然后根据结合自由能选取最好的结果进行结合模式分析。
- 相对于明确位点的对接方式来说,盲对接的花费较大。因此,建议尽量通过文献报道、实验数据确定作用残基后再进行对接,选择最好的复合物进行下一步计算;
- 以下网站(软件)也可进行蛋白多肽对接:MDockPeP、HPEPDOCK、pepATTRACT、
FlexPepDock、GalaxyPepDock、HADDOCK peptide docking、DINC和AutoDock
CrankPep (ADCP)
2. 流程图
3. 步骤
3.1 准备蛋白质三维结构
-
打开殷赋云计算平台(殷赋云 - 科研计算快易准)【处理pdb结构】小工具。
-
输入PDB ID
1ABO
(或上传PDB文件),删除多余结构,保留用于对接的蛋白肽链。处理原则:删除不参与对接的肽链、杂原子及水分子。
-
处理后,将
receptor.pdb
改名为1ABO.pdb
。
3.2 刚性对接
-
打开CABS-dock在线服务网页(CABS-dock: server for protein-peptide docking)进行刚性对接。Optional选项,请根据实际情况填写。
更多CABS-dock网站的说明,请看Help说明,我们不再解答该网站的相关内容。
-
一段时间后,网站就会得到多个复合物结果。这里,我们选择其中5个进行下一步计算。
3.3 分子动力学模拟
3.3.1 准备复合物三维结构
这里以complex1.pdb
为例,介绍分子动力学模拟及结合自由能计算的过程,其他复合物同理。
-
打开殷赋云平台【处理PDB结构(进阶版)】小工具。
-
上传复合物结构,点击【下一步】。
-
点击【生成文件】,下载文件。
- 有关PDB修复的内容,详见教程《处理PDB结构(进阶版)》
- 默认勾选
对残基重新编号
,renumber.csv
记录了残基编号的对应关系,方便后续分析,切记保存好。 - 这里虽然显示
receptor.pdb
,但实际上是复合物结构。
3.3.2 准备Amber文件
-
打开殷赋云平台【准备Amber文件】小工具。
-
上传前面准备好的
receptor.pdb
文件,点击【生成】,下载文件。
无需上传参数文件;
amber_files.zip
中的residue_numbers.csv
文件包含了蛋白、离子和水的残基编号,方便分子动力学模拟时设置约束条件。
3.3.3 开展动力学模拟
-
打开殷赋云平台【分子动力学(Amber20)】方案,创建任务,进入任务页面。
-
上传准备好的拓扑文件和坐标文件,设置模拟环境。
-
设置计算步骤,添加约束条件,点击【提交】。
简单解释下图各计算步骤目的:
Minimization:约束蛋白的原子,优化水分子、盐离子;
Minimization:约束蛋白骨架原子,优化其他原子;
Minimization:不加约束,优化整个体系;
Heat:将体系(从0)升温至300 K,为避免剧烈运动破坏体系,与第1步Minimization
一样设置约束;Equilibration:与第2步Minimization一样设置约束,使用NPT系综来控制压强,让体
系调整自身密度;Equilibration:体系密度已达平衡,不再有大变动,可使用NVT系综控制体积,并放开
所有约束,让体系继续平衡;
需要预先估计该体系总费用的用户,请参考以下指导完成分步操作。
完成第1到第6步操作后,先不设置第7步,点击【提交】,等计算完毕,到
用户中心-费用账单-账单明细
查看费用,估算出跑1ns的费用。如需完成后续模拟,可以进入
我的项目
,点击查看
,从文件列表
中下载prmtop(拓扑)和rst(坐标)文件(下拉文件列表至末尾,选取计算步骤6生成的rst文件)作为后续模拟的起点。新建分子动力学任务,上传 prmtop 和 rst 这两个文件,参考步骤7
Production
的时长分段建议设置production时长。(步骤1-6的信息已在prmtop和rst文件中,无需设置重跑)不需要估计该体系总费用的用户,可继续第7步骤。
Production:与上一步Equilibration的条件相同,从此开始长时间采样。
通常所说跑多少ns的动力学,是指生产(Production)阶段的总时长;这里设置一个20 ns。
建议将长时间模拟分成若干小段,例如,100 ns的动力学可分成连续的5个20ns Production。这样好处很多:可避免意外事故发生时毁掉整段轨迹(只要从最后一段完整的轨迹接着算),可方便后期分析时选择理想时段而不用进行轨迹分割处理,可了解进度并在合适时候提前终止任务……
约束条件中的残基编号可从
residue_numbers.csv
文件中查看。比如,残基编号总是从1开始,该体系最后一个氨基酸残基编号为68,那么复合物残基范围为1-68。
第1步Minimization和第4步Heat采用以下约束条件:
第2步Minimization和第5步Equilibration采用以下约束条件:
3.3.4 分析结果
- 待任务完成,点击【查看】,进入分析页面。查看热力学性质。
在平衡(Equilibration)或生产(Production)阶段,能量、温度、压强、体积等曲线应当平稳 。详情见《分子动力学模拟(Amber 20)》
- 查看生产(production)阶段的RMSD。该曲线在某个高度附近平稳振动,振幅大小在2Å以内,表明体系已达到平衡,可以进一步分析。
-
氢键。仅含有机小分子的体系才有氢键分析数据,纯生物大分子体系无此数据。若有分析需要,请使用《分析轨迹》方案。
-
分子结构。可以下载各个阶段的结构文件,使用PyMOL等软件查看结构。
3.4 开展结合自由能计算
- 打开殷赋云平台【结合自由能计算】方案,创建任务,进入任务页面。
注意:体系平衡后才能进一步计算结合自由能,否则,误差可能会很大。
-
从云盘文件选择拓扑文件及轨迹文件,并设置参数,然后【提交】任务。
采样范围:总共20000帧(20 ns),采样间隔为100,则采样数为20000 ÷ 100 = 200帧。可根据体系大小进行设置,一般采样200-500帧即可。
结构划分:查看【处理pdb结构(进阶版)】生成的
renumber.csv
文件,了解受体和配体蛋白的残基编号范围。如下图,受体为A链,对应残基编号1-299,配体为B链,对应残基编号300-551。
3.5 分析结果
- 待任务完成,点击【查看】,进入分析页面,查看结合自由能,并下载
mmpbsa.xlsx
文件(后面分析用到)。
MM/GB(PB)SA方法是指采用MM方法计算ΔGgas(= ΔGvdw +ΔGele),GB和PB方法计算ΔGsolv(= ΔGpolar + ΔGnonpolar)。因此,GB和PB列的ΔGvdw、ΔGele与ΔGgas数值相同,而ΔGpolar、ΔGnonpolar与ΔGsolv不同。
总结合自由能(ΔGtotal)数值越小,结合力越强,各能量项数值越小,越有利于结合。
- 同理,其他复合物也分别开展分子动力学模拟和结合自由能计算,最后将总结合自由能(ΔGtotal)数据整理如下:
Complex | Generalized Born (GB) | Poisson-Boltzmann (PB) |
---|---|---|
1 | -98.5613±4.4047 | -102.4870±5.9701 |
2 | -66.3568±2.1541 | -56.9621±2.3619 |
3 | -75.2548±3.2224 | -70.6548±3.0129 |
4 | -50.3648±5.9264 | -55.6325±4.6387 |
5 | -46.8004±7.3251 | -26.7989±7.9052 |
由此看出,复合物2的结合力最强,所以我们取其分子动力学模拟的最后一帧进行结合模式分析。
若GB和PB的结论不一致,可分析各能量项,排除$\Delta G_{\mathrm{vdw}} 、、\Delta G_{\mathrm{ele}}与与\Delta G_{\mathrm{gas}} $不利的复合物体系。
3.6 结合模式分析
- 进入分子动力学模拟结果页面,下载分子结构。
这里下载
7.prod-dry.pdb
文件。
- 用PyMOL或自己熟悉的分子图形软件打开分析复合物之间的相互作用。
分子动力学模拟会更改原子编号,请先改回原来编号再进行分析,可查看上面准备复合物三维结构时保留的
renumber.csv
文件。
从mmpbsa.xlsx
文件的gb_total
或pb_total
页可以看出,配体的LEU与PRO65与受体的TYR52和ASN51都有较强的静电力作用。将它们显示在PyMOL上,可见它们之间形成氢键(图中绿色虚线表示),氢键距离分别为2.7Å和2.9 Å。同理,再分析其他氨基酸残基。