Run docking study from python #chemoinformatics #vina #RDKit

Docking is one of the popular approach for computer aided drug design. There are lots of applications to run docking not only commercial software but also open source.

AutoDock Vina is one of the popular OSS for docking study and it was updated recently. The publication is below. https://pubs.acs.org/doi/full/10.1021/acs.jcim.1c00203

Vina has python binding and it can install with conda command. So I tried to use vina from python for self docking.

At first, I downloaded ADFRSuite and install it.

$ tar zxvf ADFRsuite_x86_64Linux_1.0.tar.gz
$ cd ADFRsuite_x86_64Linux_1.0
$ mkdir /home/iwatobipen/src/ADFRSuite1.0
$ ./install.sh -d /home/iwatobipen/src/ADFRSuite1.0 -c 0
$ export PATH=/home/iwatobipen/src/ADFR1.0/bin

Then I installed vina with following command.

$ mamba install -c ccsb-scripps vina

Now ready. To run the docking study, input files (receptor and ligand) are required as pdbqt format.

I get 1iep from PDB and saved receptor and ligand separately as 1iep_receptorH.pdb and 1iep_ligandH.sdf. There are many way to do above and I did it with pymol GUI.

Then convert pdb and sdf to pdbqt format with following command.

$ prepare_receptor -r 1iep_receptorH.pdb -o 1iep_receptor.pdbqt
$ mk_prepare_ligand.py -i 1iep_ligand.sdf -o 1iep_ligand.pdbqt --add_hydrogen --pH 7.4

After these step, I could run docking study on jupyter-notebook. To run vina dock, I need to prepare map which defines ligand-receptor binding site. To define the box, centroid of the ligand is required and get the point with rdkit.

Whole code is below.

from rdkit import Chem
from rdkit.Chem.rdMolTransforms import ComputeCentroid
from vina import Vina
v = Vina(sf_name='vina')
lig = Chem.SDMolSupplier('../1iep_ligandH.sdf', removeHs=False)[0]

#RDKit can't read pdbqt format so I use SDF for centroid calculation
centroid = ComputeCentroid(lig.GetConformer())

v.set_receptor('../generatedfile/1iep_receptor.pdbqt')
v.set_ligand_from_file('../generatedfile/1iep_ligand.pdbqt')

print(centroid.x, centroid.y, centroid.z)
> 15.6138918918919 53.38013513513513 15.45483783783784

v.compute_vina_maps(center=[centroid.x, centroid.y, centroid.z], box_size=[20, 20, 20])

energy_minimized = v.optimize()
print('Score after minimization : %.3f (kcal/mol)' % energy_minimized[0])
v.write_pose('1iep_ligand_minimized.pdbqt', overwrite=True)

# Dock the ligand
v.dock(exhaustiveness=32, n_poses=20)
v.write_poses('1iep_ligand_vina_out.pdbqt', n_poses=10, overwrite=True)

It’s basic docking code but interesting for me because now I can run docking study just using python only. It means that it’s easy to make automated docking workflow with python.

Here is a results. Pink is pose from PDB and green is pose from vina_docking. I showed top 2 pose. Following result shows vina proposed docking pose which is close to binding pose. Sure I know, pose prediction is still very difficult task and docking is not so accurate in many case. But it worth to know how to do it. ;)