Calculate volume of the binding site and molecules

Written by Trent Balius, Dec. 2016.

The method and scripts discribed here was usied in this 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.

Alternatively, 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 the receptor surface with R and inside with L
  • specifies that all points on the receptor will be used
  • distance in angstroms (0.0 avoids steric clashes), try -0.1 to completely fill the site, but some 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 modified 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).

• Calculate the volume using the following script: volume_cal_sph.py

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