蛋白与蛋白的对接计算教程
1. 前言
本教程旨在介绍如何运用ZDOCK分子对接软件、殷赋云分子动力学平台及相关辅助工具完成蛋白间相互作用的计算任务。
在本例中，在特定的研究案例下（PDB ID：3PY1）下确定了受体蛋白为CDK2，在另一特定条件下（PDB ID：1VIN）下确定配体蛋白为Cyclin-A2。该研究并未明确具体确定哪些氨基酸残基参与了相互作用过程因而采用了基于盲目对齐的结合位点确定方法为此研究建立了一个初步的基础模型随后通过ZDOCK软件实现了刚性对齐将多种可能存在的复合物体系进行了筛选在此基础上利用分子动力学模拟方法从多个候选对象中选择了具有典型特征的结果用于构建初始模型接着通过计算各复合物体系的结合自由能值来筛选出最优的匹配模式并通过热力学稳定性分析进一步验证其可行性
- 相较于明确位点的对接方式而言, 盲对接操作所消耗的时间较长. 因此, 建议尽可能通过文献综述或实验数据确定作用残基后再进行对接, 最终选择效果最优的复合物进行后续计算;
- 以下网站（软件）也可用于蛋白质间的相互作用研究: HDOCK、pyDockWEB、ClusPro2.0以及 HADDOCK.
- 本教程同样也适用于蛋白质与多肽间的相互作用研究.
2. 流程图
3. 步骤
3.1 准备蛋白质三维结构
打开殷赋云计算平台[殷赋云 - 科研计算快易准
处理pdb结构
处理pdb结构
用于处理pdb文件的模块
输入PDB ID**3PY1**（或上传PDB文件），删除多余结构，保留用于对接的蛋白肽链。
处理原则：删除不参与对接的肽链、杂原子及水分子。
处理后，将**receptor.pdb** 改名为**3py1.pdb** 。
重复上述步骤，处理配体蛋白 1VIN ，文件命名为**1vin.pdb** 。
请注意以下事项：当使用ZDOCK进行分子对接时，默认情况下受体与配体会被整合成同一个复合物结构。为了在后续分析中能够明确区分两者，在构建模型时应确保受体与配体分别位于不同的多肽链上：建议选择具有相同编号的两条肽链（如本例中受体与配体均为1号链），并将其中一条的配体分子编号设为2号链即可实现分离效果（见下文详细操作步骤）。具体操作如下：首先请调出目标PDB文件**1vin.pdb** ，在其编辑界面中选定第一条多肽链（即编号为1），然后将其替换为第二条多肽链（编号为2）。此操作需特别注意前后各留一个空格位置以避免影响后续运算结果的准确性。
3.2 刚性对接
访问ZDOCK在线服务页面[ZDOCKServer:全自动蛋白对接服务器]
该蛋白 docking 服务器是一个自动化的工具，在医学研究中具有重要应用价值。为了实现对象间的刚性对齐操作（即固定相对位置关系），系统会默认启用这一功能选项。然而，在某些特定情况下（如缺乏对参与对接氨基酸残基的具体信息——即‘盲对接’），建议将此参数设置为‘跳过残基选择’
除了zdock之外, 还有多个具备刚性对接功能的网站或软件, 请用户自行查找相关文献资料。
这里选择盲对接, 则勾选"Skip residue selection"选项框。若知道具体结合部位的氨基酸残基, 则不建议勾选此选项。
关于zdock网站的详细说明, 请参考Help文档.关于该网站的具体问题, 则不再提供解答。
待会儿将向预留的邮箱地址发送计算结果。一般情况下会有十个化合物产物。在此基础上选择五个化合物产物进行后续计算。
3.3 分子动力学模拟
3.3.1 准备复合物三维结构
例如，在本研究中采用**complex1.pdb**这一实例来阐述其分子动力学模拟及其与自由能计算方法的结合过程，并以此类推。
打开殷赋云平台【处理PDB结构（进阶版）】小工具。
上传复合物结构，点击【下一步】。
点击【生成文件】，下载文件。
关于PDB修复的内容，请参阅教程《处理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系综控制体积，并放开
所有约束，让体系继续平衡；
提交
提交
提交
3.创建分子动力学任务后，请上传指定的prmtop和rst文件，并遵循step 7关于生产时段划分的建议来设置生产时间长度。（其中包含step 1至step 6信息的内容已包含在prmtop和rst文件中）
未涉及对体系总费用进行估算的用户可以直接进入第7步操作。
生产：与平衡化阶段相同条件开始后进行长时间采样。
*一般而言，在进行分子动力学模拟时所指的时间步长（即"ns"），特指生产阶段所需的总时长；在此案例中设定为20 ns。
*建议将长时间运行任务划分为多个连续的小段来处理。例如，在进行100 ns的时间步长模拟时可将其分割为5个连续的20 ns生产阶段（Production）。这样具有诸多优点：当发生意外情况导致轨迹中断时无需担心会破坏整个轨迹（只要从最后一段完整的轨迹接着计算），在后续分析时可以选择理想的运行时段而不必进行轨迹划分处理，并能实时了解进展程度以便适时终止任务……
约束条件中的限制项可以通过访问**
residue_numbers.csv** 文件来确定。例如，在这种情况下，默认情况下所有氨基酸的限制项起始位置通常为第1位。因此在这种情况下，
复合物的整体限制范围通常是从第1位到最大的氨基酸限制项的位置（即第549位）。
第1步Minimization和第4步Heat采用以下约束条件：
第2步Minimization和第5步Equilibration采用以下约束条件：
3.3.4 分析结果
- 待任务完成，点击【查看】，进入分析页面。查看热力学性质。
在平衡或生产环节中所涉及的参数曲线应保持稳定状态。详情请参考《分子动力学模拟（Amber 20）》
热力学分析
考察生产（production）阶段的RMSD值。观察该曲线在特定高度附近呈现稳定波动状态（stable oscillation），其振幅值不超过2Å（angstrom），这表明系统已达到平衡状态（steady state），后续分析工作可据此开展
RMSD
氢作用。
仅含有机小分子类别的体系才能进行氢作用的定量分析。
而纯生物大分子体系不具备生成相关数据的能力。
如需进一步的数据支持，请参考《分析轨迹）》一文。
分子结构。可以下载各个阶段的结构文件，使用PyMOL等软件查看结构。
分子结构
3.4 开展结合自由能计算
- 打开殷赋云平台【结合自由能计算】方案，创建任务，进入任务页面。
注意：体系平衡后才能进一步计算结合自由能，否则，误差可能会很大。
从云盘文件选择拓扑文件及轨迹文件，并设置参数，然后【提交】任务。
处理pdb结构（进阶版）
采样范围：全部共计有2万帧（总计2千万秒），设置采样间隔为1秒，则能够获取到总共有2万帧的数据。基于不同的系统规模需求，在实际操作中建议设置采样数量在2千至5千帧之间。
结构划分：通过查看【处理pdb结构（进阶版）
3.5 分析结果
查看
查看
查看
该法通过计算MM方法得出气体相变自由能值（ΔGgas = ΔGvdw + ΔGele），其中利用GB和PB分别进行相应的自由能分量计算：对于气体相变自由能（ΔGsolv = ΔGpolar + ΔGnonpolar），则分别采用GB和PB方法进行求解。进一步分析可知，在GB/PB列中所对应的非极性和平极性相互作用自由能分量（即ΔGvdw 和 ΔGele）与总体气体相变自由能（ΔGgas）具有完全相同的数值特征；然而，在极性和非极性溶剂相互作用部分所对应的自由能分量（即ΔGpolar 和 ΔGnonpolar）则呈现出显著差异性。
综合自由能（Δ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 |
从结果可以看出，在此情况下复合物1 的结合能力最强。因此我们选择利用其分子动力学模拟的最后一帧的动态过程来对这一区域的结合模式进行详细分析。
当基于GB和PB的方法得出不同结论时, 可以通过分析各个能量项来判断, 从而筛选出那些具有有利能量特征的复合物体系。
3.5 结合模式分析
- 进入分子动力学模拟结果页面，下载分子结构。
这里下载**
7.prod-dry.pdb** 文件。
分子结构
该软件（如PyMOL）或个人使用的分子绘图工具能够启动对两个蛋白相互作用的分析。
分子动力学模拟可能会修改原子编号，请建议在完成分析后回到原始的原子编号以便后续分析；相关数据可以在上面准备复合物三维结构时保留的**
renumber.csv** 文件中找到。
从Excel表格文件mmpbsa.xlsx中的gb_total页面或pb_total页面可以看到，在配体中的GLU387与受体上的ARG151及ARG158之间存在较强的静电相互作用。将这些关键点在PyMOL软件中展示后发现，在图示中用黄虚线标出的位置形成了精确的离子键（即盐桥），其距离分别为3.4 Å和3.8 Å。同样地，在后续分析中我们还考虑了其他相关氨基酸残基。
