蛋白与蛋白的对接计算教程

1. 前言

本教程讲述如何使用ZDOCK分子对接软件、殷赋云平台分子动力学方案及相关小工具，实现蛋白与蛋白的对接计算。

本次示例中，受体蛋白与配体蛋白分别为CDK2（PDB ID : 3PY1）和Cyclin-A2（PDB ID : 1VIN）。该体系并不清楚具体的氨基酸残基位点，所以采用盲对接方式确定结合位点。首先采用ZDOCK进行刚性对接，挑选若干复合物进行分子动力学模拟，然后根据结合自由能选取最好的结果进行结合模式分析。

• 相对于明确位点的对接方式来说，盲对接的花费较大。因此，建议尽量通过文献报道、实验数据确定作用残基后再进行对接，选择最好的复合物进行下一步计算；
• 以下网站（软件）也可进行蛋白-蛋白对接：HDOCK、pyDockWEB、ClusPro2.0、
  HADDOCK。
• 本教程同样适用于蛋白与多肽对接。

2. 流程图

graph LR
A(准备结构) --> B(刚性对接)
B(刚性对接) --> C(选取复合物)
C(选取复合物) --> D(分子动力学模拟)
D(分子动力学模拟) --> E(结合自由能计算)

3. 步骤

3.1 准备蛋白质三维结构

1. 打开殷赋云计算平台（https://cloud.yinfotek.com/) 【处理pdb结构】小工具。

   

2. 输入PDB ID 3PY1（或上传PDB文件），删除多余结构，保留用于对接的蛋白肽链。

   处理原则：删除不参与对接的肽链、杂原子及水分子。

   ​

3. 处理后，将receptor.pdb改名为3py1.pdb。

4. 重复上述步骤，处理配体蛋白 1VIN，文件命名为1vin.pdb。

   注意，ZDOCK对接会将受体和配体合并成复合物。为便于区分，两者应采用不同链名。例如，本例的受体和配体蛋白均为A链，可将配体蛋白改为B链。修改方法：使用文本编辑器Sublime或NotePad++打开配体蛋白1vin.pdb，列选链名A，改为B，或采用字符替换方法，将A替换为B（注意前后各有一个空格，避免替换掉单词中的A）。

3.2 刚性对接

1. 打开ZDOCK在线服务网页（https://zdock.umassmed.edu/)进行刚性对接。由于不清楚具体的作用残基（“盲对接”），所以勾选Skip residue selection。

   除了zdock，还有很多可以进行刚性对接的网站或软件，请用户自行搜索或查找相关文献。

   这里选择盲对接，所以勾选"Skip residue selection"。若知道具体结合的氨基酸残基，请不要勾选。

   更多zdock网站的说明，请看Help说明，我们不再解答该网站的相关内容。

   

2. 一段时间后，计算结果将发送至预留的邮箱，通常有10个复合物结果。这里，我们选择其中5个进行下一步计算。

   

3.3 分子动力学模拟

3.3.1 准备复合物三维结构

这里以complex1.pdb为例，介绍分子动力学模拟及结合自由能计算的过程，其他复合物同理。

1. 打开殷赋云平台【处理PDB结构（进阶版）】小工具。

   

2. 上传复合物结构，点击【下一步】。

   

3. 点击【生成文件】，下载文件。

   • 有关PDB修复的内容，详见教程《处理PDB结构（进阶版）》。
   • 默认勾选对残基重新编号，renumber.csv记录了残基编号的对应关系，方便后续分析，切记保存好。
   • 这里虽然显示receptor.pdb，但实际上是复合物结构。

   

3.3.2 准备Amber文件

1. 打开殷赋云平台【准备Amber文件】小工具。

   

2. 上传前面准备好的receptor.pdb文件，点击【生成】，下载文件。

• 无需上传参数文件；

• amber_files.zip中的residue_numbers.csv文件包含了蛋白、离子和水的残基编号，方便分子动力学模拟时设置约束条件。

​

3.3.3 开展动力学模拟

1. 打开殷赋云平台【分子动力学（Amber20）】方案，创建任务，进入任务页面。

   

2. 上传准备好的拓扑文件和坐标文件，设置模拟环境。

   

3. 设置计算步骤，添加约束条件，点击【提交】。

   

简单解释上图各计算步骤目的：

1. Minimization：约束蛋白的原子，优化水分子、盐离子；

2. Minimization：约束蛋白骨架原子，优化其他原子；

3. Minimization：不加约束，优化整个体系；

4. Heat：将体系（从0）升温至300 K，为避免剧烈运动破坏体系，与第1步Minimization
   一样设置约束；

5. Equilibration：与第2步Minimization一样设置约束，使用NPT系综来控制压强，让体
   系调整自身密度；

6. Equilibration：体系密度已达平衡，不再有大变动，可使用NVT系综控制体积，并放开
   所有约束，让体系继续平衡；

