分子动力学模拟计算新冠病毒S蛋白和抗体结合自由能

参照教程,为了适合S蛋白和抗体，对原教程进行了大量适用性修改，得到结果不保证准确性（甚至是错误的），但好在流程走通，得到了一个数据。伞形采样(Umbrella Sampling)http://www.mdtutorials.com/gmx/umbrella/index.html

准备复现的是上面的论文，所以也也部分参考里面的参数。

云计算平台：北鲲云

使用的软件是云计算平台的Gromacs，下面两个是用到的版本号。

module add GROMACS/2019.6-intel-2019b

GROMACS/2020-foss-2019b

一、准备蛋白质

使用的是7K8M（Structure of the SARS-CoV-2 receptor binding domain in complex with the human neutralizing antibody Fab fragment, C102）

在这个网站下载：3D View: 6VXX (rcsb.org) 搜索然后点右上角下载就行。

蛋白质长这样，左边是C102抗体，右边是S蛋白，可以看到有一个配体。

用记事本打开pdb文件，可以看到由三条链组成：

A、B抗体蛋白，E链是S蛋白。

在记事本中删除配体部分，HETATM开头的行删去。PDB文件解读看这篇文章。PDB(Protein Data Bank)数据格式详解

删除的配体部分用Pymol打开是看不到那部分的（这里假设这个配体不影响结合自由能计算。）

后面部分在Gromacs中进行。

先把处理后的PDB上传到云平台，新建文件夹用于计算，SSH连接。

Run the structure through pdb2gmx:

gmx pdb2gmx -f 7K8M.pdb -ignh -ter -o complex.gro

选 GROMOS96 53A6

SPC water,

"None" for the N-termini,

"COO-" for the C-termini

在产生的文件topol_Protein_chain_B.itp文件最后加上下面这段话

#ifdef POSRES_B
#include "posre_Protein_chain_B.itp"
#endif在topol_Protein_chain_A.itp后面也加上

#ifdef POSRES_A

#include "posre_Protein_chain_A.itp"

#endif

可以固定A、B链。（总共会产生A\B\E六个itp文件）

第二步定义盒子

因为蛋白质大小差距较大，这里必须要改，盒子大小按论文尺寸。位置按盒子宽度的一半，x方向可以调整生成构型后再查看反复调整

gmx editconf -f complex.gro -o newbox.gro -center 15 3.275 3.625 -box 20 7.55 7.25 -princ

加-princ可以把蛋白质摆正

可以用Pymol查看盒子，显示CELL就能看到盒子。

第三步添加离子

添加水

gmx solvate -cp newbox.gro -cs spc216.gro -o solv.gro -p topol.top

添加离子按原教程mdp文件

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.1

这里蛋白会显示在盒子的尾部，但最后牵引的时候还是能拉开（好像没有地方设置了牵引力的方向）

现在的样子。

第四步能量最小化

这里的mdp文件也不需要修改

gmx grompp -f minim.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm emgmx grompp -f npt.mdp -c em.gro -p topol.top -r em.gro -o npt.tpr
gmx mdrun -deffnm npt

第五步生成采样构型

这一步在力的牵引下形成多个构型

pull_coord1_dim = Y N N
pull_coord1_groups = 1 2
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k = 1500 ; kJ mol^-1 nm^-2
pull-group1-pbcatom = 370

加粗的两行要在原来的mdp文件改一下，第一个是牵引方向，我们在x轴方向牵引

最后一个是对齐原子，370是抗体蛋白中心的原子的序号

后续的所有mdp都要这么改！！！

分组

gmx make_ndx -f npt.gro
...
> r 1-214
> name 19 Chain_A
> r 334-517
> name 20 Chain_B
> q

这两个数查看pdb文件，前214个氨基酸是抗体的，后面的是S蛋白的

然后运行MD模拟

gmx grompp -f md_pull.mdp -c npt.gro -p topol.top -r npt.gro -n index.ndx -t npt.cpt -o pull.tpr
gmx mdrun -deffnm pull -pf pullf.xvg -px pullx.xvg

分片，先把下面要用的两个文件移到新的文件夹再运行，因为会产生501个文件，按照原教程运行脚本生成distance文件，生成的文件前面一半数据可能没有第二列，数据在后面。
gmx trjconv -s pull.tpr -f pull.xtc -o conf.gro -sep

把数据复制到excel里，然后一个一个找需要用的文件的序号，记下来。

大概就这样。间隔0.2nm。

第六步伞形采样

终于到运行模拟了！注意改mdp

gmx grompp -f npt_umbrella.mdp -c conf6.gro -p topol.top -r conf6.gro -n index.ndx -o npt0.tpr

上面这个运行20遍，没次改数字

或者靠脚本完成，简单脚本，里面数字自己改。

#!/bin/bash
#################################################
j=0
for i in 0 62 96 157 186 230 268 356 433 491
dogmx grompp -f npt_umbrella.mdp -c conf${i}.gro -p topol.top -r conf${i}.gro -n index.ndx -o npt${j}.tprlet "j++"
done
exit;

然后mdrun，可以把所有npt文件复制到一个文件夹再继续。
gmx mdrun -deffnm npt0

运行20遍

或者同样脚本完成：（这个脚本只针对云计算）

#!/bin/bash
for i in `ls /home/cloudam/work/tpr/*.tpr`;do
job=${i%.*}
cat >> $job.sh << EOF
#!/bin/bash
module add GROMACS/2019.6-intel-2019b
mpiexec -v  gmx_mpi mdrun -v  -cpi $job -deffnm $job
EOF
cd /home/cloudam/jobs/mn/
sbatch -p c-16-1 -n 16 -c 1 $job.sh
done

很久很久以后。

运行md模拟：注意改mdp

gmx grompp -f md_umbrella.mdp -c npt0.gro -t npt0.cpt -p topol.top -r npt0.gro -n index.ndx -o umbrella0.tpr

*20

gmx mdrun -deffnm umbrella0

*20

同样可以脚本完成，改一下上面的就OK

第七步数据分析

在结果文件夹运行

ls umbrella*_pullf.xvg > pullf.dat

ls umbrella*_pullx.xvg > pullx.dat

ls *.tpr > tpr.dat

生成需要的三个dat文件

gmx wham -it tpr.dat -if pullf.dat -o -hist -unit kCal

会生成profile.xvg和hist.xvg,把里面的数据取出来作图。

取样太少，结果不理想。好像也没有收敛，数据和论文中的差了一半。

Done！

问题：

这样分几个文件夹进行会经常缺文件，这时候运行cp命令把需要的文件复制进去就行。

力场和论文中的不一样，中途报了几个warning，直接忽略了，肯定影响结果。

水模型和论文中不一样。

取样太少（为了省钱）

最后的MD模拟可能有问题，两个小时没有运行完被shot down了。对比了运行一个小时和两个小时的数据，一个小时的数据显然更正常。这应该是个大问题

算失败了一半，但还好。