A Jupyter Notebook to aid Docking to MurD protein¶
This notebook implements a typical protocol for docking ligands to a target protein. It uses RDKit (http://www.rdkit.org) to generate a number of reasonable conformations for each ligand and then uses SMINA (https://sourceforge.net/projects/smina/) to do the docking. Two methods of docking are implemented, the first docks into a rigid receptor, the second sets the protein side-chains around the active site to be flexible. Bear in mind flexible docking will be much, much slower. In the optional final step the resulting docked poses are rescored using a random forest model described in https://www.nature.com/articles/srep46710.
import sys
from collections import defaultdict
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import PandasTools
import pandas as pd
IPythonConsole.ipython_3d=True
%pylab inline
File location of structures for docking and file format¶
First we need get the location of the input file of structures you want to dock, replace "asinexSelection.sdf" with your file. You may want to rename the output file for conformations, and the output file containing the docked structures.
The sdf file needs to have the name included in the first line of each molecule record.
AEM 10028511
MOE2019 2D
22 24 0 0 0 0 0 0 0 0999 V2000 7.2040 -6.7290 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 6.3790 -6.7290 0.0000 N 0 0 0 0 0 0 0 0 0 0 0 0</code>
# File locations
sdfFilePath = 'asinexSelection.sdf' # The input file of structures to genrate conformations from
ConfoutputFilePath = 'asinexSelectionForDocking.sdf' # Output file containing conformations for docking
inputMols = [x for x in Chem.SDMolSupplier(sdfFilePath,removeHs=False)]
# Assign atomic chirality based on the structures:
len(inputMols) # Check how many strucures
#Check that all molecules have a name
for i, mol in enumerate(inputMols):
if mol is None:
print('Warning: Failed to read molecule %s in %s' % (i, sdfFilePath))
if not mol.GetProp('_Name'):
print('Warning: No name for molecule %s in %s' % (i, sdfFilePath))
Conformation generation¶
We next generate conformations, this uses paralellisation code from http://www.rdkit.org/docs/Cookbook.html contributed by Andrew Dalke. We don't use all cores on a desktop machine or it might be unresponsive. If running on a cluster you should modify this.
import multiprocessing
# Download this from http://pypi.python.org/pypi/futures
from concurrent import futures
# conda install progressbar
import progressbar
#Find number cores available, leave two or system might be unresponsive
numcores = multiprocessing.cpu_count()
max_workers = numcores -2
#Knowledge based torsion generator http://pubs.acs.org/doi/abs/10.1021/acs.jcim.5b00654
# This function is called in the subprocess.
# The parameters (molecule and number of conformers) are passed via a Python
ps = AllChem.ETKDG()
ps.pruneRmsThresh=0.5
ps.numThreads=0
#Edit for number of confs desired eg n = 5
n=5
def generateconformations(m, n, name):
m = Chem.AddHs(m)
ids=AllChem.EmbedMultipleConfs(m, n, ps)
for id in ids:
AllChem.UFFOptimizeMolecule(m, confId=id)
# EmbedMultipleConfs returns a Boost-wrapped type which
# cannot be pickled. Convert it to a Python list, which can.
return m, list(ids), name
smi_input_file, sdf_output_file = sys.argv[1:3]
writer = Chem.SDWriter(ConfoutputFilePath)
# suppl = [x for x in Chem.SDMolSupplier(sdfFilePath,removeHs=False)]
#suppl = Chem.SmilesMolSupplier(smi_input_file, titleLine=False)
# for mol in suppl:
# print(mol.GetPropsAsDict(includePrivate=True).get('_Name'))
with futures.ProcessPoolExecutor(max_workers=max_workers) as executor:
# Submit a set of asynchronous jobs
jobs = []
for mol in inputMols:
if mol:
name = mol.GetProp('_Name')
job = executor.submit(generateconformations, mol, n, name)
jobs.append(job)
widgets = ["Generating conformations; ", progressbar.Percentage(), " ",
progressbar.ETA(), " ", progressbar.Bar()]
pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(jobs))
for job in pbar(futures.as_completed(jobs)):
mol, ids, name = job.result()
mol.SetProp('_Name', name)
for id in ids:
writer.write(mol, confId=id)
writer.close()
ms = [x for x in Chem.SDMolSupplier(ConfoutputFilePath,removeHs=False)]
# Assign atomic chirality based on the structures:
for m in ms: Chem.AssignAtomChiralTagsFromStructure(m)
len(ms) # check how many conformations
Docking to Protein¶
After generating the conformations we can now do the docking. In this example we use smina which can be downloaded from https://sourceforge.net/projects/smina/ you will need to know where smina has been installed. The protein and ligand examples provided are taken from https://fragalysis.diamond.ac.uk/viewer/react/preview/target/MURD MURD-x0373.
Docking using smina
Need protein minus the ligand in pdb format,
the ligand extracted from binding site in pdb format,
Conformations to be docked as sdf from conformation generation above
DockedFilePath = 'All_Docked.sdf.gz' is the File for the Docked structures
ProteinForDocking = 'protein_minus_ligand.pdb'
LigandFromProtein = '373ligand_only.pdb'
DockedFilePath = 'All_Docked.sdf.gz'
FlexibleDockedFilePath = 'FlexDocked.sdf.gz'
!'/usr/local/bin/smina.osx' --cpu 10 --seed 0 --autobox_ligand '{LigandFromProtein}' -r '{ProteinForDocking}' -l '{ConfoutputFilePath}' -o '{DockedFilePath}'
Rescoring using Random Forest Model¶
Optional, Rescore using a random forest model described in https://www.nature.com/articles/srep46710
Download from https://github.com/oddt/rfscorevs You will need the path to the binary
Path to protein containing ligand in pdb format
protein_plus_373ligand from Diamond
File to store rescored resukts
TargetProtein = 'protein_plus_373ligand.pdb'
scoreResults = 'DockedRescored.csv'
!/usr/local/bin/rf-score-vs --receptor '{TargetProtein}' '{DockedFilePath}' -o csv -O '{scoreResults}' --field name --field RFScoreVS_v2
docked_df = PandasTools.LoadSDF(DockedFilePath,molColName='Molecule')
#docked_df.head(n=5)
scores_df = pd.read_csv(scoreResults)
#scores_df.head(n=5)
results_df = pd.concat([docked_df, scores_df], axis=1)
results_df.head(5)
| ID | Molecule | minimizedAffinity | name | RFScoreVS_v2 | |
|---|---|---|---|---|---|
| 0 | ART 15342099 |  |
-5.26767 | ART 15342099 | 6.127652 |
| 1 | ART 15342099 |  |
-5.05578 | ART 15342099 | 6.112676 |
| 2 | ART 15342099 |  |
-4.83187 | ART 15342099 | 6.085660 |
| 3 | ART 15342099 |  |
-4.74720 | ART 15342099 | 6.063921 |
| 4 | ART 15342099 |  |
-4.72279 | ART 15342099 | 6.109961 |
Saving the results¶
Now combine rescored file with docked structure file and export to "Alldata.sdf.gz" this is a big file so export compressed
combinedResults = 'Alldata.sdf.gz'
PandasTools.WriteSDF(results_df, combinedResults, molColName="Molecule", idName="ID", properties=list(results_df.columns))