1: Basic Full Tutorial#
Authors: Mateusz K Bieniek, Ben Cree, Rachael Pirie, Joshua T. Horton, Natalie J. Tatum, Daniel J. Cole
- Add chemical functional group (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 simplified FEgrow workflow for generating a new molecule with a common core for a specific binding site, via the addition of a user-defined R-group.
This notebook introduces FEgrow conceptually but better tools are introduced in the later notebooks.
These de novo ligand is then subjected to ADMET analysis. Valid conformers of the added R-group 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
The history saving thread hit an unexpected error (OperationalError('attempt to write a readonly database')).History will not be written to the database.
<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.
Prepare the ligand scaffold#
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. In this case, we are selecting the hydrogen atom labelled H:8
attachment_index = 8
Add RGroup to your scaffold#
In this tutorial, we show how one can create an R-group from Smiles
R_group_methanol = Chem.AddHs(Chem.MolFromSmiles('*CO'))
R_group_methanol
Build a new molecules#
# we have to specify where the R-group should be attached using the attachment index
rmol = fegrow.build_molecule(scaffold, R_group_methanol, attachment_index)
The R-Group lacks initial coordinates. Defaulting to Chem.rdDistGeom.EmbedMolecule.
[11:23:06] UFFTYPER: Unrecognized atom type: *_ (0)
rmol
| Smiles | Molecule | |
|---|---|---|
| ID | ||
| None | [H]OC([H])([H])c1c([H])nc([H])c([H])c1[H] |  |
rmol.rep2D()
rmol.rep3D()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
<py3Dmol.view at 0x78df756dbb10>
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:23:06] DEPRECATION WARNING: please use MorganGenerator
| MW | HBA | HBD | LogP | Pass_Ro5 | has_pains | has_unwanted_subs | has_prob_fgs | synthetic_accessibility | Smiles | |
|---|---|---|---|---|---|---|---|---|---|---|
| ID | ||||||||||
| None | 109.128 | 2 | 1 | 0.574 | True | False | False | False | 7.299 | [H]OC([H])([H])c1c([H])nc([H])c([H])c1[H] |
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 2 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.03s.
@> 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 0x78df6fc1ffd0>
Any conformers that clash with the protein (any atom-atom distance less than 1 Angstrom), are removed.
rmol.remove_clashing_confs(rec_final)
Removed 0 conformers.
| Smiles | Molecule | |
|---|---|---|
| ID | ||
| None | [H]OC([H])([H])c1c([H])nc([H])c([H])c1[H] | |
rmol.rep3D(prody=rec_final)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
<py3Dmol.view at 0x78df6fdb7810>
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%|███████████████████████| 2/2 [00:02<00:00, 1.11s/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()
2
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 0x78df6a9f9d50>
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"1_mini_rmol_best_conformers.pdb")
print(final_energies)
ID Conformer Energy
0 None 0 0.000
1 None 1 0.628
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.009 | 978768.5363404725 |
| 1 | 0 | 1 | 2.945 | 1133913.693489917 |
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 978768.5363404725
1 1133913.693489917
Name: Kd, dtype: pint[nanomolar][Float64]