cMD

1.蛋白前处理

1.1补全侧链

PDB网站获得的结构往往发生残基缺失，需要进行补全。在使用同源建模。前者的好处是使用的晶体结构，更准确。后者虽然操作简便，但尽管同源建模和晶体结构RMSD很小，实际上有很多残基都会发生一定角度的偏离，有时会造成结果的定性错误。因此需要根据实际情况，选择使用的方法，有些问题可以使用同源建模（比如根本没解析晶体结构），有些问题建议使用晶体结构。

1.1.1同源建模

以当前PDB结构作为模板，生成残基完整的结构。
在PDB下载结构，FASTA序列。
删去蛋白结构中所有配体，cofactor，水（如果有重要水可将其另存为新文件，后面会用），离子，多聚体可只保留一个。
进入SWISS-MODEL Interactive Workspace (expasy.org)
右侧选user template（如果没有晶体结构则不选这个），上传PDB文件，FASTA序列。网站自动以上传的结构为模板，对序列建模。
完成后下载，命名为3RX4_swissmodel.pdb。
使用该方法可能会导致建模的坐标和原始坐标不一样，可以进行下一步。

1.1.2补全晶体

使用pymol手动补全晶体结构缺失的残基和原子。pymol加载建模结构及晶体结构，建模结构align to 晶体结构，保存。文本打开晶体结构，检索missing，找到缺失信息，在建模结构中找到对应的，复制到晶体文件对应部分。

1.2质子化

从PDB网站提供的文献中获得pH，预测此时蛋白的质子化情况。

方法1

cpu256上需要先激活（有的直接pdb2pqr直接就能调用）：
  • conda activate pdb2pqr
使用pdb2pqr：
  • pdb2pqr30 input.pdb out.pqr --ff AMBER --ffout AMBER --with-ph 7.4 --pdb-output out.pdb
input和out改为自己的输入输出文件，设置输入输出力场是amber，—chain为保留原始的原子编号（不从1开始，如果不设置，默认从1开始），--with-ph为结晶状态下的pH值，并自主处理末端氨基酸，即在C末端增加OXT（对于中间缺失的loop，不会增加OXT）。输出为*.pqr。注意代码中的横线是两个减号。

方法2

使用H++网站预测（没用过，amber官网用的这个）
H++ (web-based computational prediction of protonation states and pK of ionizable groups in macromolecules) (vt.edu)
登录按提示操作即可
最后，质子化完成后检查文件中的HIS是否被自动改为HID或者HIE

2.小分子前处理

如果是从PDB下载的晶体结构，可以再下载下面的小分子结构，坐标和复合物中的是一致的，后续步骤也能正确显示。其他情况可以使用GV打开原始PDB文件，tools-atom list选择小分子，rows-invert select反选，edits-delete atoms删去其他的原子，保存为pdb文件。

2.1加氢

如果是从对接得到的结构，小分子已经有氢，则不需要此步骤。
文件读入maestro，然后Tasks – Ligand Preparation – 修改其中Generate possible states at target pH（改为结晶实验的pH，如7.0+-0），run。结束后左侧任务名下任选一个结构。重新用GV打开配体文件，询问是否加氢时选是，对照着maestro的结果在GV里为小分子加氢或删去氢，保存为gjf文件。
注意：若有重要的水分子，提取出来，加氢之后，保存，再和蛋白质及小分子共同存成一个pdb。最后用GW打开，生成优化文件。

2.2参数计算

因为gromacs没有运行小分子所需的静电势（RESP），力场信息等，需要先用量化计算波函数，进行参数拟合，得到这些信息。
使用记事本打开gjf文件，将开头%，#，后面的内容替换为下面的。注意下面第8行，第一个数字表示电荷，第二个数字表示自旋多重度。电荷改为实际电荷，自旋多重度在没有过渡金属的情况下写1。在文件最后添加ante的两段。确保整个文件最后有至少两个空行。
  • %mem=4GB
    • 

  • %nproc=4
    • 

  • # opt b3lyp/6-311g(d,p) geom=connectivity pop=mk iop(6/33=2,6/42=6,6/50=1) 
    • 


    • 

  • 3RX4_ligand
    • 


    • 

  • -1 1
    • 

  •  F                 -5.74000000   10.85800000   15.05100000
    • 

  •  S                -14.53700000   12.33400000   19.06700000
    • 

  •  C                 -6.89000000    7.26900000   18.19600000 
    • 

  •  
    • 

  •  1 15 1.0
    • 

  •  2 23 1.0 24 1.0
    • 

  •  3 5 1.0 7 2.0 10 1.0
    • 


    • 

  •  antechamber-ini.esp
    • 


    • 

  •  antechamber.esp
    • 

  •  
    • 

  •  
