Calculate volume of the binding site and molecules

From DISI
Jump to navigation Jump to search

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

Calculating the volume of a binding site using scripts

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