Source code for schrodinger.livedesign.setup_reaction
"""
Users commonly sketch reactions that are not directly useable in
reaction enumeration. This module includes tools that are needed to convert
a reaction to a suitable format for use in reaction enumeration.
Here's how LiveDesign calls the reaction enumeration tool:
seurat/TaskEngineTasks/python/schrodinger_enumerate/schrodinger_enumerate.py
This requires jaguar to be available to work.
"""
from rdkit import Chem
from rdkit.Chem import rdmolops
from rdkit.Chem import rdqueries
from rdkit.Chem.rdChemReactions import ChemicalReaction
from schrodinger.livedesign import convert
from schrodinger.livedesign import preprocessor
from schrodinger.thirdparty import rdkit_adapter
from schrodinger.structure._structure import NO_SPECIFIC_ISOTOPE
try:
from schrodinger.application.jaguar.packages import reaction_mapping
except ImportError as e:
# in build without jaguar
reaction_mapping = None
RDK_LABEL_PROP = 'i_tmp_rdkmol'
ST_MAP_NUMBER_PROP = 'i_rdkit_molAtomMapNumber'
def _setup_rgroup_rxn_map_nums(mol):
"""
Turn indexed R groups into dummy atoms with atom map numbers > 100
Reactions are commonly drawn with mapped R groups, and nothing else mapped. The R
groups should be turned into star (wildcard) atoms with the correct atom mapping.
to be tested:
- collisions if some existing atoms are mapped
"""
for a in mol.GetAtoms():
try:
rlabel = a.GetUnsignedProp(preprocessor.MOL_PROP_R_LABEL)
except KeyError:
continue
a.SetAtomMapNum(rlabel + 100)
a.ClearProp(preprocessor.MOL_PROP_R_LABEL)
a.ClearProp(preprocessor.ATOM_PROP_DUMMY_LABEL)
def _lock_attachments(mol: Chem.Mol) -> Chem.Mol:
"""
If there are R groups, include all hydrogens in the query so that connections can
only be formed at the R group.
for example:
CCCR -> [CH3][CH2][CH2]R
"""
mol = Chem.Mol(mol)
Chem.SanitizeMol(mol, sanitizeOps=Chem.SANITIZE_ADJUSTHS)
mol.UpdatePropertyCache(strict=False)
mol_with_hydrogens = Chem.AddHs(mol)
return rdmolops.MergeQueryHs(mol_with_hydrogens)
def _has_rgroup_or_star(mol: Chem.Mol) -> bool:
for a in mol.GetAtoms():
if a.HasProp(preprocessor.MOL_PROP_R_LABEL):
return True
if a.DescribeQuery().strip() == 'AtomNull':
return True
return False
def _get_rdk_atom(st_atom, rdk_mols):
rdk_mol_idx = st_atom.structure.property[RDK_LABEL_PROP]
rdk_atom_idx = st_atom.property[rdkit_adapter.adapter.RDK_INDEX]
return rdk_mols[rdk_mol_idx].GetAtomWithIdx(rdk_atom_idx)
def _apply_rxn_mapping(reactants, products, rxn_mapping):
"""
Add the map indices to the atoms in the reaction.
:return: list of pairs of atoms that correspond
"""
map_num = 0
reacting_atoms = []
for reactant_atom, product_atom in rxn_mapping.items():
if reactant_atom is None or product_atom is None:
# atom appears/disappears in the reaction
continue
try:
reactant_rdk_atom = _get_rdk_atom(reactant_atom, reactants)
product_rdk_atom = _get_rdk_atom(product_atom, products)
except KeyError:
# Atoms can be added/removed in the reaction, and won't be mapped.
# Just skip them.
continue
reactant_map_number = reactant_rdk_atom.GetAtomMapNum()
if reactant_map_number >= 100:
if reactant_map_number == product_rdk_atom.GetAtomMapNum():
# Leave R Groups alone (we already mapped them).
continue
elif product_rdk_atom.GetAtomMapNum() < 100:
raise RuntimeError('R group remapped to non- R group atom')
else:
raise RuntimeError('R group remapped different R group')
elif product_rdk_atom.GetAtomMapNum() >= 100:
raise RuntimeError('R group remapped to non- R group atom')
map_num += 1
reactant_rdk_atom.SetAtomMapNum(map_num)
product_rdk_atom.SetAtomMapNum(map_num)
reacting_atoms.append((reactant_rdk_atom, product_rdk_atom))
return reacting_atoms
def _has_query_charge(atom):
"""Does the atom have a charge or a charge query?"""
return bool(atom.GetFormalCharge() or 'Charge' in atom.DescribeQuery())
def _replace_with_0charge(atom):
"""Replace an "implicit" charge 0 with a query for exactly 0 charge."""
if atom.HasQuery():
atom.ExpandQuery(
rdqueries.FormalChargeEqualsQueryAtom(atom.GetFormalCharge()))
else:
qa = rdqueries.AtomNumEqualsQueryAtom(atom.GetAtomicNum())
qa.ExpandQuery(
rdqueries.FormalChargeEqualsQueryAtom(atom.GetFormalCharge()))
mol = atom.GetOwningMol()
idx = atom.GetIdx()
mol.ReplaceAtom(idx, qa)
def _fix_charge_modifing_atoms(reacting_atoms: list):
"""
Change [N+:1]>>[N:1] to [N+:1]>>[N+0:1]
If charge isn't specified in a reaction product, then the charge is
propagated from the reactant - so the input reaction would have done
nothing! Our sketcher doesn't allow drawing an "explicit 0 charge",
so we need to guess if that's what the user meant.
"""
for reactant_atom, product_atom in reacting_atoms:
if _has_query_charge(reactant_atom) ^ _has_query_charge(product_atom):
if not _has_query_charge(reactant_atom):
_replace_with_0charge(reactant_atom)
if not _has_query_charge(product_atom):
_replace_with_0charge(product_atom)
def _is_removable_hydrogen(atom):
return atom.atomic_number == 1 and atom.isotope == NO_SPECIFIC_ISOTOPE and ST_MAP_NUMBER_PROP not in atom.property
def _remove_unnecessary_hs(st):
"""
Removes Hydrogens that are not necessary in reaction mapping to prevent
unnecessary computation when mapping reaction with jaguar.
:param: The structure to be mapped in a reaction.
:return: The input Structure with all removable hydrogens deleted from it.
"""
remove_indices = [at.index for at in st.atom if _is_removable_hydrogen(at)]
st.deleteAtoms(remove_indices)
return st
[docs]def map_reaction(reactants, products):
"""
Call Jaguar's reaction mapper.
:param reactants: list of schrodinger Structures
:param products: list of schrodinger Structures
:return: list of pairs of atoms that correspond
"""
# sets this property to tell the mapper which atoms match.
for a in (a for st in (reactants + products) for a in st.atom):
try:
map_number = a.property[ST_MAP_NUMBER_PROP]
a.property[reaction_mapping.AUTOTS_ATOM_CLASS] = map_number
except KeyError:
pass
return reaction_mapping.get_rp_map(reactants, products)
[docs]def setup_reaction(rxn: ChemicalReaction) -> ChemicalReaction:
"""
Converts display reaction into a reaction useable by core suite's enumeration tool.
Users commonly sketch reactions without atom mappings, and that assume
certain things.
Reaction enumeration requires that:
- In reactants, hydrogens attached to aromatic nitrogens must be explicit
- Reactants be aromatized rather than kekulized
- Each R-group must be converted to '*' (any atom) with a mapping number >100
- If there are R groups, that's the only place the fragments should be allowed to connect
- Atoms must be mapped from the reactant to the product.
"""
params = Chem.AdjustQueryParameters.NoAdjustments()
params.aromatizeIfPossible = True
# May want some of these, let's see what users draw:
# params.adjustSingleBondsToDegreeOneNeighbors = True
# params.adjustSingleBondsBetweenAromaticAtoms = True
has_rlabel = []
reactants = []
st_reactants = []
for i, m in enumerate(rxn.GetReactants()):
has_rlabel.append(_has_rgroup_or_star(m))
m = convert.add_hs_to_aromatic_nitrogen(m)
m = Chem.AdjustQueryProperties(m, params)
_setup_rgroup_rxn_map_nums(m)
m.UpdatePropertyCache(False)
Chem.rdCoordGen.AddCoords(m)
reactants.append(m)
# Jaguar reaction mapper requires coordinates to be present
st = rdkit_adapter.from_rdkit(m, conformer=-1)
st.property[RDK_LABEL_PROP] = i
st = _remove_unnecessary_hs(st)
st_reactants.append(st)
products = []
st_products = []
for i, m in enumerate(rxn.GetProducts()):
_setup_rgroup_rxn_map_nums(m)
m.UpdatePropertyCache(False)
Chem.rdCoordGen.AddCoords(m)
products.append(m)
# Jaguar reaction mapper requires coordinates to be present
st = rdkit_adapter.from_rdkit(m, conformer=-1)
st.property[RDK_LABEL_PROP] = i
st = _remove_unnecessary_hs(st)
st_products.append(st)
# This happens at the end in case atoms are added/removed from the reaction during processing
rxn_mapping = map_reaction(st_reactants, st_products)
if rxn_mapping is None:
raise RuntimeError(
'Reaction atoms could not be mapped between reactants and products')
# Always returns only one mapping
reacting_atoms = _apply_rxn_mapping(reactants, products, rxn_mapping[0])
_fix_charge_modifing_atoms(reacting_atoms)
rxn_mapped = ChemicalReaction()
for m in reactants:
if has_rlabel[i]:
m = _lock_attachments(m)
rxn_mapped.AddReactantTemplate(m)
for m in products:
rxn_mapped.AddProductTemplate(m)
return rxn_mapped