opt：几何优化。以前有教程是在真空中使用HF/6-31G(d)，是因为此级别会高估偶极矩10%-15%，正好等效反映出水的极化作用。但这种高估效果并不稳定。因此使用b3lyp/6-311g(d,p)在隐式溶剂模型scrf=(solvent=water)中计算显然更好，不要使用HF。geom确保原子之间键连正确。
iop(6/33=2)是进行RESP Fitting并输出到Gaussian的.log文件。
IOp(6/42)：Density of points per unit area in esp fit. 0：Default (1) N：Points per unit area。
IOp(6/50)：Whether to write Antechamber file during ESP charge fitting.
pop=mk：Produce charges fit to the electrostatic potential at points selected according to the Merz-Singh-Kollman scheme 
提交Gaussian计算。
  • nohup g16 3RX4_ligand.gjf &
如果是g09就写g09
完成后，antechamber提取小分子静电势信息，得到prep和NEWPDB文件。
  • antechamber -i 3RX4_ligand.log -fi gout -o 3RX4_ligand.prep -fo prepi -c resp
-i   input file name
-fi  input file format
-o   output file name
-fo  output file format
-c charge method
parmchk2提取小分子力场参数信息，得到frcmod文件
  • parmchk2 -i 3RX4_ligand.prep -f prepi -o 3RX4_ligand.frcmod -a y
-a print out all force field parameters including those in the parmfile

2.3位置调整

如果在1.1使用了方法1，得到的坐标和原始PDB坐标会有细微差别，如果可以接受，跳过此步骤。
2.3通过对接将小分子重新放到口袋中，或者提供多个起始构象。

2.3.1明确口袋

知道小分子所在口袋的情况使用该方法。
进入pymol，加载带小分子的晶体结构，选中小分子，为sele。
  • centerofmass sele
记录小分子质心，作为后续对接的盒子中心。
  • obabel NEWPDB.PDB -osdf -O ligand.sdf
将NEWPDB.PDB转成sdf文件进行对接。上传文件至服务器。
  • smina --seed 1 --num_modes 10 -r complex.pdb -l 0415_g16.sdf --center_x  81.581 --center_y 29.646 --center_z 12.226 --size_x 10 --size_y 10 --size_z 10 -o docking.pdb
得到对接结果，选择其中之一或全部。

2.3.2口袋未知

不知道分子在哪个口袋，需要使用盲对接。
下载pymol脚本，用来划分格子。
drawgridbox.py
将该脚本放在pymol的工作路径中。
  • load drawgridbox
    • 

  • drawgridbox repc0, nx=5, ny=5, nz=5,  lw=1, g=0, b=0
    • 

  • centerofmass
可以得到Box dimensions和Center of Mass，记录。

2.4原子名调整

之前计算参数时，由于Gaussian的原子名与直接抠出的小分子的原子名不一样，在后续分别读取两个文件的参数和小分子坐标时，会因为原子名不一样而无法匹配，因此需要将原子名修改一致。但是经Gaussian和antachamber处理后的原子顺序和抠出的pdb文件的不一样。例如，
pdb文件：
HETATM    1  F   SFI A 317      -5.740  10.858  15.051                       F
HETATM    2  S   SFI A 317     -14.537  12.334  19.067                       S
处理后的NEWPDB.PDB文件：
HETATM   15  F1  MOL     1      -3.136   4.050   0.212                       F
HETATM   23  S1  MOL     1       6.251   0.712  -0.140                       S
可以看到原子序号（第二列），原子名（第三列）等都不一样。
老办法是打开可视化软件，展示原子序号，眼睛观察，找到两个文件中的同一个原子，例如上面的两个F原子，然后将NEWPDB的‘HETATM   15  F1  MOL   1’复制取代pdb的‘HETATM    1  F   SFI A 317’。对所有原子进行这一操作。
现找到了新的方法能更方便的处理。删除NEWPDB.PDB中的氢原子（各种方法都可以，如打开GV tools-atom list-鼠标点击symbol列-rows-sort selected -拖动鼠标选择所有H-edit-delete-selected atoms），接着仍在atom list中，edit-z matrix-standardize，保存。对抠出小分子的pdb文件同样标准化。文本打开两文件，可以发现每个原子的index是对应的，文本按列操作，复制NEWPDB所有原子的‘F1  MOL     1’三列到pdb文件对应位置覆盖即可。
新方法的原理在于GV有一套对原子编号的算法，确保每个化合物有唯一原子编号，因此只要两个分子结构相同，就可以获得相同的编号，也就使得了两个文件中所有原子按顺序一一对应了。

3.准备运行文件

3.1方法1：Amber转

3.1.1Amber生成文件

由于使用gmx生成参数文件流程繁琐，一种方法是先生成Amber参数，再转为gmx参数。
打开amber自带的leap程序，准备运行文件
  • tleap
加载蛋白力场，小分子力场。oldff/leaprc.ff14SB也可以是leaprc.protein.ff14SB。
  • source oldff/leaprc.ff14SB
    • 

  • source leaprc.gaff
