Calculate volume of the binding site and molecules
Written by Trent Balius, Dec. 2016.
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.
- 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:
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).
- 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