Calculate volume of the binding site and molecules

From DISI
Revision as of 18:26, 14 September 2018 by TBalius (talk | contribs)
Jump to navigation Jump to search

Written by Trent Balius, Dec. 2016.

The method and scripts discribed here was usied in the following paper [1].

The script to calculate volume is available here: volume_cal_sph.py

syntax: python volume_cal_sph.py input.sph spacing output_prefix
example: python volume_cal_sph.py binding_site.sph 0.5 binding_site 

download using curl as follows:

curl http://docking.org/~tbalius/code/for_dock_3.7/volume_cal_sph.py > volume_cal_sph.py
curl http://docking.org/~tbalius/code/for_dock_3.7/mol2.py > mol2.py
curl http://docking.org/~tbalius/code/for_dock_3.7/sph_lib.py > sph_lib.py
curl http://docking.org/~tbalius/code/for_dock_3.7/pdb_lib.py > pdb_lib.py
curl http://docking.org/~tbalius/code/for_dock_3.7/close_sph.py > close_sph.py
curl http://docking.org/~tbalius/code/for_dock_3.7/mol2toSPH_radius.py > mol2toSPH_radius.py

how the volume calculation works.

  • First, Lay a grid over the spheres.
  • Count the number or points contained in the spheres (Ns).
  • Count the number of points in the grid box (Ng).
  • Calculate the volume of the grid box (Vb).
 Vs ~= Ns/Ng * Vb
  • This method also produces dx files, so you can visualize (using chimera or VMD)the volume as grids.

Calculating the volume of a binding site.

You can run blastermaster.py which is distributed with DOCK3.7 and then use the all_spheres.sph or lowdielectric.sph to define the pocket.

mkdir cal_vol
cp ../working/lowdielectric.sph . 
cp ../working/all_spheres.sph .

You should visualize these sphere in UCSF Chimera to make sure that they file the site and do not go outside DOCK_3.7_2014/09/25_FXa_Tutorial#Receptor_Preparation. If they do you can use a text editor to remove the necessary spheres or add new ones.


Alternavively, this command will find all ligands close to your the ligand.

python close_sph.py all_spheres.sph ../xtal-lig.pdb delphi_close.sph 2.0
head delphi_close.sph

If you do not have a ligand, you may use a central binding site residue, or place a carbon atom at the area of insterest.


Run the python program as follows:

>> python ~/zzz.scripts/volume_cal_sph.py lowdielectric_mod.sph 0.5 out

Here is the output:

input file =  lowdielectric_mod.sph
scale = 0.5
outputprefix = out
max corner =  30.49316 10.03153 8.27566
min corner =  23.27718 1.93082 -2.45935
0.5 15 17 22 0.5 0.5 0.5 [23.277180000000001, 1.9308199999999998, -2.4593500000000001]
molN= 1092   boxN= 5610   boxV= 701.25
molV= 136.5

This script also produces a dx file so that you can visualize (in chimera) the points which are overlapping with the spheres.


Alternatively, you can calculate the spheres as follows:

  • run dms (or you can also generated the molecular surface with Chimera) to generate a molecular surface.
$DOCKBASE/proteins/dms/bin/dms rec.pdb -a -g dms.log -p -n -o rec.ms
  • Use the sphgen program(distributed with all versions of DOCK) to flood the surface of the protein with spheres, which are then cluster by distance.
vi INSPH
  • This file should contain:
    • specifies the input file
    • spheres generated will be outside of teh receptor surface
    • specifies that all points won the receptor will be used
    • distance in angstroms (avoids steric clashes)
    • max surface radius of the spheres in angstroms
    • min surface radius of the spheres in angstroms
    • the specified outfile containing all generated spheres
rec.ms 
R            
X            
0.0          
4.0          
1.4          
rec.sph 
  • Run the Sphgen using the input file INSPH with the command:
$DOCKBASE/proteins/sphgen/bin/sphgen 


INSPH is input file
OUTSPH is the file containing the information about sphere genereation
rec.sph contains the spheres

Information modifed from rizzo group wiki

  • Select the cluster that defines the binding site of interest by visualization in Chimera.
    • Copy the sphere file. Using a text editor (vim) remove all clusters except the one of interest. (or see above for script that does this).

Calculating the volume of a small molecule.

Convert ligands to spheres using the following script: mol2toSPH_radius.py

Say, you would like to calculate how much two docking poses overlap in volume:

Convert the mol2 files to sph files:

python ~/zzz.scripts/mol2toSPH_radius.py molone.mol2 molone.sph
python ~/zzz.scripts/mol2toSPH_radius.py moltwo.mol2 moltwo.sph

Create a combined sphere file:

cat molone.sph > bothmol.sph
sed -e 's/cluster     1 /cluster     2 /g' moltwo.sph | grep -v DOCK >> bothmol.sph
>> python ~/zzz.scripts/volume_cal_sph.py molone.sph 0.5 molone
input file =  molone.sph
scale = 0.5
outputprefix = molone
max corner =  51.964 41.34 37.041
min corner =  39.322 31.367 28.854
0.5 26 20 17 0.5 0.5 0.5 [39.321999999999996, 31.366999999999997, 28.853999999999999]
molN= 1502   boxN= 8840   boxV= 1105.0
molV= 187.75
>> python ~/zzz.scripts/volume_cal_sph.py moltwo.sph 0.5 moltwo
input file =  moltwo.sph
scale = 0.5
outputprefix = moltwo
max corner =  49.176 38.767 36.529
min corner =  39.059 31.2 27.753
0.5 21 16 18 0.5 0.5 0.5 [39.058999999999997, 31.199999999999996, 27.753]
molN= 1257   boxN= 6048   boxV= 756.0
molV= 157.125
>> python ~/zzz.scripts/volume_cal_sph.py bothmol.sph 0.5 bothmol
input file =  bothmol.sph
scale = 0.5
outputprefix = bothmol
max corner =  51.964 41.34 37.041
min corner =  39.059 31.2 27.753
0.5 26 21 19 0.5 0.5 0.5 [39.058999999999997, 31.199999999999996, 27.753]
molN= 1868   boxN= 10374   boxV= 1296.75
molV= 233.5

The overlap region may be calculated as follows:

157.125 + 187.75 - 233.5 = 111.375