oldff/leaprc.ff14SB中带了核酸，水等一系列参数，就不需要额外加载水了。如果使用后者，需要额外加载
  • loadamberparams frcmod.ionsjc_tip3p
    • 

  • loadamberparams frcmod.ions234lm_1264_tip3p
  • loadamberparams  3RX4_ligand.frcmod 
    • 

  • loadamberprep 3RX4_ligand.prep   
  • mol=loadpdb 3RX4_ligand.pdb    
    • 

  • pro=loadpdb 3RX4_complex.pdb
  • bond ramp.27.SG ramp.82.SG
有二硫键需要进行连接。检索原PDB文件中的SSBOND可以知道。
（小分子，蛋白质分开存储，输出可视化pdb）
  • com=combine {pro mol}
    • 

  • saveamberparm pro pro.prmtop pro.inpcrd
    • 

  • saveamberparm mol lig.prmtop lig.inpcrd
    • 

  • saveamberparm com com.prmtop com.inpcrd
    • 

  • savepdb com com-look.pdb
  • charge com
体系距离盒子边界至少10埃
  • solvatebox com TIP3PBOX 10
Total unperturbed charge:  -3.000001
Total perturbed charge:    -3.000001
> solvatebox com TIP3PBOX 10
  Solute vdw bounding box:              57.758 122.585 111.395
  Total bounding box for atom centers:  77.758 142.585 131.395
  Solvent unit box:                     18.774 18.774 18.774
  Total vdw box size:                   81.057 145.508 134.008 angstroms.
  Volume: 1580540.610 A^3 
  Total mass 823480.202 amu,  Density 0.865 g/cc
  Added 43687 residues.
  • addionsrand com Na+ 83 Cl- 74
（检查，有浮点数无妨）
  • charge com
可视化检查
  • saveamberparm com complex.prmtop complex.inpcrd 
    • 

  • savepdb com com.pdb
退出
quit

3.1.2 转gromacs参数

  • python
  • import parmed as pmd
    • 

  • amber = pmd.load_file('complex.prmtop','complex.inpcrd')
    • 

  • amber.save('topol.top')
    • 

  • amber.save('gromacs.gro')
WAT，system1是gmx不能识别的标识符，需替换为SOL，Protein，修改topol.top文件。在vim里先输入冒号，再粘贴后面的内容（未知原因每次输入到第二个/都会退出，但复制粘贴可以）。
  • vi topol.top
    • 

  • :1,$s/WAT/SOL/g
    • 

  • :1,$s/system1/Protein/g
后续有步骤需要对体系进行限制动力学，这里在top文件中添加限制代码。文本打开top文件，检索moleculetype，第二个应该是小分子（MOL），在蛋白质坐标结束和小分子坐标开始前添加以下内容:
同理，在小分子坐标结束和离子开始前添加以下内容：
后面每个步骤都需要使用该top文件，但不用担心每一步都会限制，在后面使用的mdp文件中，参数define控制是否进行限制。define = -DPOSRES  ; position restrain the protein
修改gromac.gro文件
  • :1,$s/WAT/SOL/g
生成索引文件
  • gmx_mpi make_ndx -f gromacs.gro -o index.ndx
输入 1 |  17 （中间要有空格）
生成posre.itp文件（规定用于限制体系运动的谐正势参数）
  • gmx_mpi genrestr -f gromacs.gro -n index.ndx -o posre.itp
限制蛋白，输入 1
  • gmx_mpi genrestr -f gromacs.gro -n index.ndx -o posre2.itp 
限制小分子，输入 17
posre2.itp要从1开始编号，比如可以把posre.itp的序号按列复制到这里。

3.2方法2：Gmx直接生成

4.  MD

后续使用的所有sh，mdp文件：
动力学文件.7z

4.1  更改环境变量

  • source lywu.sh

4.2  能量最小化 

目的是消除体系中可能存在的过大斥力，避免模拟一开始，某些原子速度过快出现不稳定或崩溃。
  • gmx_mpi grompp -f minim.mdp -c gromacs.gro -p topol.top -o em.tpr
    • 

  • nohup mpirun -np 20 gmx_mpi mdrun -v -deffnm em & 
np：使用的核数。mdrun是gmx程序，负责管理进程。

4.3  平衡 

对蛋白原子做限制性动力学，使得水弛豫开之前蛋白不会瞎动，避免不测
  • gmx_mpi grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr 
    • 

  • nohup mpirun -np 20 gmx_mpi mdrun -v -deffnm nvt &

4.4  升温

  • gmx_mpi grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
    • 

  • nohup mpirun -np 20 gmx_mpi mdrun -v -deffnm npt &

4.5  常规MD测试

  • gmx_mpi grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index -o md_0_1.tpr
    • 

  • nohup mpirun -np 150 gmx_mpi mdrun -v -deffnm md_0_1 &
一键准备体系
傻瓜式操作网站https://charmm-gui.org/
构建膜蛋白体系可以使用。不足之处是如果报错，不知道报错原因，难以处理，因此简单的体系可以不用。当然用也是可以的。