Self docking study workflow with vina #chemoinformatics #vina #RDKit #pdb-tools

I posted about how to run vina from python. But I split receptor and ligand with pymol GUI at previous post, Hmm…. it’s not automated process.

I tried to write code for full auto self docking with vina. It will work only very limited option and case but It’ll be first step for Virtual screening study at home ;).

To make receptor and ligand file from pdb file, I used pdb-tools which is one of my favorite package for handing pdb.

OK, let’s write code. To use ADFR code and pdb_tools ‘subprocess’ is used. Prepare pdb file with pdb_tools step by step.

  1. pdb_selchain command selects A chain of the input pdb.
  2. pdb_delhetatm command deletes hetatom in PDB so, it will generate receptor pdb.
  3. pdb_selresname command selects used defined residue, so I used STI(imatinib) as a target ligand name was given from ArgumentParser()
  4. Most of code expect for above described step is same as previous post.
import os
import argparse
import subprocess
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.rdMolTransforms import ComputeCentroid
from vina import Vina
v = Vina(sf_name='vina')


def getParser():
    parser = argparse.ArgumentParser()
    parser.add_argument('input_pdb', help='input pdb file')
    parser.add_argument('--chain', default='A')
    parser.add_argument('--ligand', default=None)
    parser.add_argument('--output_dir', default='output')
    return parser

def check_dir(outputdir):
    if not os.path.isdir(outputdir):
        os.mkdir(outputdir)
    

def prepare_ligand_and_protein(inputfile, chain, ligand, outputdir):
    out_r_pdb = inputfile.replace('.pdb', f'_{chain}_receptor.pdb')
    out_r_pdb = os.path.join(outputdir, out_r_pdb)
    out_r_pdbqt = out_r_pdb.replace('.pdb', '.pdbqt')
    subprocess.call(f'pdb_selchain -{chain} {inputfile} | pdb_delhetatm | pdb_tidy > {out_r_pdb}', shell=True)
    subprocess.call(f'prepare_receptor -r {out_r_pdb} -o {out_r_pdbqt} -A checkhydrogens', shell=True)

    if ligand is not None:
        out_l_pdb = inputfile.replace('.pdb', f'_{chain}_ligand.pdb')
        out_l_pdb = os.path.join(outputdir, out_l_pdb)
        out_l_pdbqt = out_l_pdb.replace('.pdb', '.pdbqt')
        out_l_smi = out_l_pdb.replace('.pdb', '.smi')
        subprocess.call(f'pdb_selchain -{chain} {inputfile} | pdb_selresname -{ligand} > {out_l_pdb}', shell=True)
        subprocess.call(f'mk_prepare_ligand.py -i {out_l_pdb} -o {out_l_pdbqt} --add_hydrogen --pH 7.4', shell=True)
        return {'receptor_pdbqt':out_r_pdbqt,
                'ligand_pdbqt':out_l_pdbqt,
                'ligand_pdb':out_l_pdb,
                }
    return {'receptor_pdbqt':out_r_pdbqt}



if __name__=='__main__':
    parser = getParser()
    args = parser.parse_args()
    check_dir(args.output_dir)
    res = prepare_ligand_and_protein(args.input_pdb, args.chain, args.ligand, args.output_dir)
    if len(res) == 1:
        print('receptor PDBQT file is prepared')
        exit()

    ligand = Chem.MolFromPDBFile(res['ligand_pdb'])
    centroid = ComputeCentroid(ligand.GetConformer())
    # The all bond order of RDMol from the PDB is 1 but I think it doesn't problem to compute center of mass because it depends on location of atoms.
    # So I didnt call AssignBondOrdersFromTemplate
    print(f'Computed centroid {centroid.x:.03f}, {centroid.y:.03f}, {centroid.z:.03f}')
    v.set_receptor(res['receptor_pdbqt'])
    v.set_ligand_from_file(res['ligand_pdbqt'])
    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(os.path.join(args.output_dir, 'receptor_ligand_minimized.pdbqt'), overwrite=True)

    # Dock the ligand
    v.dock(exhaustiveness=32, n_poses=20)
    v.write_poses(os.path.join(args.output_dir, 'dock_vina_out.pdbqt'), n_poses=10, overwrite=True)

Now Ready.

Let’s test my code. It’s just type ‘python auto_dock_wf.py pdbfilename –ligand ligand name’

$ pdb_fetch 1iep > 1iep.pdb
# If ligand name isn't provided, only receptor.pdbqt file will be generated for cross docking
$ python auto_dock_wf.py 1iep.pdb --ligand STI

After waiting few min, I could get same pose ;)

Tips; This pdb file has no gap so I used the file directly but modification a/o preparation(protonation, minimization, add/remove water etc…) of input pdb file is required. To do it, other tools includes commercial packages will be required.

I updated today’s code on to my github repo. I ‘ll update and expand the code for cross docking.

https://github.com/iwatobipen/vs_vina

Advertisement

Published by iwatobipen

I'm medicinal chemist in mid size of pharmaceutical company. I love chemoinfo, cording, organic synthesis, my family.

One thought on “Self docking study workflow with vina #chemoinformatics #vina #RDKit #pdb-tools

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d bloggers like this: