import enum
from functools import partial
from typing import Generator
from typing import Union
from rdkit import Chem
from rdkit.Chem import rdChemReactions
from rdkit.Chem import rdmolops
from schrodinger.thirdparty import rdkit_adapter
ATOM_PROP_ATOM_LABEL = "atomLabel"
BOND_PROP_BOND_STEREO = "_MolFileBondStereo"
BOND_PROP_BOND_CONFIG = "_MolFileBondCfg"
MOL_PROP_ATTACHPT = "molAttchpt"
MOL_PROP_R_LABEL = "_MolFileRLabel"
[docs]def convert(data: str, input_format: Format, output_format: Format) -> str:
"""
Main entrypoint for converting one serialization to another. A few
convience functions are provided below for common LD conversions.
:param data: input text string
:param input_format: expected format for input string
:param output_format: desired format for output string
:return: converted text string
"""
if _is_reaction(data, input_format):
return _from_rdkit_reaction(_to_rdkit_reaction(data,
format=input_format),
format=output_format)
else:
return _from_rdkit_mol(_to_rdkit_mol(data, format=input_format),
format=output_format)
[docs]def sdf_to_rdkit(
molblock: str
) -> Union[Chem.rdchem.Mol, Chem.rdChemReactions.ChemicalReaction]:
"""
:param molblock: given SDF or RXN molblock
:return: corresponding RDKit mol or reaction
"""
if _is_reaction(molblock, Format.SDF):
return _to_rdkit_reaction(molblock, format=Format.SDF)
else:
return _to_rdkit_mol(molblock, format=Format.SDF)
[docs]def rdkit_to_sdf(mol: Chem.rdchem.Mol, kekulize: bool = True) -> str:
"""
:param input: given RDKit mol
:param kekulize: whether to kekulize mol
:return: corresponding SDF molblock
"""
return _from_rdkit_mol(mol, format=Format.SDF, kekulize_sdf=kekulize)
[docs]def rdkit_to_cxsmiles(mol: Chem.rdchem.Mol) -> str:
"""
:param mol: given RDKit mol
:return: corresponding CXSMILES stripped of coordinates
"""
return _from_rdkit_mol(mol, format=Format.CXSMILES)
[docs]def get_sd_reader(molblock: str) -> Generator[Chem.rdchem.Mol, None, None]:
"""
:param molblock: string of SDF molblocks
:return: generator that provides mols from SDF molblcoks
"""
reader = Chem.SDMolSupplier()
reader.SetData(molblock,
sanitize=False,
removeHs=False,
strictParsing=False)
for mol in reader:
yield _replace_attachment_points(mol)
[docs]def atom_has_attachment_point(atom):
return atom.HasProp(ATOM_PROP_ATOM_LABEL) and atom.GetProp(
ATOM_PROP_ATOM_LABEL).startswith("_AP")
def _is_reaction(data: str, format: Format) -> bool:
if format == Format.SDF:
return data.startswith("$RXN")
elif format in (Format.SMILES, Format.CXSMILES, Format.SMARTS):
return ">>" in data
else:
return False # Unsupported reaction format
def _to_rdkit_mol(data: str, format: Format) -> Chem.rdchem.Mol:
def _from_sdf(molblock: str) -> Chem.rdchem.Mol:
reader = get_sd_reader(molblock)
mol = next(reader)
if mol is None:
raise RuntimeError("Single molblock required; 0 present")
for m in reader:
# there is another molblock if this can loop
raise RuntimeError("Single molblock required; multiple present")
return mol
_from_smiles = partial(Chem.MolFromSmiles, sanitize=False)
str_to_mol = {
Format.SDF: _from_sdf,
Format.SMILES: _from_smiles,
Format.CXSMILES: _from_smiles,
Format.SMARTS: Chem.MolFromSmarts,
Format.INCHI: partial(Chem.MolFromInchi, sanitize=False),
}
with rdkit_adapter.convert_log_to_exception():
mol = str_to_mol[format](data)
_fix_r0_rgroup(mol)
# partial sanitization to preserve input molecule
ops = (Chem.SANITIZE_ALL ^ Chem.SANITIZE_CLEANUP ^
Chem.SANITIZE_PROPERTIES ^ Chem.SANITIZE_KEKULIZE ^
Chem.SANITIZE_FINDRADICALS ^ Chem.SANITIZE_CLEANUPCHIRALITY)
Chem.SanitizeMol(mol, sanitizeOps=ops)
mol.UpdatePropertyCache(strict=False)
return mol
def _from_rdkit_mol(mol: Chem.rdchem.Mol,
format: Format,
kekulize_sdf: bool = True) -> str:
def _to_cxsmiles(mol: Chem.rdchem.Mol) -> str:
mol.RemoveAllConformers() # remove coordinates from CXSMILES
return Chem.MolToCXSmiles(mol)
def _to_sdf(mol: Chem.rdchem.Mol, kekulize: bool) -> str:
mol = _replace_dummy_atoms_with_attchpts(mol)
# SHARED-8282, SHARED-8498: there's cases that can't be kekulized,
# so just kekulize if possible
if kekulize:
Chem.rdmolops.KekulizeIfPossible(mol, clearAromaticFlags=True)
return Chem.SDWriter.GetText(mol, kekulize=False, force_v3000=True)
mol_to_str = {
Format.SDF: partial(_to_sdf, kekulize=kekulize_sdf),
Format.SMILES: Chem.MolToSmiles,
Format.SMARTS: Chem.MolToSmarts,
Format.CXSMILES: _to_cxsmiles,
Format.INCHI: Chem.MolToInchi,
}
return mol_to_str[format](Chem.Mol(mol))
def _to_rdkit_reaction(data: str,
format: Format) -> Chem.rdChemReactions.ChemicalReaction:
str_to_rxn = {
Format.SDF: rdChemReactions.ReactionFromRxnBlock,
Format.SMILES: partial(rdChemReactions.ReactionFromSmarts,
useSmiles=True),
Format.SMARTS: rdChemReactions.ReactionFromSmarts,
Format.CXSMILES: partial(rdChemReactions.ReactionFromSmarts,
useSmiles=True),
}
# RDKit doesn't support SMARTS extensions in reactions
if format in (Format.SMARTS, Format.SMILES, Format.CXSMILES):
data = data.split(" ")[0]
if format == Format.INCHI:
raise RuntimeError("Cannot convert INCHI format to reaction")
rxn = str_to_rxn[format](data)
for reactant in rxn.GetReactants():
_fix_r0_rgroup(reactant)
for product in rxn.GetProducts():
_fix_r0_rgroup(product)
return rxn
def _from_rdkit_reaction(rxn: Chem.rdChemReactions.ChemicalReaction,
format: Format) -> str:
rxn_to_str = {
Format.SDF: rdChemReactions.ReactionToRxnBlock,
Format.SMILES: rdChemReactions.ReactionToSmiles,
Format.SMARTS: rdChemReactions.ReactionToSmarts,
Format.CXSMILES: rdChemReactions.ReactionToSmiles,
}
if format == Format.INCHI:
raise RuntimeError("Cannot convert reaction to INCHI format")
return rxn_to_str[format](rxn)
def _set_attachment_point_coordinates(mol, edit_mol, attch_pts):
"""
Calculates coordinate of recently added attachment points by using
coordinates of implicit hydrogens in the mol and assigns chiral
tags using bond directions and new coordinates
"""
mol.UpdatePropertyCache(strict=False)
temp_mol = Chem.AddHs(mol, addCoords=True, onlyOnAtoms=attch_pts.keys())
temp_mol_positions = temp_mol.GetConformer().GetPositions()
for parent, attch_pt_list in attch_pts.items():
parent_atom = temp_mol.GetAtomWithIdx(parent)
added_hydrogens = (h for h in parent_atom.GetNeighbors()
if h.GetIdx() >= mol.GetNumAtoms())
for at_idx, h in zip(attch_pt_list, added_hydrogens):
# h is a graph hydrogen whose coordinates we can use for
# the recently added attachment point
xyz = temp_mol_positions[h.GetIdx()]
edit_mol.GetConformer().SetAtomPosition(at_idx, xyz)
mol = reapply_molblock_wedging(edit_mol)
Chem.AssignChiralTypesFromBondDirs(mol)
return Chem.Mol(mol)
def _replace_attachment_points(mol):
"""
Replace attachment points with dummy atoms to preserve chirality
"""
if mol is None:
return None
attch_pts = {
a.GetIdx(): [] for a in mol.GetAtoms() if a.HasProp(MOL_PROP_ATTACHPT)
}
if not attch_pts:
return mol
# Replace attachment points with dummy atoms
edit_mol = Chem.RWMol(mol)
attch_pt_num = 1
for atom_idx in attch_pts:
num_attch_pts = 1
parent_atom = edit_mol.GetAtomWithIdx(atom_idx)
attch_num = parent_atom.GetIntProp(MOL_PROP_ATTACHPT)
if attch_num == -1:
# this atom has two attachment points
num_attch_pts = 2
attch_num = 1
exp_hs = parent_atom.GetNumExplicitHs()
for i in range(num_attch_pts):
new_at = edit_mol.AddAtom(Chem.Atom(0))
attch_pts[atom_idx].append(new_at)
edit_mol.AddBond(atom_idx, new_at, Chem.BondType.SINGLE)
edit_mol.GetAtomWithIdx(new_at).SetProp(ATOM_PROP_ATOM_LABEL,
f'_AP{attch_num}')
attch_num += 1
parent_atom.ClearProp(MOL_PROP_ATTACHPT)
if exp_hs:
parent_atom.SetNumExplicitHs(exp_hs - num_attch_pts)
# Dummy atoms should be query atoms. This ensures the dummy atom is written
# as '*' instead of 'R'
params = rdmolops.AdjustQueryParameters.NoAdjustments()
params.makeDummiesQueries = True
edit_mol = rdmolops.AdjustQueryProperties(edit_mol, params)
return _set_attachment_point_coordinates(mol, edit_mol, attch_pts)
def _replace_dummy_atoms_with_attchpts(mol):
# According to the MDL specification, ATTCHPT=1 is the first attachment point,
# ATTCHPT=2 is the second attachment point, and ATTCHPT=-1 means that this
# atom has the first and second attachment points.
remove_attchpts = []
for at in mol.GetAtoms():
if atom_has_attachment_point(at):
remove_attchpts.append(at.GetIdx())
parent = at.GetNeighbors()[0]
exp_hs = parent.GetNumExplicitHs()
if parent.HasProp(MOL_PROP_ATTACHPT):
parent.SetIntProp(MOL_PROP_ATTACHPT, -1)
parent.SetNumExplicitHs(exp_hs + 2)
else:
parent.SetIntProp(MOL_PROP_ATTACHPT,
int(at.GetProp(ATOM_PROP_ATOM_LABEL)[3]))
parent.SetNumExplicitHs(exp_hs + 1)
if not remove_attchpts:
return mol
# remove attachment points
edit_mol = Chem.RWMol(mol)
edit_mol.BeginBatchEdit()
for at_idx in remove_attchpts:
edit_mol.RemoveAtom(at_idx)
edit_mol.CommitBatchEdit()
return edit_mol.GetMol()
def _detect_unkekulized_atoms(m):
opts = Chem.SANITIZE_ALL ^ Chem.SANITIZE_PROPERTIES ^ Chem.SANITIZE_CLEANUP ^ Chem.SANITIZE_CLEANUPCHIRALITY ^ Chem.SANITIZE_FINDRADICALS
errs = Chem.DetectChemistryProblems(m, opts)
for err in errs:
if err.GetType() == 'KekulizeException':
return sorted(err.GetAtomIndices())
return None
[docs]def reapply_molblock_wedging(mol):
for bnd in mol.GetBonds():
# only change the wedging if the bond was wedged in the input data (we
# recognize that by looking for the properties the mol file parser
# sets):
if bnd.HasProp(BOND_PROP_BOND_STEREO):
# V2000
val = bnd.GetProp(BOND_PROP_BOND_STEREO)
if val == "1":
bnd.SetBondDir(Chem.BondDir.BEGINWEDGE)
elif val == "6":
bnd.SetBondDir(Chem.BondDir.BEGINDASH)
elif val == "4":
bnd.SetBondDir(Chem.BondDir.UNKNOWN)
elif bnd.HasProp(BOND_PROP_BOND_CONFIG):
# V3000
val = bnd.GetProp(BOND_PROP_BOND_CONFIG)
if val == "1":
bnd.SetBondDir(Chem.BondDir.BEGINWEDGE)
elif val == "3":
bnd.SetBondDir(Chem.BondDir.BEGINDASH)
elif val == "2":
bnd.SetBondDir(Chem.BondDir.UNKNOWN)
elif bnd.GetBondDir() in (Chem.BondDir.BEGINDASH,
Chem.BondDir.BEGINWEDGE,
Chem.BondDir.UNKNOWN):
bnd.SetBondDir(Chem.BondDir.NONE)
return mol
[docs]def add_hs_to_aromatic_nitrogen(mol):
"""
Intended to be used with molecules which have kekulization failures due to
aromatic system(s) containing Ns where the user hasn't provided the location of
the implicit Hs in the system.
This picks an arbitrary (but canonical) aromatic N to add an H to in each
aromatic system.
Returns the original mol if it can't fix it.
"""
Chem.GetSymmSSSR(mol)
tm = Chem.Mol(mol)
unkekulized_atoms = _detect_unkekulized_atoms(tm)
while unkekulized_atoms:
tm, new_unkekulized_atoms = _fix_ring_system(tm, unkekulized_atoms)
if new_unkekulized_atoms == unkekulized_atoms:
# got stuck, return the original
return mol
unkekulized_atoms = new_unkekulized_atoms
return tm
def _fix_ring_system(m, unkekulized_atoms):
m.UpdatePropertyCache(False)
# we want to always generate the same result for the same molecule
# Which result to return will cause arguments, yay, but we'll
# prefer smaller rings, then use the canonical atom ranking to break ties
ranks = Chem.CanonicalRankAtoms(m)
ri = m.GetRingInfo()
atoms_to_check = unkekulized_atoms[:]
atoms_to_check.sort(key=lambda x: (ri.MinAtomRingSize(x), ranks[x]))
while atoms_to_check:
ai = atoms_to_check.pop(0)
atom = m.GetAtomWithIdx(ai)
if atom.GetAtomicNum() != 7:
continue
if atom.GetDegree() != 2 or atom.GetNumExplicitHs() != 0:
continue
tm = Chem.Mol(m)
atom = tm.GetAtomWithIdx(ai)
atom.SetNumExplicitHs(1)
atom.SetNoImplicit(True)
# did we fix something?
new_unkekulized_atoms = _detect_unkekulized_atoms(tm)
if new_unkekulized_atoms != unkekulized_atoms:
return tm, new_unkekulized_atoms
return m, unkekulized_atoms
def _fix_r0_rgroup(mol):
highest_rlabel = 0
atoms_with_r0 = []
for atom in mol.GetAtoms():
try:
label = atom.GetIntProp(MOL_PROP_R_LABEL)
except KeyError:
continue
highest_rlabel = max(highest_rlabel, label)
if label == 0:
atoms_with_r0.append(atom)
if atoms_with_r0:
highest_rlabel += 1
for atom in atoms_with_r0:
atom.SetIntProp(MOL_PROP_R_LABEL, highest_rlabel)