Difference between revisions of "2019 AMBER tutorial with PDBID 2BXF"
From Rizzo_Lab
Stonybrook (talk | contribs) (→TLeap) |
Stonybrook (talk | contribs) (→TLeap) |
||
| Line 41: | Line 41: | ||
###Needed so we can use igb=8 model | ###Needed so we can use igb=8 model | ||
set default PBradii mbondi3 | set default PBradii mbondi3 | ||
| − | |||
###Load Protein pdb file | ###Load Protein pdb file | ||
Revision as of 15:03, 19 April 2019
2BXF with an explicit solvent model
Prepare the files
Convert 2BXF.lig.withH.charged.mol2 to pdb in chimera
Convert 2BXF.rec.withH.charged.mol2 to pdb in chimera
Copy into zzz.master
Parameters
The system we are working on has two main components (Protein receptor & ligand). The usual forcefield "ff14SB" contains all the parameters needed for calculations of the protein. However, the ligand is a non-protein component. Therefore, ff14SB forcefield does not contain the parameters needed for the calculations regarding the ligand. Therefore, we need to generate parameters needed for the ligand. The following steps will be taken using antechamber in order to generate the ligand parameters.
Move into 000.programs
Paramaterize the ligand
antechamber -i ../zzz.master/2BXF.lig.withH.charged.pdb -fi pdb -o 2BXF_lig.am1bcc.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc -1
Check for missing force field parameters
parmchk2 -i 2BXF_lig.am1bcc.mol2 -f mol2 -o 2BXF_lig.am1bcc.frcmod
TLeap
Move into 001.tleap_build
Create tleap.build.in file
#!/usr/bin/sh
###Load Protein force field source leaprc.protein.ff14SB ###Load GAFF force field (for our ligand) source leaprc.gaff ###Load TIP3P (water) force field source leaprc.water.tip3p ####Load Ions frcmod for the tip3p model loadamberparams frcmod.ionsjc_tip3p ###Needed so we can use igb=8 model set default PBradii mbondi3
###Load Protein pdb file rec=loadpdb ../zzz.master/backup.pdb bond rec.49.SG rec.58.SG bond rec.241.SG rec.249.SG bond rec.71.SG rec.87.SG bond rec.86.SG rec.97.SG bond rec.120.SG rec.165.SG bond rec.164.SG rec.173.SG bond rec.510.SG rec.555.SG bond rec.472.SG rec.482.SG bond rec.457.SG rec.473.SG bond rec.196.SG rec.242.SG bond rec.261.SG rec.275.SG bond rec.312.SG rec.357.SG bond rec.274.SG rec.285.SG bond rec.388.SG rec.434.SG bond rec.433.SG rec.444.SG bond rec.554.SG rec.563.SG bond rec.356.SG rec.365.SG
###Load Ligand frcmod/mol2 loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
###Create gas-phase complex
gascomplex= combine {rec lig}
###Write gas-phase pdb savepdb gascomplex 2BXF.gas.complex.pdb
###Write gas-phase toplogy and coord files for MMGBSA calc saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7 saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7 saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
###Create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}
###Solvate the system solvateoct solvcomplex TIP3PBOX 12.0
###Neutralize system (it will add either Na or Cl depending on net charge) addions solvcomplex Cl- 0 addions solvcomplex Na+ 0
###Write solvated pdb file savepdb solvcomplex 2BXF.wet.complex.pdb
###Check the system charge solvcomplex check solvcomplex
###Write Solvated topology and coord file
saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
Create Amber topology and coordinates files for the MD simulation
tleap -f tleap.build.in
#!/usr/bin/sh
###Load Protein force field source leaprc.protein.ff14SB ###Load GAFF force field (for our ligand) source leaprc.gaff ###Load TIP3P (water) force field source leaprc.water.tip3p ####Load Ions frcmod for the tip3p model loadamberparams frcmod.ionsjc_tip3p ###Needed so we can use igb=8 model set default PBradii mbondi3
###Load Protein pdb file rec=loadpdb ../zzz.master/backup.pdb bond rec.49.SG rec.58.SG bond rec.241.SG rec.249.SG bond rec.71.SG rec.87.SG bond rec.86.SG rec.97.SG bond rec.120.SG rec.165.SG bond rec.164.SG rec.173.SG bond rec.510.SG rec.555.SG bond rec.472.SG rec.482.SG bond rec.457.SG rec.473.SG bond rec.196.SG rec.242.SG bond rec.261.SG rec.275.SG bond rec.312.SG rec.357.SG bond rec.274.SG rec.285.SG bond rec.388.SG rec.434.SG bond rec.433.SG rec.444.SG bond rec.554.SG rec.563.SG bond rec.356.SG rec.365.SG
###Load Ligand frcmod/mol2 loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
###Create gas-phase complex
gascomplex= combine {rec lig}
###Write gas-phase pdb savepdb gascomplex 2BXF.gas.complex.pdb
###Write gas-phase toplogy and coord files for MMGBSA calc saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7 saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7 saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
###Create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}
###Solvate the system solvateoct solvcomplex TIP3PBOX 12.0
###Neutralize system (it will add either Na or Cl depending on net charge) addions solvcomplex Cl- 0 addions solvcomplex Na+ 0
###Write solvated pdb file savepdb solvcomplex 2BXF.wet.complex.pdb
###Check the system charge solvcomplex check solvcomplex
###Write Solvated topology and coord file saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
#!/usr/bin/sh
###Load Protein force field source leaprc.protein.ff14SB ###Load GAFF force field (for our ligand) source leaprc.gaff ###Load TIP3P (water) force field source leaprc.water.tip3p ####Load Ions frcmod for the tip3p model loadamberparams frcmod.ionsjc_tip3p ###Needed so we can use igb=8 model set default PBradii mbondi3
###Load Protein pdb file rec=loadpdb ../zzz.master/backup.pdb bond rec.49.SG rec.58.SG bond rec.241.SG rec.249.SG bond rec.71.SG rec.87.SG bond rec.86.SG rec.97.SG bond rec.120.SG rec.165.SG bond rec.164.SG rec.173.SG bond rec.510.SG rec.555.SG bond rec.472.SG rec.482.SG bond rec.457.SG rec.473.SG bond rec.196.SG rec.242.SG bond rec.261.SG rec.275.SG bond rec.312.SG rec.357.SG bond rec.274.SG rec.285.SG bond rec.388.SG rec.434.SG bond rec.433.SG rec.444.SG bond rec.554.SG rec.563.SG bond rec.356.SG rec.365.SG
###Load Ligand frcmod/mol2 loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
###Create gas-phase complex
gascomplex= combine {rec lig}
###Write gas-phase pdb savepdb gascomplex 2BXF.gas.complex.pdb
###Write gas-phase toplogy and coord files for MMGBSA calc saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7 saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7 saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
###Create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}
###Solvate the system solvateoct solvcomplex TIP3PBOX 12.0
###Neutralize system (it will add either Na or Cl depending on net charge) addions solvcomplex Cl- 0 addions solvcomplex Na+ 0
###Write solvated pdb file savepdb solvcomplex 2BXF.wet.complex.pdb
###Check the system charge solvcomplex check solvcomplex
###Write Solvated topology and coord file saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
Analysis
RMSD
How far our ligand moved and how far our receptor moved
Create cpptraj.strip.wat.in
#!/usr/bin/sh parm ../../001.tleap_build/2BXF.wet.complex.prmtop #read in trajectory trajin ../../003.production/001.restrained/10.prod.trj #read in reference reference ../../001.tleap_build/2BXF.wet.complex.rst7 #compute rmsd and align CA to the crystal structure rmsd rms1 reference :1-131@CA #strip Solvent strip :WAT:Na+:Cl- #create gas-phase trajectory trajout 2BXF.stripfit.restrained.gas.trj nobox
Run cpptraj
cpptraj -i cpptraj.strip.wat.in
Create cpptraj.rmsd.lig.in
#!/usr/bin/sh #trajin the trajectory trajin 2BXF.stripfit.restrained.gas.trj #read in the reference reference ../../001.tleap_build/2BXF.gas.complex.rst7 #compute the RMSD (do not fit the internal geometries first, included rigid body motions #and convert the frames to ns (framenum*.005) rmsd rms1 ":132&!(@H=)" nofit mass out 2BXF.lig.restrained.rmsd.nofit.dat time .005 #histogram the nofit rmsd histogram rms1,*,*,.1,* norm out 2BXF.lig.restrained.rmsd.nofit.histogram.dat
Run cpptraj
cpptraj -p ../../001.tleap_build/2BXF.gas.complex.prmtop -i cpptraj.rmsd.lig.in
Create cpptraj.rmsd.rec.in
#!/usr/bin/sh #trajin the trajectory trajin 2BXF.stripfit.restrained.gas.trj #read in the reference reference ../../001.tleap_build/2BXF.gas.complex.rst7 #compute the RMSD (do not fit the internal geometries first, included rigid body motions #and convert the frames to ns (framenum*.005) rmsd rms1 ":1-131&!(@H=)" nofit mass out 2BXF.rec.restrained.rmsd.nofit.dat time .005 #histogram the nofit rmsd histogram rms1,*,*,.1,* norm out 2BXF.rec.restrained.rmsd.nofit.histogram.dat
Run cpptraj
cpptraj -p ../../001.tleap_build/2BXF.gas.complex.prmtop -i cpptraj.rmsd.rec.in
H Bond
Hydrogen bonding
Create cpptraj.hbond.in
#!/usr/bin/sh #read in trajectory trajin ../../003.production/001.restrained/10.prod.trj #wrap everything into one periodic cell #autoimage #compute intra and water mediated hydrogen bonds hbond hb1 :1-288 out 2BXF_sunitinib.hbond.out avgout 2BXF_sunitinib.hbond.avg.dat solventdonor :WAT solventacceptor :WAT@O nointramol brid\ geout 2BXF_sunitinib.bridge-water.dat dist 3.0 angle 140
Run cpptraj
cpptraj -p ../../001.tleap_build/2BXF.wet.complex.prmtop -i cpptraj.hbond.in
MMGBSA
This will analyze how strongly our small molecule binds to our receptor.
Create mmgbsa.in
mmgbsa 2BXF analysis
&general
interval=1, netcdf=1,
keep_files=0,
/
&gb
igb=8,
saltcon=0.0, surften=0.0072,
surfoff=0.0, molsurf=0,
/
&nmode
drms=0.001, maxcyc=10000,
nminterval=250, nmendframe=2000,
nmode_igb=1,
/
Create mmgbsa.sh
#!/bin/bash
#PBS -l walltime=35:00:00
#PBS -l nodes=1:ppn=24
#PBS -N 4qmz_mmgbsa
#PBS -V
#PBS -j oe
#PBS -q long-24core
cd $PBS_O_WORKDIR
#Define topology files
solv_prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
complex_prmtop="../../001.tleap_build/2BXF.gas.complex.prmtop"
receptor_prmtop="../../001.tleap_build/2BXF.gas.receptor.prmtop"
ligand_prmtop="../../001.tleap_build/2BXF.gas.ligand.prmtop"
trajectory="../../003.production/001.restrained/10.prod.trj"
#create mmgbsa input file
cat >mmgbsa.in<<EOF
mmgbsa HIVgp41 analysis
&general
interval=1, netcdf=1,
keep_files=0,
/
&gb
igb=8,
saltcon=0.0, surften=0.0072,
surfoff=0.0, molsurf=0,
/
&nmode
drms=0.001, maxcyc=10000,
nminterval=250, nmendframe=2000,
nmode_igb=1,
/
EOF
MMPBSA.py -O -i mmgbsa.in \
-o 2BXF.mmgbsa.results.dat \
-eo 2BXF.mmgbsa.per-frame.dat \
-sp ${solv_prmtop} \
-cp ${complex_prmtop} \
-rp ${receptor_prmtop} \
-lp ${ligand_prmtop} \
-y ${trajectory}
Submit your script
qsub mmgbsa.sh
The last line of our 2BXF.mmgbsa.results.dat file shows that our delta G binding is -24.4431 +/- 5.4851.