• 需要预先估计该体系总费用的用户，请参考以下指导完成分步操作。

  1. 完成第1到第6步操作后，先不设置第7步，点击【提交】，等计算完毕，到 用户中心-费用账单-账单明细查看费用，估算出跑1ns的费用。

  2. 如需完成后续模拟，可以进入我的项目，点击查看，从文件列表中下载prmtop（拓扑）和rst（坐标）文件（下拉文件列表至末尾，选取计算步骤6生成的rst文件）作为后续模拟的起点。

     

  3. 新建分子动力学任务，上传 prmtop 和 rst 这两个文件，参考步骤7Production的时长分段建议设置production时长。（步骤1-6的信息已在prmtop和rst文件中，无需设置重跑)

• 不需要估计该体系总费用的用户，可继续第7步骤。

7. Production：与上一步Equilibration的条件相同，从此开始长时间采样。

   • 通常所说跑多少ns的动力学，是指生产（Production）阶段的总时长；这里设置一个20 ns。

   • 建议将长时间模拟分成若干小段，例如，100 ns的动力学可分成连续的5个20ns Production。这样好处很多：可避免意外事故发生时毁掉整段轨迹（只要从最后一段完整的轨迹接着算），可方便后期分析时选择理想时段而不用进行轨迹分割处理，可了解进度并在合适时候提前终止任务……

     

约束条件中的残基编号可从residue_numbers.csv文件中查看。比如，残基编号总是从1开始，该体系最后一个氨基酸残基编号为551，那么复合物残基范围为1-551。

第1步Minimization和第4步Heat采用以下约束条件：

第2步Minimization和第5步Equilibration采用以下约束条件：

3.3.4 分析结果

1. 待任务完成，点击【查看】，进入分析页面。查看热力学性质。

在平衡（Equilibration）或生产（Production）阶段，能量、温度、压强、体积等曲线应当平稳 。详情见《分子动力学模拟（Amber 20）》

2. 查看生产（production）阶段的RMSD。该曲线在某个高度附近平稳振动，振幅大小在2Å以内，表明体系已达到平衡，可以进一步分析。

3. 氢键。仅含有机小分子的体系才有氢键分析数据，纯生物大分子体系无此数据。若有分析需要，请使用《分析轨迹）》方案。

4. 分子结构。可以下载各个阶段的结构文件，使用PyMOL等软件查看结构。

3.4 开展结合自由能计算

1. 打开殷赋云平台【结合自由能计算】方案，创建任务，进入任务页面。

注意：体系平衡后才能进一步计算结合自由能，否则，误差可能会很大。

2. 从云盘文件选择拓扑文件及轨迹文件，并设置参数，然后【提交】任务。

   

• 采样范围：总共20000帧（20 ns），采样间隔为100，则采样数为20000 ÷ 100 = 200帧。可根据体系大小进行设置，一般采样200-500帧即可。

• 结构划分：查看【处理pdb结构（进阶版）】生成的renumber.csv文件，了解受体和配体蛋白的残基编号范围。如下图，受体为A链，对应残基编号1-299，配体为B链，对应残基编号300-551。

  

3.5 分析结果

1. 待任务完成，点击【查看】，进入分析页面，查看结合自由能，并下载mmpbsa.xlsx文件（后面分析用到）。

• MM/GB(PB)SA方法是指采用MM方法计算$\Delta G_{\mathrm{gas}}$（= $\Delta G_{\mathrm{vdw}}$ + $\Delta G_{\mathrm{ele}}$），GB和PB方法计算$\Delta G_{\mathrm{solv}}$（= $\Delta G_{\mathrm{polar}}$ + $\Delta G_{\mathrm{nonpolar}}$）。因此，GB和PB列的$\Delta G_{\mathrm{vdw}}$、$\Delta G_{\mathrm{ele}}$与$\Delta G_{\mathrm{gas}}$数值相同，而$\Delta G_{\mathrm{polar}}$、$\Delta G_{\mathrm{nonpolar}}$与$\Delta G_{\mathrm{solv}}$不同。

• 总结合自由能（$\Delta G_{\mathrm{total}}$）数值越小，结合力越强，各能量项数值越小，越有利于结合。

2. 同理，其他复合物也分别开展分子动力学模拟和结合自由能计算，最后将总结合自由能（$\Delta G_{\mathrm{total}}$）数据整理如下：

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

由此看出，复合物1的结合力最强，所以我们取其分子动力学模拟的最后一帧进行结合模式分析。

若GB和PB的结论不一致，可分析各能量项，排除$\Delta G_{\mathrm{vdw}} $、$\Delta G_{\mathrm{ele}}$与$\Delta G_{\mathrm{gas}} $不利的复合物体系。

3.5 结合模式分析

1. 进入分子动力学模拟结果页面，下载分子结构。

这里下载7.prod-dry.pdb文件。

​

2. 用PyMOL或自己熟悉的分子图形软件打开分析两个蛋白之间的相互作用。

分子动力学模拟会更改原子编号，请先改回原来编号再进行分析，可查看上面准备复合物三维结构时保留的renumber.csv文件。

从mmpbsa.xlsx文件的gb_total或pb_total页可以看出，配体的GLU387与受体的ARG151和ARG158都有非常强的静电力作用。将它们显示在PyMOL上，可见它们之间刚好形成盐桥（图中黄色虚线表示），盐桥距离分别为3.4 Å和3.8 Å。同理，再分析其他氨基酸残基。