Mengjuei - Blog

擬真 :: Displaying Electrostatic Surface with AMBER/PBSA (AMBER 9) and PyMol

The pyMol part of this document is based on Dr. James Stroud's howto. Therefore the license of this document also has to be in the public domain. Also please be noted that we have only optimized and tested this page on Safari and Firefox.

AMBER/PBSA, PyMol and VMD

AMBER/PBSA under AMBER 9 package

PyMol is a full featured open source molecular graphics program by Dr. Warren Delano of Delano Scientific. The PyMol website is http://pymol.sf.net. Also be sure to download recent version of pymol to make sure that no endian issue is involved (byte-order differences between MIPS, PPC and x86 machines).

VMD: We are not sure at this point why VMD's grd plugin is no working on our example. (probably due to the OpenGL issues)

This page has been tested with PyMol 0.98 on Mac OS X 10.4.3, on Irix 6.5, on Fedora Core 4 Linux and PyMol 0.96 on Windows XP.

The Instructions

1. Get a PDB file, for example, 1HVR.pdb from the PDB website to your desktop.
2. Generate Delphi formated binary electrostatics map from AMBER/PDB as 1hvr.phi, save it to your desktop.
   1. Getting rid of water molecules, hydrogen atoms and the XK2 ligand.
      You probably need this: (if you don't understand shell script, please use your favorite text editor to do the job)

      #!/bin/sh
      cat $1 | grep -v -e 'H1  WAT' -e 'H2  WAT' \
            -e 'O   WAT' -e '^ATOM ...... .H' -e '^ATOM ...... .LP' \
            -e '^ATOM ...... H.' -e '^HETATM' -e '^CONECT'
      #Please note that getting rid of CONECT is not for every case.
      #"grep -v -e" trick was provided by isomer(clive).

      % ./removeH.sh 1HVR.pdb > 1hvr_noH.pdb

      Please read the fine manual of any kind of unix to find out why.
   2. Getting prmtop/inpcrd files with LEaP, (It doesn't need any periodic boundary or box).

      % tleap -f $AMBERHOME/dat/cmd/leaprc.ff99 -f tleap.script

      where tleap.script is

      tmp = loadpdb 1hvr_noH.pdb
      saveAmberParm tmp 1hvr.top 1hvr.crd
      quit

      Please refer to your fine AMBER manual.
   3. Using SANDER's PBSA module and following input file to generate the map.

       Test of PB reaction field on protein 1HVR and output phimap
       &cntrl 
         ntx=1, irest=0,
         imin=1, ntmin=2, maxcyc=0,
         ntpr=1, igb=10, ntb=0,
         ntc=1, ntf=1, tol=0.000001,
       /   
       &pb
         npbverb=1, istrng=0, epsout=80.0, epsin=1.0, space=1.,
         accept=0.001, sprob=1.4, cutnb=9, phiout=1, phiform=0
       /
      

      Run:

      $ pbsa -O -i pbsa.in -c 1hvr.crd -p 1hvr.top -o 1hvr.out

      or

      $ sander -O -i pbsa.in -c 1hvr.crd -p 1hvr.top -o 1hvr.out

      for current version of sander in AMBER 9.
3. open the 1hvr_noH.pdb from PyMol

   PyMol> load 1hvr_noH.pdb

   or use the menu "File" -> "Open..." to locate the file (screenshot). Before you doing that, you might want to change the current directory to the path you want to store those data files, like,

   PyMol> cd /Users/mjhsieh/Desktop

   that saves you a lot of typing.
4. If you want to load-in AMBER parameter(prmtop) files and coordinate(inpcrd/restart) file, that's fine. Do this:

   PyMol> load 1hvr.top
   PyMol> load 1hvr.crd

   Please note that parameter filename has to be *.top and coordinate filename has to be *.crd for PyMol to be able to recognize it.
5. Show the surface of your molecule in the PyMol menu : click "S" and then select surface. (screenshot)
   or

   PyMol> show surface, 1HVR

6. If you want to change the radii of surface probe to some value different from 1.4Å in PyMol, e.g. 2.0Å, enter the command in the command line field.

   PyMol> set solvent_radius, 2.0

   Of course this command can be done through GUI, check the menu "Setting --> Edit All..."
7. Load the electrostatic grid in PyMol:

   PyMol> load pbsa.phi

   This grid is now an object in your PyMol menu called "pbsa". Of course you can use the menu to open this file. (screenshot, screenshot)
8. Create a color ramp in PyMol:

   PyMOL> ramp_new e_lvl, pbsa, [-7, 0, 7]

   This color ramp is now an object in your PyMol menu called "e_lvl". (screenshot)
9. Color the surface according to the grid and map:

   PyMOL> set surface_color, e_lvl, 1HVR

   (screenshot)
10. You can also adjust the color scale to [-5, 0, 5] if you like.

    PyMOL> ramp_new e_lvl, pbsa, [-5, 0, 5]

    If the color doesn't update upon the procedure, just redo

    PyMOL> set surface_color, e_lvl, 1HVR

11. In case that some setting is not taking effect, a command can be inputed to reassure.

    PyMol> rebuild

12. Raytrace in PyMol with the "ray" command (this can take a while):

    PyMOL> ray

    or press the "Ray" button.
13. Before you click on the PyMol drawing area, save the raytraced file as a .png file in PyMol with the "png" command:

    PyMOL> png my-picture.png

    or "File" -> "Save image..."

There you have it.

By mjhsieh at March 8, 2006 01:06 PM | Monthly Archives