Hiv-B409

背景

体系参数准备

受体准备

# PDB 1EC1
# 75.1 /home/dddc/zzy/project/koff/hiv/hiv-B409/sys-pre/pro
# 安装pdb2pqr
pip install --upgrade pip
pip install pdb2pqr


conda activate pdb2pqr
# 加氢，去水
pdb2pqr30 1ec1.pdb 1ec1.pqr --ff AMBER --ffout AMBER --drop-water --with-ph 7.4 --pdb-output 1ec1-H.pdb
# 去除小分子
grep -v "HETATM" 1ec1-H.pdb >1ec1-pro.pdb
# 文献中测定koff的是 突变体（Q7K),故删去1ec1-pro.pdb中Gln7的侧链原子(只保留CA,C,O,N)，修改主链名称为Lys,tleap会自动补齐
# !!! 两条链都要操作！！！

# 4090 服务器上ambertools安装不下来，在85.24上运行（23上没找到 ff19SB)
# 85.24 /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/pro
scp dddc@172.21.75.1:/home/dddc/zzy/project/koff/hiv/hiv-B409/sys-pre/pro/* ./
tleap -f tleap-pdbclean.in
# 得到 1ec1-pdb2pqr-tleap.pdb
# 无二硫键，无离子


# 检查和师姐的序列是否一致
# /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/pro/check
cp /home/databank/lywu/vsremd_plus/protein-ligand/hiv-1-protease/a047/tleap/apo-H-align.pdb ./
cp ../1ec1-pdb2pqr-tleap.pdb ./
awk '{print $4}' 1ec1-pdb2pqr-tleap.pdb |uniq >zzy
awk '{print $4}' apo-H-align.pdb |uniq >wly
vi -d zzy wly
# 184个氨基酸

配体处理

# 85.24 /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/ligand
# 2018版薛定谔（默认力场 OPLS_2005),保留晶体结构的手性信息，PH与蛋白一致
/home/lywu/software/Schrodinger_Suites_2018/ligprep -WAIT -epik -ph 7.4 -g -isd B409.sdf -osd B409-ligpred.sdf


# 生成gjf文件
antechamber -i B409-ligpred.sdf -fi sdf -o lig.gjf -fo gcrt -at gaff2 -gn "%nproc=4" -gm "%mem=4GB" -gk "#B3LYP/6-31G* scrf=(solvent=water) pop=MK iop(6/33=2,6/42=6) opt" -rn MOL  -nc 0
# 经纠正，不能在水溶剂里优化，高斯关键词如下：

# 高斯优化
nohup g16 lig.gjf &



# /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/ligand/tleap-input
cp ../b3lyp-16/lig.log ./
antechamber -i lig.log -fi gout -o lig.prep -fo prepi -c resp
parmchk2 -i lig.prep -f prepi -o lig.frcmod -a y


# gv里分别打开 NEWPDB.PDB 与 B409.sdf（晶体结构配体文件）
# 打开GV tools-atom list-鼠标点击symbol列-rows-sort selected -拖动鼠标选择所有H-edit-delete-selected atoms），接着仍在atom list中，edit-z matrix-standardize，分别保存<步骤一>。
# 得到 g16-gv.pdb crystal-gv.pdb
# GV里 检查 g16-gv.pdb crystal-gv.pdb，更改crystal-gv.pdb （atom list - Tag)，使其Tag 与g16-gv.pdb 完全一致，另存为 crystal-gv2.pdb。


# /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/ligand/tleap-input/
dos2unix *
grep -v "CONECT" g16-gv.pdb |cut -b 1-26 > symbol
cat crystal-gv2.pdb | cut -b 29-78|sed '1,2d'|sed '1i Combineing the coordinates from the crystal structure with the atomic numbers from the Gaussian optimized structure'> coordinate
paste -d " " symbol coordinate|sed '1i TITLE      lig.pdb'> lig.pdb
rm symbol coordinate

复合物准备

# 85.24 /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/tleap
cp ../ligand/tleap-input/lig.pdb ./
cp ../pro/1ec1-pdb2pqr-tleap.pdb ./


# tleap
source leaprc.protein.ff19SB
source leaprc.gaff2
source leaprc.water.opc
loadamberparams /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/ligand/tleap-input/antechamber/lig.frcmod
loadamberprep /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/ligand/tleap-input/antechamber/lig.prep
mol=loadpdb lig.pdb 
pro=loadpdb 1ec1-pdb2pqr-tleap.pdb
com=combine {pro mol}
saveamberparm com native.prmtop native.inpcrd
savepdb com com-dry.pdb
charge com
solvatebox com OPCBOX 12.0
addionsrand com Cl- 6
charge com
saveamberparm com complex.prmtop complex.inpcrd
savepdb com com.pdb

# amber参数转为gmx参数，并略作修改为gmx认的得的标准化版本，添加位置限制信息
mkdir gmx
cd gmx 


python 
import parmed as pmd
amber = pmd.load_file('../complex.prmtop','../complex.inpcrd')
amber.save('gmx.top')
amber.save('gmx.gro')


# 添加位置限制信息用于使得在预平衡时溶质位置基本不变
# ！！！我的体系，gmx.top中两条链各自为 protein1,protein2
# 0908 已解决，出现该问题主要是因为两条链不完全一致
# 在蛋白质结束和小分子开始前添加位置限制（moleculetype）
sed -i '17597 a ; Include Position restraint file\n#ifdef  POSRES\n#include "posre1.itp\n#endif\n' gmx.top
# 在小分子结束和离子开始前添加位置限制（moleculetype）
sed -i '18603 a ; Include Position restraint file\n#ifdef  POSRES\n#include "posre2.itp\n#endif\n' gmx.top


# 位置限制文件只是个索引，规定了哪些原子被限制
gmx_mpi genrestr -f gmx.gro -o posre1.itp 
# 看起来只能写一条链的index (1-1596)
# Group    13 (MOL) has   104 elements 
sed -n '1,108p' posre1.itp > posre2.itp

体系平衡

# 85.24 /home/databank/zzy/project/MD/koff/hiv/hiv-B409/pre-equilibrium


mkdir em nvt npt mdp


cp /home/databank/lywu/vsremd_plus/protein-ligand/hiv-1-protease/a047/pre-equilibrium/mdp-file/* ./mdp/
cp ../sys-pre/tleap/gmx/* ./


### 能量最小化
gmx_mpi grompp -f ../mdp/em.mdp -c ../gmx.gro -r ../gmx.gro -p ../gmx.top -o em.tpr 
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm em &


### 等温等容预平衡
# 修改为1fs
gmx_mpi grompp -f ../mdp/nvt.mdp -c ../em/em.gro -r ../em/em.gro -p ../gmx.top -o nvt.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm nvt &
# 修改为2fs
gmx_mpi grompp -f ../mdp/nvt.mdp -c nvt.gro -r nvt.gro -p ../gmx.top -o nvt2.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm nvt2 &


### 等温等压预平衡
# 修改为1fs
gmx_mpi grompp -f ../mdp/npt.mdp -c ../nvt/nvt2.gro -r ../nvt/nvt2.gro -p ../gmx.top -o npt.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm npt &
# 修改为2fs
gmx_mpi grompp -f ../mdp/npt.mdp -c npt.gro -r npt.gro -p ../gmx.top -o npt2.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm npt2 &


mkdir smd-input
cd smd-input
cp ../gmx.top ./
cp ../npt/npt2.gro ./npt.gro

SMD

# SMD要到85.23上跑
# 85.23 /home/databank_70t/zzy/project/koff/hiv/hiv-B409/smd 
scp lywu@172.21.85.24:/home/databank/zzy/project/MD/koff/hiv/hiv-B409/pre-equilibrium/smd-input/* ./


# 初始文件：gmx.top  npt.gro md.mdp  md.mdp  smd.sh  splumed.in（dat文件模板）
cp /home/databank_70t/zzy/project/koff/trypsin-31/smd/splumed.in ./
cp /home/databank_70t/zzy/project/koff/trypsin-31/smd/md.mdp ./
cp /home/databank_70t/zzy/project/koff/trypsin-31/smd/smd.sh ./

# 在npt.gro里找到对应的原子index,修改 splumed.in 文件里的index
# ligand: COM ATOMS=421,1990
 25ASP     CG  421   3.211   4.129   3.994 -0.8174  0.4479  0.2517
 25ASP     CG 1990   3.237   3.904   3.543  1.1701 -0.2484  0.1551
# active-residue: COM ATOMS=3139-3242
  100MOL    C34 3139   3.335   3.146   4.472  0.8458 -0.0524  0.3882
  100MOL    H49 3242   2.181   3.630   3.872 -0.3864  1.3420  1.1905

# 以splumed.in为模板生成等同于副本数的dat文件,更改其中CV的值（即将配体拉出去的距离）
for i in {0..39};do m=`echo $i|awk '{print 0.20+$i*0.10}'`; sed "s/temp1/${m}/g" splumed.in > smd.${i}.dat;sed -i "s/COLVAR/COLVAR.${i}/g" smd.${i}.dat; done


# 检查smd.sh的核数，mdp的步数
nohup sh smd.sh &


# pymol 里看SMD的结果，挑选20个构象
# 85.23 /home/databank_70t/zzy/project/koff/hiv/hiv-B409/pre-equilibrium/smd/gro-check


# 85.23 /home/databank_70t/zzy/project/koff/hiv/hiv-B409/smd/gro-check/18replica-select
# 将挑选好的构象放入 vsremd-select
# npt 也算一个（结合态）
cp ../npt.gro ./
for i in {0,6,13,14,15,17,18,19,20,22,23,24,26,27,29,30,39};do cp ../smd${i}.gro ./;done
# 重命名为0-20，便于vsREMD操作
ls > index
for i in {0..17};do a=`expr ${i} + 1`; b=`sed -n "${a}p" index`;mv $b smd-select${i}.gro;done

vsREMD

# 到24 上跑vsREMD 100个核
# 85.24 /home/databank/zzy/project/MD/koff/hiv/hiv-B409/vsremd


# 将SMD结果中挑选出的20个构象从23上传输过来
scp -r dddc@172.21.85.23:/home/databank_70t/zzy/project/koff/hiv/hiv-B409/smd/gro-check/18replica-select ./smd-select 
# 从师姐那儿复制来
scp dddc@172.21.85.23:/home/databank_70t/lywu/hiv-1/a047/vsremd/test-100ns/temp ./
scp dddc@172.21.85.23:/home/databank_70t/lywu/hiv-1/a047/vsremd/test-100ns/md.mdp ./


# 修改mdp，跑1ns的测试
# 制作index：Protein_MOL Cl-_Water
gmx_mpi make_ndx -f ../smd-select/smd-select0.gro -o index.ndx
> 1|13
> 14|15


# 创建并进入1ns测试路径 85.24 /home/databank/zzy/project/MD/koff/hiv/hiv-B409/vsremd/1ns-test
source ~/.gmx2022.5-vsremd.sh
nohup sh -c 'for i in {0..17};do a=`expr ${i} + 1`; b=`sed -n "${a}p" temp` ; mkdir md$i; sed "s/300.00/${b}/g" md.mdp > md$i/md.mdp; cd md$i; gmx_mpi grompp -f md.mdp -p ../gmx.top -n ../index.ndx -r ../../smd-select/smd-select${i}.gro -c ../../smd-select/smd-select${i}.gro -o md.tpr; cd ../; done'


# 85.24 /home/databank/zzy/project/MD/koff/hiv/hiv-B409/vsremd/1ns-test
(base) [lywu@localhost 1ns-test]$ nohup mpirun -np 108 gmx_mpi mdrun -v -deffnm md -multidir md0 md1 md2 md3 md4 md5 md6 md7 md8 md9 md10 md11 md12 md13 md14 md15 md16 md17 -replex 1000 &
[1] 236353