Tutorial on running Molecular Dynamics for GIST grid generation

From DISI
Revision as of 02:27, 22 January 2017 by TBalius (talk | contribs) (→‎Run AMBER)
Jump to navigation Jump to search

Tutorial written by Trent Balius (Jan. 9, 2017).

Disclaimer

This is foremost for training Shoichet lab members. But we hope that the community finds this useful.

To use this tutorial you will need the following:

  • AMBER12 or higher.
  • AMBERTOOLS.
  • GPUs to run AMBER on. In our case we uses an SGE queuing system to manage these jobs.
  • It is helpful if you have DOCK3.7 the latest version.

Introduction

Grid inhomogeneous Solvation Theory (GIST) is a method for calculating thermodynamic properties of water on GIST lattice.

Here, we will use these grids for DOCKing (although, there are many other uses).

Set up environment

you can add the following to your environment file (.cshrc) for just issue them on the command line.

source python with bio_python package.

source /nfs/soft/python/envs/complete/latest/env.csh

set DOCKBASE

setenv DOCKBASE "/nfs/home/tbalius/zzz.github/DOCK"

set AMBERHOME

 setenv AMBERHOME /nfs/soft/amber/amber14

Put msms in your path

set path = ( /nfs/home/tbalius/zzz.programs/msms $path )

Set up directories

cd /nfs/work/yourusername/
mkdir ambergisttutorial
cd ambergisttutorial 

Prepare for AMBER

  • download 4NVE form the PDB brake it into receptor and ligand pieces.
  • download pdb 4NVA form the PDB align it to 4NVE and out put waters that are close to the ligand.
  • here is one way to do it:

(1) Break up the protein from the ligand. run be_blasti.py

 $DOCKBASE/proteins/pdb_breaker/be_blasti.py --pdbcode 4NVE nocarbohydrate original_numbers | tee -a pdbinfo_using_biopython.log
 $DOCKBASE/proteins/pdb_breaker/be_blasti.py --pdbcode 4NVA nocarbohydrate original_numbers | tee -a pdbinfo_using_biopython.log
 cd 4NVE 
 sed -e "s/HETATM/ATOM  /g" lig.pdb > xtal-lig.pdb
 cd ../
 cd 4NVA
 grep HOH $mountdir/workingdir/4NVA/4NVA_A.pdb > $workdir/water.pdb
 cd ../

(2) Next align the two receptors (4NVE to 4NVA) and find those waters close to the ligand.

 mkdir aligned
 cd aligned

create file named "chimera.com" and write the following into the file:

 # Gist template #0
 open ../4NVA/rec.pdb
 # 4NVE rec
 open ../4NVE/rec.pdb
 # 4NVE lig
 open ../4NVE/lig.pdb
 # 4NVA waters
 open ../4NVA/water.pdb
 mmaker #0 #1
 matrixcopy #1 #2
 sel #3:HOH & #2 z<8
 del #3 & ~sel
 write format pdb  0 template.pdb
 write format pdb  1 aligned.rec.pdb
 write format pdb  2 aligned.lig.pdb
 write format pdb  3 close.waters.pdb
 /nfs/soft/chimera/current/bin/chimera --nogui chimera.com > & chimera.com.out
 cd ../

(3) Now lets run tleap.

(3.1) run leap to renumber residues.

mkdir tleap
cd tleap
cp ../aligned/template.pdb rec.pdb
cp ../aligned/close.waters.pdb wat.pdb

produces tleap input file named tleap.rec.1.in

set default PBradii mbondi2
# load the protein force field
source leaprc.ff14SB
loadamberparams frcmod_teb_mf_mod.hemall
loadamberprep heme_all_plus1.in
REC = loadpdb rec.pdb
# draw bond between HIS and HEME
# 5 coordenated iron.
bond REC.175.NE2 REC.293.FE
saveamberparm REC rec.1.leap.prm7 rec.1.leap.rst7
quit
$AMBERHOME/bin/tleap -s -f tleap.rec.1.in > ! tleap.rec.1.out
$AMBERHOME/bin/ambpdb -p rec.1.leap.prm7 < rec.1.leap.rst7 >! rec.1.leap.pdb

3.2 requires script 0000 to be run first and align/ dir to be present

reduce -HIS -FLIPs rec.1.leap.pdb >! rec.nowat.H_mod.pdb
sed -i 's/HETATM/ATOM  /g' rec.nowat.H_mod.pdb
grep "^ATOM  " rec.nowat.H_mod.pdb | sed -e 's/   new//g' | sed 's/   flip//g' | grep -v "OXT" | grep -v " 0......HEM" >! rec.nowat.H_mod1.pdb
curl docking.org/~tbalius/code/waterpaper2017/scripts/replace_his_with_hie_hid_hip.py > replace_his_with_hie_hid_hip.py
curl docking.org/~tbalius/code/waterpaper2017/scripts/replace_cys_to_cyx.py > replace_cys_to_cyx.py
curl docking.org/~tbalius/code/waterpaper2017/scripts/add.ters.py > add.ters.py
python replace_his_with_hie_hid_hip.py rec.nowat.H_mod1.pdb rec.nowat.H_mod2.pdb

The following script is useful if your protein has disulfides:

python replace_cys_to_cyx.py rec.nowat.H_mod2.pdb rec.nowat.H_mod3.pdb

output may be interted into the tleap input below to draw bonds. But this is not need here.

python add.ters.py rec.nowat.H_mod3.pdb rec.nowat.H_mod4.pdb
mv rec.nowat.H_mod4.pdb rec.nowat.H_final.pdb


(3.3) produces tleap input file (past on line at a time into your terminal, be sure not to have a space before the bottom EOF.

download HEME parameters:

curl docking.org/~tbalius/code/waterpaper2017/parms/frcmod_teb_mf_mod.hemall > frcmod_teb_mf_mod.hemall
curl docking.org/~tbalius/code/waterpaper2017/parms/heme_all_plus1.in > heme_all_plus1.in
cat << EOF >! tleap.rec.in
set default PBradii mbondi2
# load the protein force field
source leaprc.ff14SB
loadamberparams frcmod_teb_mf_mod.hemall
loadamberprep heme_all_plus1.in
REC = loadpdb rec.nowat.H_final.pdb
# draw bond between HIS and HEME
# 5 coordenated iron.
bond REC.175.NE2 REC.293.FE
XTALWAT = loadpdb wat.pdb
RECWAT  = combine {REC XTALWAT}
saveamberparm RECWAT rec.leap.prm7 rec.leap.rst7
solvateBox RECWAT TIP3PBOX 10.0
saveamberparm RECWAT rec.watbox.leap.prm7 rec.watbox.leap.rst7
quit
EOF


$AMBERHOME/bin/tleap -s -f tleap.rec.in > ! tleap.rec.out

Run AMBER

(1) Equilibrate

(2) Production

(3) Restart

Run GIST post processing

Combining GIST grids

Run convergence analysis