Difference between revisions of "2019 AMBER tutorial with PDBID 2BXF"

From Rizzo_Lab
Jump to: navigation, search
(→‎Production)
(→‎Analysis)
Line 230: Line 230:
 
  saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
 
  saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
  
−
= Analysis =
+
  #!/usr/bin/sh
−
 
+
  ###Load Protein force field
−
== RMSD ==
+
source leaprc.protein.ff14SB
−
How far our ligand moved and how far our receptor moved
+
  ###Load GAFF force field (for our ligand)
−
 
+
  source leaprc.gaff
−
Create cpptraj.strip.wat.in
+
  ###Load TIP3P (water) force field
−
 
+
  source leaprc.water.tip3p
−
  #!/usr/bin/sh                                                              
+
  ####Load Ions frcmod for the tip3p model
−
  parm ../../001.tleap_build/2BXF.wet.complex.prmtop
+
  loadamberparams frcmod.ionsjc_tip3p
−
  #read in trajectory                                                       
+
  ###Needed so we can use igb=8 model
−
  trajin ../../003.production/001.restrained/10.prod.trj
+
  set default PBradii mbondi3
−
  #read in reference                                                         
+
  ###Load Protein pdb file
−
  reference ../../001.tleap_build/2BXF.wet.complex.rst7
+
  rec=loadpdb ../zzz.master/backup.pdb
−
  #compute rmsd and align CA to the crystal structure                       
+
  bond rec.49.SG  rec.58.SG
−
  rmsd rms1 reference :1-131@CA
+
  bond rec.241.SG rec.249.SG
−
  #strip Solvent                                                             
+
  bond rec.71.SG  rec.87.SG
−
  strip :WAT:Na+:Cl-
+
bond rec.86.SG  rec.97.SG
−
  #create gas-phase trajectory                                               
+
  bond rec.120.SG rec.165.SG
−
  trajout 2BXF.stripfit.restrained.gas.trj nobox
+
  bond rec.164.SG rec.173.SG
−
 
+
  bond rec.510.SG rec.555.SG
−
Run cpptraj
+
bond rec.472.SG rec.482.SG
−
  cpptraj -i cpptraj.strip.wat.in
+
bond rec.457.SG rec.473.SG
−
 
+
  bond rec.196.SG rec.242.SG
−
Create cpptraj.rmsd.lig.in
+
  bond rec.261.SG rec.275.SG
−
  #!/usr/bin/sh                                                             
+
bond rec.312.SG rec.357.SG
−
#trajin the trajectory                                                     
+
  bond rec.274.SG rec.285.SG
−
trajin 2BXF.stripfit.restrained.gas.trj
+
  bond rec.388.SG rec.434.SG
−
  #read in the reference                                                     
+
  bond rec.433.SG rec.444.SG
−
reference ../../001.tleap_build/2BXF.gas.complex.rst7
+
bond rec.554.SG rec.563.SG
−
  #compute the RMSD (do not fit the internal geometries first, included rigid body motions                                                             
+
  bond rec.356.SG rec.365.SG
−
#and convert the frames to ns (framenum*.005)                             
+
  ###Load Ligand frcmod/mol2
−
rmsd rms1 ":132&!(@H=)" nofit mass out 2BXF.lig.restrained.rmsd.nofit.dat time .005
+
  loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod
−
  #histogram the nofit rmsd                                                 
+
lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
−
histogram rms1,*,*,.1,* norm out 2BXF.lig.restrained.rmsd.nofit.histogram.dat
+
  ###Create gas-phase complex
−
 
+
  gascomplex= combine {rec lig}
−
 
+
  ###Write gas-phase pdb
−
Run cpptraj
+
  savepdb gascomplex 2BXF.gas.complex.pdb
−
  cpptraj -p ../../001.tleap_build/2BXF.gas.complex.prmtop -i cpptraj.rmsd.lig.in
+
  ###Write gas-phase toplogy and coord files for MMGBSA calc
−
 
+
  saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7
−
Create cpptraj.rmsd.rec.in
+
  saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7
−
  #!/usr/bin/sh                                                                                   
+
  saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
−
#trajin the trajectory                                                                         
+
  ###Create solvated complex (albeit redundant)
−
trajin 2BXF.stripfit.restrained.gas.trj
+
  solvcomplex= combine {rec lig}
−
  #read in the reference                                                                         
+
  ###Solvate the system
−
reference ../../001.tleap_build/2BXF.gas.complex.rst7
+
  solvateoct solvcomplex TIP3PBOX 12.0
−
  #compute the RMSD (do not fit the internal geometries first, included rigid body motions       
+
  ###Neutralize system (it will add either Na or Cl depending on net charge)
−
#and convert the frames to ns (framenum*.005)                                                   
+
  addions solvcomplex Cl- 0
−
rmsd rms1 ":1-131&!(@H=)" nofit mass out 2BXF.rec.restrained.rmsd.nofit.dat time .005
+
  addions solvcomplex Na+ 0
−
  #histogram the nofit rmsd                                                                       
+
  ###Write solvated pdb file
−
histogram rms1,*,*,.1,* norm out 2BXF.rec.restrained.rmsd.nofit.histogram.dat
+
  savepdb solvcomplex 2BXF.wet.complex.pdb
−
 
+
###Check the system
−
Run cpptraj
+
charge solvcomplex
−
  cpptraj -p ../../001.tleap_build/2BXF.gas.complex.prmtop -i cpptraj.rmsd.rec.in
+
check solvcomplex
−
 
+
  ###Write Solvated topology and coord file
−
==H Bond==
+
  saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
−
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.
 

Revision as of 15:01, 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
  1. !/usr/bin/sh
      1. Load Protein force field

source leaprc.protein.ff14SB

      1. Load GAFF force field (for our ligand)

source leaprc.gaff

      1. Load TIP3P (water) force field

source leaprc.water.tip3p

        1. Load Ions frcmod for the tip3p model

loadamberparams frcmod.ionsjc_tip3p

      1. Needed so we can use igb=8 model

set default PBradii mbondi3


      1. 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

      1. Load Ligand frcmod/mol2

loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2

      1. Create gas-phase complex

gascomplex= combine {rec lig}

      1. Write gas-phase pdb

savepdb gascomplex 2BXF.gas.complex.pdb

      1. 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

      1. Create solvated complex (albeit redundant)

solvcomplex= combine {rec lig}

      1. Solvate the system

solvateoct solvcomplex TIP3PBOX 12.0

      1. Neutralize system (it will add either Na or Cl depending on net charge)

addions solvcomplex Cl- 0 addions solvcomplex Na+ 0

      1. Write solvated pdb file

savepdb solvcomplex 2BXF.wet.complex.pdb

      1. Check the system

charge solvcomplex check solvcomplex

      1. 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

#!/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
Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=2019_AMBER_tutorial_with_PDBID_2BXF&oldid=11802"