Difference between revisions of "2016 AMBER tutorial with Thrombin"

From Rizzo_Lab
Jump to: navigation, search
(MM-GBSA Energy Calculation)
Line 20: Line 20:
===MM-GBSA Energy Calculation===
===MM-GBSA Energy Calculation===
Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) is a great method to calculate or estimate relative binding affinity of a ligand(s) to a receptor. The binding energy calculated from this method are also known as free energies of binding, where the more negative values indicate stronger binding. For this section, the topology files for the ligand, receptor and complex are needed.
Create a new directory:
  mkdir 005.MMGBSA
Create an input file name
  vim gb.rescore.in
Enter the following into the input file:
  Single point GB energy calc
  ntf    = 1,        ntb    = 0,        ntc    = 2,
  idecomp= 0,
  igb    = 5,        saltcon= 0.00,
  gbsa  = 2,        surften= 1.0,
  offset = 0.09,    extdiel= 78.5,
  cut    = 99999.0,  nsnb  = 99999,
  imin  = 5,        maxcyc = 1,        ncyc  = 0,
Create a tcsh/bash/csh script (run.sander.rescore.csh) with the following information:
  #! /bin/tcsh
  #PBS -l nodes=1:ppn=1
  #PBS -l walltime=48:00:00
  #PBS -o zzz.qsub.out
  #PBS -e zzz.qsub.err
  #PBS -V
  #PBS -N mmgbsa
  set workdir = /nfs/user03/kbelfon/amber_tutorial/005.mmgbsa
  cd $workdir
  sander -O -i gb.rescore.in \
  -o gb.rescore.out.com \
  -p ../002.tleap/4TKG.com.gas.leap.prm7 \
  -c ../002.tleap/4TKG.com.gas.leap.rst7 \
  -y ../004.ptraj/4TKG.com.trj.stripfit \
  -r restrt.com \
  -ref ../002.tleap/4TKG.com.gas.leap.rst7 \
  -x mdcrd.com \
  -inf mdinfo.com
  sander -O -i gb.rescore.in \
  -o gb.rescore.out.lig \
  -p ../002.tleap/4TKG.lig.gas.leap.prm7 \
  -c ../002.tleap/4TKG.lig.gas.leap.rst7 \
  -y ../004.ptraj/4TKG.lig.trj.stripfit \
  -r restrt.lig \
  -ref ../002.tleap/4TKG.lig.gas.leap.rst7 \
  -x mdcrd.lig \
  -inf mdinfo.lig
  sander -O -i gb.rescore.in \
  -o gb.rescore.out.test.rec \
  -p ../002.tleap/4TKG.rec.gas.leap.prm7 \
  -c ../002.tleap/4TKG.rec.gas.leap.rst7 \
  -y ../004.ptraj/4TKG.rec.trj.stripfit \
  -r restrt.rec \
  -ref ../002.tleap/4TKG.rec.gas.leap.rst7 \
  -x mdcrd.rec \
  -inf mdinfo.rec
Execute this script on the seawulf cluster or machine(s) of your preference
  qsub run.sander.rescore.csh
Three output files will be generated once the job is completed:             
'''gb.rescore.out.com''', '''gb.rescore.out.lig''', and '''gb.rescore.out.rec'''
These files represent the single point energy calculation results for the complex (.com), the ligand (.lig) and the receptor (.rec). The energy will be output by the program Sander for each frame specified in the input file. The final results for one frame in one of the three files should look as the following:
                                    FINAL RESULTS
    NSTEP      ENERGY          RMS            GMAX        NAME    NUMBER
      1      5.9132E+03    2.0005E+01    1.2640E+02    C        159
  BOND    =      661.8980  ANGLE  =    1751.7992  DIHED      =    2581.7692
  VDWAALS =    -1696.6585  EEL    =  -13958.9335  EGB        =    -3125.9524
  1-4 VDW =      747.0185  1-4 EEL =    7750.8118  RESTRAINT  =        0.0000
  ESURF  =    11201.4791
minimization completed, ENE= 0.59132314E+04 RMS= 0.200047E+02
'''Extracting Data from MM-GBSA calculation and calculating Free energy of Binding'''
From the output files above:
EELS = ΔGcoul
EGB = ΔGpolar
With this information ΔGnonpolar can be solved using equation(1):
ΔGnonpolar = SASA*0.00542 + 0.92                                                      (1)
Once ΔGnonpolar is solved then ΔGmmgbsa can be determined by equation(2):
ΔGmmgbsa = ΔGvdw + ΔGcoul + ΔGpolar + ΔGnonpolar                                      (2)
Solve equation 2 and 3 using the extracted information from all three output files. So therefore you should have ΔGmmgbsa for the complex, ligand and receptor                               
Finally ΔΔGbind can be calculated using equation (3):
ΔΔGbind = ΔGmmgbsa,complex – (ΔGmmgbsa,lig + ΔGmmgbsa,rec)                       (3)
Plot your ΔΔGbind and examine the plot for changes in the ligand position and the ΔΔGbind. Also, you should calculate the mean and standard deviation for your ΔΔGbind.
The following script (get.mmgbsa.sh) can be used to extract the energy from the three output files obtained above and calculate ΔΔGbind:
  #! /bin/bash
  # by Haoquan
  echo com lig rec > namelist
  LIST=`cat namelist`
  for i in $LIST ; do
  grep VDWAALS gb.rescore.out.$i | awk '{print $3}' > $i.vdw
  grep EGB    gb.rescore.out.$i | awk '{print $9}' > $i.polar
  grep EELS    gb.rescore.out.$i | awk '{print $6}' > $i.coul
  grep ESURF  gb.rescore.out.$i | awk '{print $3 * 0.00542 + 0.92}' > $i.surf
  paste -d " " $i.vdw $i.polar $i.surf $i.coul | awk '{print $1 + $2 + $3 + $4}' > data.$i
  rm $i.*
  paste -d " " data.com data.lig data.rec | awk '{print $1 - $2 - $3}' > data.all
  for ((j=1; j<=`wc -l data.all | awk '{print $1}'`; j+=1)) do
  echo $j , >> time
  paste -d " " time data.all > MMGBSA_vs_time.dat 
  rm namelist time data.*
Execute this script:
  bash get.mmgbsa.sh
A text file called MMGBSA_vs_time.dat with x and y values separated by a space and comma should be created. Use XMGRACE to plot this dat file using the following command in Linux:
  xmgrace MMGBSA_vs_time.dat
==V. Frequently Encountered Problems==
==V. Frequently Encountered Problems==

Revision as of 14:14, 6 April 2016

For additional Rizzo Lab tutorials see AMBER Tutorials.

In this tutorial, we will learn how to run a molecular dynamics simulation of a protein-ligand complex. We will then post-process that simulation by calculating structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).

I. Introduction

II. Structural Preparation

Antechamber, Parmchk, tLeap

III. Simulation using pmemd


IV. Simulation Analysis


MM-GBSA Energy Calculation

V. Frequently Encountered Problems