【殷赋云教程】蛋白与多肽的对接计算教程
本教程讲述如何使用CABS-dock网站、殷赋云平台分子动力学方案及相关小工具，实现蛋白与多肽的对接计算。
本次示例中，受体蛋白分别为3BP1（PDB ID : 1ABO），配体多肽序列为APTMPPPLPP。该体系并不
清楚具体的氨基酸残基位点，所以采用盲对接方式确定结合位点。首先采用CABS-dock进行刚性对接，
挑选若干复合物进行分子动力学模拟，然后根据结合自由能选取最好的结果进行结合模式分析。
相对于明确位点的对接方式来说，盲对接的花费较大。因此，建议尽量通过文献报道、实验数据确定作用残基后再进行对接，选择最好的复合物进行下一步计算；
以下网站（软件）也可进行蛋白多肽对接：MDockPeP、HPEPDOCK、pepATTRACT、
FlexPepDock、GalaxyPepDock、HADDOCK peptide docking、DINC和AutoDock
CrankPep (ADCP)
打开殷赋云计算平台（https://cloud.yinfotek.com/)【处理pdb结构】小工具。
输入PDB ID 1ABO（或上传PDB文件），删除多余结构，保留用于对接的蛋白肽链。
处理原则：删除不参与对接的肽链、杂原子及水分子。
处理后，将receptor.pdb改名为1ABO.pdb。
打开CABS-dock在线服务网页（http://biocomp.chem.uw.edu.pl/CABSdock/）进行刚性对接。Optional选项，请根据实际情况填写。
更多CABS-dock网站的说明，请看Help说明，我们不再解答该网站的相关内容。
一段时间后，网站就会得到多个复合物结果。这里，我们选择其中5个进行下一步计算。
这里以complex1.pdb为例，介绍分子动力学模拟及结合自由能计算的过程，其他复合物同理。
打开殷赋云平台【处理PDB结构（进阶版）】小工具。
上传复合物结构，点击【下一步】。
点击【生成文件】，下载文件。
有关PDB修复的内容，详见教程《处理PDB结构（进阶版）》
默认勾选对残基重新编号，renumber.csv记录了残基编号的对应关系，方便后续分析，切记保存好。
这里虽然显示receptor.pdb，但实际上是复合物结构。
打开殷赋云平台【准备Amber文件】小工具。
上传前面准备好的receptor.pdb文件，点击【生成】，下载文件。
无需上传参数文件；
amber_files.zip中的residue_numbers.csv文件包含了蛋白、离子和水的残基编号，方便分子动力学模拟时设置约束条件。
打开殷赋云平台【分子动力学（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 这两个文件，参考步骤7Production的时长分段建议设置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采用以下约束条件：
待任务完成，点击【查看】，进入分析页面。查看热力学性质。
在平衡（Equilibration）或生产（Production）阶段，能量、温度、压强、体积等曲线应当平稳 。详情见《分子动力学模拟（Amber 20）》
查看生产（production）阶段的RMSD。该曲线在某个高度附近平稳振动，振幅大小在2Å以内，表明体系已达到平衡，可以进一步分析。
分子结构。可以下载各个阶段的结构文件，使用PyMOL等软件查看结构。
打开殷赋云平台【结合自由能计算】方案，创建任务，进入任务页面。
注意：体系平衡后才能进一步计算结合自由能，否则，误差可能会很大。
从云盘文件选择拓扑文件及轨迹文件，并设置参数，然后【提交】任务。
采样范围：总共20000帧（20 ns），采样间隔为100，则采样数为20000 ÷ 100 = 200帧。可根据体系大小进行设置，一般采样200-500帧即可。
结构划分：查看【处理pdb结构（进阶版）】生成的renumber.csv文件，了解受体和配体蛋白的残基编号范围。如下图，受体为A链，对应残基编号1-299，配体为B链，对应残基编号300-551。
待任务完成，点击【查看】，进入分析页面，查看结合自由能，并下载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}} $不利的复合物体系。
进入分子动力学模拟结果页面，下载分子结构。
这里下载7.prod-dry.pdb文件。
用PyMOL或自己熟悉的分子图形软件打开分析复合物之间的相互作用。
分子动力学模拟会更改原子编号，请先改回原来编号再进行分析，可查看上面准备复合物三维结构时保留的renumber.csv文件。
从mmpbsa.xlsx文件的gb_total或pb_total页可以看出，配体的LEU与PRO65与受体的TYR52和ASN51都有较强的静电力作用。将它们显示在PyMOL上，可见它们之间形成氢键（图中绿色虚线表示），氢键距离分别为2.7Å和2.9 Å。同理，再分析其他氨基酸残基。
