2.1: Single Molecule#
Authors: Mateusz K Bieniek, Ben Cree, Rachael Pirie, Joshua T. Horton, Natalie J. Tatum, Daniel J. Cole
- Add chemical functional groups (R-groups) in user-defined positions
- Output ADMET properties
- Perform constrained optimisation
- Score poses
- Output structures for free energy calculations
Overview#
This notebook demonstrates the entire FEgrow workflow for generating a series of ligands with a common core for a specific binding site, via the addition of a user-defined set of R-groups.
These de novo ligands are then subjected to ADMET analysis. Valid conformers of the added R-groups are enumerated, and optimised in the context of the receptor binding pocket, optionally using hybrid machine learning / molecular mechanics potentials (ML/MM).
An ensemble of low energy conformers is generated for each ligand, and scored using the gnina convolutional neural network (CNN). Output structures are saved as pdb files ready for use in free energy calculations.
The target for this tutorial is the main protease (Mpro) of SARS-CoV-2, and the core and receptor structures are taken from a recent study by Jorgensen & co-workers.
import prody
from rdkit import Chem
import fegrow
from fegrow import RGroups, Linkers, ChemSpace
# Initialise libraries
rgroups = RGroups()
linkers = Linkers()
<frozen importlib._bootstrap>:241: RuntimeWarning: to-Python converter for boost::shared_ptr<RDKit::FilterHierarchyMatcher> already registered; second conversion method ignored.
<frozen importlib._bootstrap>:241: RuntimeWarning: to-Python converter for boost::shared_ptr<RDKit::FilterCatalogEntry> already registered; second conversion method ignored.
MolGridWidget(grid_id='m2')
MolGridWidget(grid_id='m1')
Prepare the ligand template#
The provided core structure lig.pdb has been extracted from a crystal structure of Mpro in complex with compound 4 from the Jorgensen study (PDB: 7L10), and a Cl atom has been removed to allow growth into the S3/S4 pocket. The template structure of the ligand is protonated with Open Babel:
!obabel sarscov2/lig.pdb -O sarscov2/coreh.sdf -p 7
1 molecule converted
Load the protonated ligand into FEgrow:
init_mol = Chem.SDMolSupplier('sarscov2/mini.sdf', removeHs=False)[0]
# get the FEgrow representation of the rdkit Mol
scaffold = fegrow.RMol(init_mol)
Show the 2D (with indices) representation of the core. This is used to select the desired growth vector.
scaffold.rep2D(idx=True, size=(500, 500))
Using the 2D drawing, select an index for the growth vector. Note that it is currently only possible to grow from hydrogen atom positions. In this case, we are selecting the hydrogen atom labelled H:40 to enable growth into the S3/S4 pocket of Mpro.
attachment_index = 8
Optional: insert a linker#
We have added a library of linkers suggested by Ertl et al. If you wish to extend your R groups selection via a linker, select them below. :1 is defined to be attached to the core (there exists a mirror image of each linker i.e. :1 & *:2 swapped).
Linkers combinatorially augment chosen R groups, so if you choose 2 linkers and 3 R groups, this will result in 6 molecules being built.
Note : If you want to use linkers make sure that you use the correct function below, in cell [11].#
linkers
l = linkers[linkers.Name == 'R1CR2'].Mol.item()
l
| SmileIndex | 0 |
|---|---|
| ID | |
| SMILES | [H]C([H])([*:1])[*:2] |
| display_smiles | R1CR2 |
# get linkers programmatically
rcr_linker = linkers[linkers.Name == 'R1CR2'].Mol.item()
rocr_linker = linkers.Mol[6], # use the linker table index directly, e.g. index 6 is "R2COR1"
# or pick linkers from the grid
grid_linkers = linkers.get_selected()
# select only one
# grid_linker = grid_linkers.pop()
# create one new molecule merged with a linker
scaffold_with_linker = fegrow.build_molecule(scaffold, rcr_linker, attachment_index)
# visualise:
# note that the linker leaves the second attachement point prespecified (* character)
scaffold_with_linker.rep2D(idx=True, size=(500, 500))
scaffold_with_linker
| Smiles | Molecule | |
|---|---|---|
| ID | ||
| None | *C([H])([H])c1c([H])nc([H])c([H])c1[H] |  |
Select RGroups for your template#
R-groups can be selected interactively or programmaticaly.
We have provided a set of common R-groups (see fegrow/data/rgroups/library), which can be browsed and selected interactively below.
Molecules from the library can alternatively be selected by name, as demonstrated below.
Finally, user-defined R-groups may be provided as .mol files. In this case, the hydrogen atom selected for attachment should be replaced by the element symbol R. See the directory manual_rgroups for examples.
rgroups
# retrieve the first interactively selected group
interactive_rgroups = rgroups.get_selected()
# interactive_rgroup = interactive_rgroups.pop()
# you can also directly access the built-in dataframe programmatically
R_group_ethanol = rgroups[rgroups.Name == '*CCO'].Mol.item()
# select the R-group using the index
R_group_cyclopropane = rgroups.Mol[69]
# use SMILES
R_group_methanol = Chem.AddHs(Chem.MolFromSmiles('*CO'))
# add your R-groups from files
R_group_propanol = Chem.MolFromMolFile('manual_rgroups/propan-1-ol-r.mol', removeHs=False)
Build a congeneric series#
Now that the R-groups have been selected, we merge them with the ligand core:
Note : Use rmols = fegrow.build_molecules(template_with_linker, selected_rgroups) if using a linker, by commenting out the first fegrow.build_molecules function.#
# we can either use the original template (so no linker)
# in this case we have to specify the attachment index
# rmol = fegrow.build_molecule(scaffold, R_group_methanol, attachment_index)
# or we can use the template merged with the linker
# in which case the attachement point is not needed (R* atom is used)
rmol = fegrow.build_molecule(scaffold_with_linker, R_group_ethanol)
rmol
| Smiles | Molecule | |
|---|---|---|
| ID | ||
| None | [H]OC([H])([H])C([H])([H])C([H])([H])c1c([H])n... |  |
rmol
| Smiles | Molecule | |
|---|---|---|
| ID | ||
| None | [H]OC([H])([H])C([H])([H])C([H])([H])c1c([H])n... | |
The R-group library can also be viewed as a 2D grid, or individual molecules can be selected for 3D view (note that the conformation of the R-group has not yet been optimised):
rmol.rep2D()
rmol.rep3D()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
<py3Dmol.view at 0x797251865b10>
Once the ligands have been generated, they can be assessed for various ADMET properties, including Lipinksi rule of 5 properties, the presence of unwanted substructures or problematic functional groups, and synthetic accessibility.
rmol.toxicity()
[11:28:55] DEPRECATION WARNING: please use MorganGenerator
| MW | HBA | HBD | LogP | Pass_Ro5 | has_pains | has_unwanted_subs | has_prob_fgs | synthetic_accessibility | Smiles | |
|---|---|---|---|---|---|---|---|---|---|---|
| ID | ||||||||||
| None | 137.182 | 2 | 1 | 1.007 | True | False | False | False | 7.71 | [H]OC([H])([H])C([H])([H])C([H])([H])c1c([H])n... |
For each ligand, a specified number of conformers (num_conf) is generated by using the RDKit ETKDG algorithm. Conformers that are too similar to an existing structure are discarded. Empirically, we have found that num_conf=200 gives an exhaustive search, and num_conf=50 gives a reasonable, fast search, in most cases.
If required, a third argument can be added flexible=[0,1,...], which provides a list of additional atoms in the core that are allowed to be flexible. This is useful, for example, if growing from a methyl group and you would like the added R-group to freely rotate.
rmol.generate_conformers(num_conf=50,
minimum_conf_rms=0.5,
# flexible=[3, 18, 20])
)
Generated 16 conformers.
Prepare the protein#
The protein-ligand complex structure is downloaded, and PDBFixer is used to protonate the protein, and perform other simple repair:
# get the protein-ligand complex structure
!wget -nc https://files.rcsb.org/download/7L10.pdb
# load the complex with the ligand
sys = prody.parsePDB('7L10.pdb')
# remove any unwanted molecules
rec = sys.select('not (nucleic or hetatm or water)')
# save the processed protein
prody.writePDB('rec.pdb', rec)
# fix the receptor file (missing residues, protonation, etc)
fegrow.fix_receptor("rec.pdb", "rec_final.pdb")
# load back into prody
rec_final = prody.parsePDB("rec_final.pdb")
File ‘7L10.pdb’ already there; not retrieving.
@> 2609 atoms and 1 coordinate set(s) were parsed in 0.02s.
@> 4638 atoms and 1 coordinate set(s) were parsed in 0.03s.
View enumerated conformers in complex with protein:
rmol.rep3D(prody=rec_final)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
<py3Dmol.view at 0x79724a4e0b10>
Any conformers that clash with the protein (any atom-atom distance less than 1 Angstrom), are removed.
rmol.remove_clashing_confs(rec_final)
Removed 7 conformers.
| Smiles | Molecule | |
|---|---|---|
| ID | ||
| None | [H]OC([H])([H])C([H])([H])C([H])([H])c1c([H])n... | |
rmol.rep3D(prody=rec_final)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
<py3Dmol.view at 0x79724aed9690>
Optimise conformers in context of protein#
The remaining conformers are optimised using hybrid machine learning / molecular mechanics (ML/MM), using the ANI2x neural nework potential for the ligand energetics (as long as it contains only the atoms H, C, N, O, F, S, Cl). Note that the Open Force Field Parsley force field is used for intermolecular interactions with the receptor.
sigma_scale_factor: is used to scale the Lennard-Jones radii of the atoms.
relative_permittivity: is used to scale the electrostatic interactions with the protein.
water_model: can be used to set the force field for any water molecules present in the binding site.
# opt_mol, energies
energies = rmol.optimise_in_receptor(
receptor_file="rec_final.pdb",
ligand_force_field="openff",
use_ani=True,
sigma_scale_factor=0.8,
relative_permittivity=4,
water_model = None,
platform_name = 'CPU', # or e.g. 'CUDA'
)
/home/dresio/software/mambaforge/envs/fegrow-onechannel/lib/python3.11/site-packages/parmed/structure.py:1799: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
coords = np.asanyarray(value, dtype=np.float64)
using ani2x
/home/dresio/software/mambaforge/envs/fegrow-onechannel/lib/python3.11/site-packages/torchani/aev.py:16: UserWarning: cuaev not installed
warnings.warn("cuaev not installed")
/home/dresio/software/mambaforge/envs/fegrow-onechannel/lib/python3.11/site-packages/torchani/__init__.py:59: UserWarning: Dependency not satisfied, torchani.ase will not be available
warnings.warn("Dependency not satisfied, torchani.ase will not be available")
Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.
/home/dresio/software/mambaforge/envs/fegrow-onechannel/lib/python3.11/site-packages/torchani/resources/
failed to equip `nnpops` with error: No module named 'NNPOps'
Optimising conformer: 100%|███████████████████████| 9/9 [00:28<00:00, 3.13s/it]
The rmol might have no available conformers due to unresolvable steric clashes with the protein. This can be checked using the RDKit's function:
rmol.GetNumConformers()
7
Optionally, display the final optimised conformers. Note that, unlike classical force fields, ANI allows bond breaking. You may occasionally see ligands with distorted structures and very long bonds, but in our experience these are rarely amongst the low energy structures and can be ignored.
Conformers are now sorted by energy, only retaining those within 5 kcal/mol of the lowest energy structure:
final_energies = rmol.sort_conformers(energy_range=5)
rmol.rep3D()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
<py3Dmol.view at 0x797251879590>
Save all of the lowest energy conformers to files and print the sorted energies in kcal/mol (shifted so that the lowest energy conformer is zero).
rmol.to_file(f"rmol_best_conformers.pdb")
print(final_energies)
ID Conformer Energy
0 None 0 0.000
1 None 1 3.274
2 None 2 4.162
3 None 3 4.320
The conformers are scored using the Gnina molecular docking program and convolutional neural network scoring function. [Note that this step is not supported on macOS]. If unavailable, the Gnina executable is downloaded during the first time it is used. The CNNscores may also be converted to predicted Kd (nM) (see column "Kd").
affinities = rmol.gnina(receptor_file="rec_final.pdb")
affinities
| ID | Conformer | CNNaffinity | Kd | |
|---|---|---|---|---|
| 0 | 0 | 0 | 3.190 | 645505.578769006 |
| 1 | 0 | 1 | 3.183 | 655918.6808887193 |
| 2 | 0 | 2 | 3.256 | 554549.0938767607 |
| 3 | 0 | 3 | 3.193 | 640722.5365758869 |
Predicted binding affinities may be further refined using the structures output by FEgrow, using your favourite free energy calculation engine. See our paper for an example using SOMD to calculate the relative binding free energies of 13 Mpro inhibitors.
# display units
affinities.Kd
0 645505.578769006
1 655918.6808887193
2 554549.0938767607
3 640722.5365758869
Name: Kd, dtype: pint[nanomolar][Float64]