🔄 Ring Perception
Problem
You want to determine whether two atoms belong to the same ring. See examples in Table 1.
Ring perception is one of the fundamental algorithms when handling chemical structures. While OEChem TK does not provide a solution for SSSR (Smallest Set of Smallest Rings) it includes a wide range of functions (see Table 2) that answer related questions. This recipe will illustrate how to utilize these functions to determine whether two atoms belong to the same ring.
This recipe will also include examples:
Identifying Spiro Atoms of a Molecule to identify spiro atoms of a molecule
Identifying Macro-cycle Atoms of a Molecule to identify macro-cycle atoms
See also
Smallest Set of Smallest Rings (SSSR) Considered Harmful section in OEChem TK manual
|
|
|
|
API |
See subsection |
|---|---|
Ring Perception section in OEChem TK |
|
Ingredients
|
Difficulty level
🌶️ 🌶️
Download
Source Code
ring_perception
#!/usr/bin/env python3
# (C) 2023 Cadence Design Systems, Inc. (Cadence)
# All rights reserved.
# TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
# provided to current licensees or subscribers of Cadence products or
# SaaS offerings (each a "Customer").
# Customer is hereby permitted to use, copy, and modify the Sample Code,
# subject to these terms. Cadence claims no rights to Customer's
# modifications. Modification of Sample Code is at Customer's sole and
# exclusive risk. Sample Code may require Customer to have a then
# current license or subscription to the applicable Cadence offering.
# THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
# NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
# PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
# liable for any damages or liability in connection with the Sample Code
# or its use.
"""Identify rings in molecule."""
import argparse
import os
import pathlib
import sys
import rich.console
from openeye import oechem
from rich_argparse import HelpPreviewAction, RichHelpFormatter
__SCRIPT_NAME__ = pathlib.Path(__file__).absolute().stem
__SCRIPT_DESC__ = "Find atoms in same ring."
__SCRIPT_TOOLKITS__ = ["oechem"]
__SCRIPT_CATEGORIES__ = ["cheminformatics"]
def parse_options() -> argparse.Namespace:
"""Set up command line options."""
parser = argparse.ArgumentParser(
add_help=True,
formatter_class=RichHelpFormatter,
description="[yellow]" + __SCRIPT_DESC__ + "[/yellow]",
)
# input options
input_group = parser.add_argument_group("Input molecule options")
exclusive_image_group = input_group.add_mutually_exclusive_group(required=True)
exclusive_image_group.add_argument(
"-in",
"--input",
type=str,
required=False,
metavar="MOL-FILE",
help="Input molecule file.",
)
exclusive_image_group.add_argument(
"-s",
"--smiles",
type=str,
metavar="SMILES",
help="Input molecule smiles",
)
parser.add_argument("--help-image", action=HelpPreviewAction)
return parser.parse_args()
def main() -> int:
"""Identify atoms in same ring."""
args = parse_options()
# initialize molecule
mol = oechem.OEGraphMol()
if args.input:
ifs = oechem.oemolistream()
if not ifs.open(args.input):
oechem.OEThrow.Fatal("Cannot open input file!")
if not oechem.OEReadMolecule(ifs, mol):
oechem.OEThrow.Fatal("Cannot read input file!")
elif args.smiles:
if not oechem.OESmilesToMol(mol, args.smiles):
oechem.OEThrow.Fatal("Cannot initialize molecule from smiles!")
if not mol.IsValid():
oechem.OEThrow.Fatal("Input molecule can not be initialized!")
console = rich.console.Console()
for atom in mol.GetAtoms(oechem.OEAtomIsInRing()):
same_ring_atoms = [
a
for a in mol.GetAtoms(oechem.OEAtomIsInRing())
if atoms_in_same_ring(atom, a)
]
if len(same_ring_atoms) != 0:
same_ring_str = " ".join([str(a) for a in same_ring_atoms])
console.print(f"atom {atom!s} in the same ring as {same_ring_str}")
return os.EX_OK
class ChainAtomOrAlreadyTraversed(oechem.OEUnaryAtomPred):
"""Predicate used to find atoms in same ring system."""
def __init__(self, exclude: list[oechem.OEAtomBase]) -> None:
"""Initialize predicate."""
oechem.OEUnaryAtomPred.__init__(self)
self._exclude = exclude
def __call__(self, atom: oechem.OEAtomBase) -> bool:
"""Evaluate atom."""
if not atom.IsInRing():
return False
return atom in self._exclude
def CreateCopy(self): # noqa: ANN201, N802
"""Copy constructor."""
return ChainAtomOrAlreadyTraversed(self._exclude).__disown__()
def atoms_in_same_ring( # noqa: PLR0911
atom_one: oechem.OEAtomBase, atom_two: oechem.OEAtomBase
) -> bool:
"""Determine whether two atoms belong to the same ring system."""
if not atom_one.IsInRing() or not atom_two.IsInRing():
# any of them is a chain atom
return False
if atom_one == atom_two:
return True
if (bond := atom_one.GetBond(atom_two)) is not None and not bond.IsInRing():
return False
first_path = list(
oechem.OEShortestPath(atom_one, atom_two, oechem.OEAtomIsInChain())
)
first_path_length = len(first_path)
if first_path_length == 2: # noqa: PLR2004
return True # neighbors
if first_path_length == 0:
return False # not is same ring system
smallest_one = oechem.OEAtomGetSmallestRingSize(atom_one)
smallest_two = oechem.OEAtomGetSmallestRingSize(atom_two)
if first_path_length > smallest_one and first_path_length > smallest_two:
return False # too far away
# try to find the second shortest different path
exclude_pred = ChainAtomOrAlreadyTraversed(first_path[1:-1])
second_path = list(oechem.OEShortestPath(atom_one, atom_two, exclude_pred))
second_path_length = len(second_path)
if second_path_length == 0:
return False # can not be in the same ring
if second_path_length > smallest_one and second_path_length > smallest_two:
return False # too far away
sum_ring_size = len(first_path) + len(second_path) - 2
if sum_ring_size > smallest_one and sum_ring_size > smallest_two:
return False
in_ring_one = oechem.OEAtomIsInRingSize(atom_one, sum_ring_size)
in_ring_two = oechem.OEAtomIsInRingSize(atom_two, sum_ring_size)
return in_ring_one and in_ring_two
setattr(main, "__SCRIPT_NAME__", __SCRIPT_NAME__)
setattr(main, "__SCRIPT_DESC__", __SCRIPT_DESC__)
setattr(main, "__SCRIPT_TOOLKITS__", __SCRIPT_TOOLKITS__)
setattr(main, "__SCRIPT_CATEGORIES__", __SCRIPT_CATEGORIES__)
if __name__ == "__main__":
sys.exit(main())
Solution
The algorithm implemented in atoms_in_same_ring is based on a simple concept: if two atoms belong to the same ring then there must be at least two alternative ring paths between them. These paths are identified by using the OEShortestPath function.
def atoms_in_same_ring( # noqa: PLR0911
atom_one: oechem.OEAtomBase, atom_two: oechem.OEAtomBase
) -> bool:
"""Determine whether two atoms belong to the same ring system."""
if not atom_one.IsInRing() or not atom_two.IsInRing():
# any of them is a chain atom
return False
if atom_one == atom_two:
return True
if (bond := atom_one.GetBond(atom_two)) is not None and not bond.IsInRing():
return False
first_path = list(
oechem.OEShortestPath(atom_one, atom_two, oechem.OEAtomIsInChain())
)
first_path_length = len(first_path)
if first_path_length == 2: # noqa: PLR2004
return True # neighbors
if first_path_length == 0:
return False # not is same ring system
smallest_one = oechem.OEAtomGetSmallestRingSize(atom_one)
smallest_two = oechem.OEAtomGetSmallestRingSize(atom_two)
if first_path_length > smallest_one and first_path_length > smallest_two:
return False # too far away
# try to find the second shortest different path
exclude_pred = ChainAtomOrAlreadyTraversed(first_path[1:-1])
second_path = list(oechem.OEShortestPath(atom_one, atom_two, exclude_pred))
second_path_length = len(second_path)
if second_path_length == 0:
return False # can not be in the same ring
if second_path_length > smallest_one and second_path_length > smallest_two:
return False # too far away
sum_ring_size = len(first_path) + len(second_path) - 2
if sum_ring_size > smallest_one and sum_ring_size > smallest_two:
return False
in_ring_one = oechem.OEAtomIsInRingSize(atom_one, sum_ring_size)
in_ring_two = oechem.OEAtomIsInRingSize(atom_two, sum_ring_size)
return in_ring_one and in_ring_two
ChainAtomOrAlreadyTraversed is the atom predicate that is used with OEShortestPath to identify a second alternative path between two atoms.
1class ChainAtomOrAlreadyTraversed(oechem.OEUnaryAtomPred):
2 """Predicate used to find atoms in same ring system."""
3
4 def __init__(self, exclude: list[oechem.OEAtomBase]) -> None:
5 """Initialize predicate."""
6 oechem.OEUnaryAtomPred.__init__(self)
7 self._exclude = exclude
8
9 def __call__(self, atom: oechem.OEAtomBase) -> bool:
10 """Evaluate atom."""
11 if not atom.IsInRing():
12 return False
13 return atom in self._exclude
14
15 def CreateCopy(self): # noqa: ANN201, N802
16 """Copy constructor."""
17 return ChainAtomOrAlreadyTraversed(self._exclude).__disown__()
Usage
See Download section to download the script.
> ring_perception --help
> ring_perception --smiles 'c1ccc2c(c1)cc[nH]2'
atom 0 C in the same ring as 0 C 1 C 2 C 3 C 4 C 5 C
atom 1 C in the same ring as 0 C 1 C 2 C 3 C 4 C 5 C
atom 2 C in the same ring as 0 C 1 C 2 C 3 C 4 C 5 C
atom 3 C in the same ring as 0 C 1 C 2 C 3 C 4 C 5 C 6 C 7 C 8 N
atom 4 C in the same ring as 0 C 1 C 2 C 3 C 4 C 5 C 6 C 7 C 8 N
atom 5 C in the same ring as 0 C 1 C 2 C 3 C 4 C 5 C
atom 6 C in the same ring as 3 C 4 C 6 C 7 C 8 N
atom 7 C in the same ring as 3 C 4 C 6 C 7 C 8 N
atom 8 N in the same ring as 3 C 4 C 6 C 7 C 8 N
Discussion
The following examples show how to use ring perception functions available in OEChem TK.
Identifying Spiro Atoms of a Molecule
The atoms_in_same_ring can be used to identify spiro atoms, which are the single, central atom that connects two distinct molecular rings.
class IsSpiroAtom(oechem.OEUnaryAtomPred):
"""
Predicate class to identify spiro atoms.
Spiro atoms are the single, central atom that connects more distinct molecular rings.
"""
def __call__(self, atom: oechem.OEAtomBase) -> bool:
"""Evaluate the atom."""
spiro_atom_hvy_degree = 4
if not atom.IsInRing() or atom.GetHvyDegree() != spiro_atom_hvy_degree:
return False
ring_neighs = [nbr for nbr in atom.GetAtoms() if nbr.IsInRing()]
if len(ring_neighs) != spiro_atom_hvy_degree:
# spiro atom have to have 4 ring neighbors
return False
# spiro atom has to have to 2 pairs of neighbors in 2 different rings)
for neigh in ring_neighs:
num_neighs_in_same_ring = 0
for other_neigh in ring_neighs:
if neigh.GetIdx() == other_neigh.GetIdx():
continue
if atoms_in_same_ring(neigh, other_neigh):
num_neighs_in_same_ring += 1
if num_neighs_in_same_ring != 1:
# bridgehead atom not a spiro atom
return False
return True
scale = oedepict.OEScale_AutoScale
opts = oedepict.OE2DMolDisplayOptions(width, height, scale)
disp = oedepict.OE2DMolDisplay(mol, opts)
highlight = oedepict.OEHighlightByBallAndStick(oechem.OEBlueTint)
oedepict.OEAddHighlighting(disp, highlight, IsSpiroAtom())
oedepict.OERenderMolecule(image, disp)
oedepict.OEWriteImage("depict_spiro_atoms.svg", image)
Download code |
|
Identifying Macro-cycle Atoms of a Molecule
The atoms_in_same_ring can also be used to identify macro-cycle atoms, which are the atoms that are part of large rings in a molecule.
def get_macro_cycle_atoms(
mol: oechem.OEMolBase, min_ring_size: int = 10
) -> list[oechem.OEAtomBase]:
"""Find atoms in rings larger than a specified size."""
oechem.OEFindRingAtomsAndBonds(mol)
if oechem.OECount(mol, oechem.OEAtomIsInRing()) == 0:
return [] # no ring atoms
small_ring_atoms: set[oechem.OEAtomBase] = set()
for atom in mol.GetAtoms(oechem.OEAtomIsInRing()):
if oechem.OEAtomGetSmallestRingSize(atom) <= min_ring_size:
small_ring_atoms.add(atom)
macro_cyclic_ring_seeds: set[oechem.OEAtomBase] = (
set(mol.GetAtoms(oechem.OEAtomIsInRing())) - small_ring_atoms
)
if macro_cyclic_ring_seeds == set():
return [] # no macro-cyclic atoms
macro_cycle_atoms: set[oechem.OEAtomBase] = set()
macro_cycle_atoms.update(macro_cyclic_ring_seeds)
for atom_one in macro_cyclic_ring_seeds:
for atom_two in small_ring_atoms:
if atom_two in macro_cycle_atoms:
continue # already identified as macro-cyclic
if atoms_in_same_ring(atom_one, atom_two):
macro_cycle_atoms.add(atom_two)
return list(macro_cycle_atoms)
macro_cycle_atoms = get_macro_cycle_atoms(mol, min_ring_size=10)
abset = oechem.OEAtomBondSet()
for atom in macro_cycle_atoms:
abset.AddAtom(atom)
scale = oedepict.OEScale_AutoScale
opts = oedepict.OE2DMolDisplayOptions(width, height, scale)
disp = oedepict.OE2DMolDisplay(mol, opts)
highlight = oedepict.OEHighlightByBallAndStick(oechem.OEBlueTint)
oedepict.OEAddHighlighting(disp, highlight, abset)
oedepict.OERenderMolecule(image, disp)
oedepict.OEWriteImage("depict_macro_cycle_atoms.svg", image)
Download code |
|
Identifying Ring Systems of a Molecule
num_rings, ring_list = oechem.OEDetermineRingSystems(mol)
disp = oedepict.OE2DMolDisplay(mol, opts)
highlight = oedepict.OEHighlightByLasso(oechem.OEBlack)
highlight.SetConsiderAtomLabelBoundingBox(True)
ring_pred = oechem.OEPartPredAtom(ring_list)
for ring_idx, color in zip(
range(1, num_rings + 1), oechem.OEGetVividColors(), strict=False
):
ring_pred.SelectPart(ring_idx)
ring_set = oechem.OEAtomBondSet(mol.GetAtoms(ring_pred))
highlight.SetColor(color)
oedepict.OEAddHighlighting(disp, highlight, ring_set)
oedepict.OERenderMolecule(image, disp)
oedepict.OEWriteImage("depict_ring_system.svg", image)
Download code |
|
See also
OEDetermineRingSystems method
Identifying Aromatic Ring Systems of a Molecule
num_rings, ring_list = oechem.OEDetermineAromaticRingSystems(mol)
disp = oedepict.OE2DMolDisplay(mol, opts)
highlight = oedepict.OEHighlightByLasso(oechem.OEBlack)
highlight.SetConsiderAtomLabelBoundingBox(True)
ring_pred = oechem.OEPartPredAtom(ring_list)
for ring_idx, color in zip(
range(1, num_rings + 1), oechem.OEGetVividColors(), strict=False
):
ring_pred.SelectPart(ring_idx)
ring_set = oechem.OEAtomBondSet(mol.GetAtoms(ring_pred))
highlight.SetColor(color)
oedepict.OEAddHighlighting(disp, highlight, ring_set)
oedepict.OERenderMolecule(image, disp)
oedepict.OEWriteImage("depict_aromatic_ring_system.svg", image)
Download code |
|
See also
Identifying Atom in Certain Ring Size
class LabelRingSize(oedepict.OEDisplayAtomPropBase):
"""Functor to label atoms by their ring size."""
def __init__(self, max_ring_size: int) -> None:
"""Initialize predicate."""
self._max_ring_size = max_ring_size
oedepict.OEDisplayAtomPropBase.__init__(self)
def __call__(self, atom: oechem.OEAtomBase) -> str:
"""Return atom label to be displayed."""
if not atom.IsInRing():
return ""
rings = [
r
for r in range(3, self._max_ring_size)
if oechem.OEAtomIsInRingSize(atom, r)
]
if len(rings) == 0:
return ""
return "(" + ",".join([str(r) for r in rings]) + ")"
def CreateCopy(self): # noqa: ANN201, N802
"""Copy constructor."""
return LabelRingSize(self._max_ring_size).__disown__()
scale = oedepict.OEScale_AutoScale
opts = oedepict.OE2DMolDisplayOptions(width, height, scale)
opts.SetAtomColorStyle(oedepict.OEAtomColorStyle_WhiteMonochrome)
opts.SetAtomPropertyFunctor(LabelRingSize(max_ring_size=10))
disp = oedepict.OE2DMolDisplay(mol, opts)
oedepict.OERenderMolecule(image, disp)
oedepict.OEWriteImage("depict_ring_size.svg", image)
Download code |
|
See also
Identifying Atom in Certain Aromatic Ring Size
class LabelAromaticRingSize(oedepict.OEDisplayAtomPropBase):
"""Functor to label atoms by their aromatic ring size."""
def __init__(self, max_ring_size: int) -> None:
"""Initialize predicate."""
self._max_ring_size = max_ring_size
oedepict.OEDisplayAtomPropBase.__init__(self)
def __call__(self, atom: oechem.OEAtomBase) -> str:
"""Return atom label to be displayed."""
if not atom.IsInRing():
return ""
rings = [
r
for r in range(3, self._max_ring_size)
if oechem.OEAtomIsInAromaticRingSize(atom, r)
]
if len(rings) == 0:
return ""
return "(" + ",".join([str(r) for r in rings]) + ")"
def CreateCopy(self): # noqa: ANN201, N802
"""Copy constructor."""
return LabelAromaticRingSize(self._max_ring_size).__disown__()
scale = oedepict.OEScale_AutoScale
opts = oedepict.OE2DMolDisplayOptions(width, height, scale)
opts.SetAtomColorStyle(oedepict.OEAtomColorStyle_WhiteMonochrome)
opts.SetAtomPropertyFunctor(LabelAromaticRingSize(max_ring_size=10))
disp = oedepict.OE2DMolDisplay(mol, opts)
oedepict.OERenderMolecule(image, disp)
oedepict.OEWriteImage("depict_aromatic_ring_size.svg", image)
Download code |
|
Identifying Atoms’ Smallest Ring Size
"""Functor to label atoms by their smallest ring size."""
def __init__(self) -> None:
"""Initialize predicate."""
oedepict.OEDisplayAtomPropBase.__init__(self)
def __call__(self, atom: oechem.OEAtomBase) -> str:
"""Return atom label to be displayed."""
if not atom.IsInRing():
return ""
smallest = oechem.OEAtomGetSmallestRingSize(atom)
return "(" + str(smallest) + ")"
def CreateCopy(self): # noqa: ANN201, N802
"""Copy constructor."""
return LabelSmallestRingSize().__disown__()
scale = oedepict.OEScale_AutoScale
opts = oedepict.OE2DMolDisplayOptions(width, height, scale)
opts.SetAtomColorStyle(oedepict.OEAtomColorStyle_WhiteMonochrome)
opts.SetAtomPropertyFunctor(LabelSmallestRingSize())
disp = oedepict.OE2DMolDisplay(mol, opts)
oedepict.OERenderMolecule(image, disp)
oedepict.OEWriteImage("depict_smallest_ring_size.svg", image)
Download code |
|
See also in OEChem TK manual
Theory
API
OEAtomBondSet class
OEPartPredAtom predicate