diff --git a/chebi_utils/extract_properties.py b/chebi_utils/extract_properties.py index c96f351..a012219 100644 --- a/chebi_utils/extract_properties.py +++ b/chebi_utils/extract_properties.py @@ -238,12 +238,281 @@ def get_rings(mol: Chem.Mol) -> dict[str, list]: return atom_extensions +def _add_angular_methyls( + mol: Chem.Mol, + atom_extensions: dict[str, list], + iupac_to_atom: dict[int, int], +) -> None: + """Label angular methyls C18 (on C13) and C19 (on C10) when present. + + Looks for the unique carbon neighbour of the attachment atom that is not + part of the gonane core (C1–C17). Estrogens without C19 simply yield no + ``steroid_19`` predicate. + """ + core_atom_indices = set(iupac_to_atom.values()) + for attachment_position, methyl_position in ((13, 18), (10, 19)): + attachment_atom_idx = iupac_to_atom.get(attachment_position) + if attachment_atom_idx is None: + continue + methyl_candidates = [ + neighbor.GetIdx() + for neighbor in mol.GetAtomWithIdx(attachment_atom_idx).GetNeighbors() + if neighbor.GetAtomicNum() == 6 and neighbor.GetIdx() not in core_atom_indices + ] + if len(methyl_candidates) == 1: + atom_extensions.setdefault(f"steroid_{methyl_position}", []).append( + methyl_candidates[0] + ) + + + +# --- NEW: C17 side chain numbering (C20, C21, …) --- + + +def _carbon_neighbor_indices(mol: Chem.Mol, atom_idx: int, excluded: set[int]) -> list[int]: + return [ + neighbor.GetIdx() + for neighbor in mol.GetAtomWithIdx(atom_idx).GetNeighbors() + if neighbor.GetAtomicNum() == 6 and neighbor.GetIdx() not in excluded + ] + + +def _reachable_carbon_indices(mol: Chem.Mol, start_idx: int, blocked: set[int]) -> set[int]: + """Carbon atoms reachable from ``start_idx`` without entering ``blocked``.""" + reachable = set() + stack = [start_idx] + while stack: + atom_idx = stack.pop() + if atom_idx in reachable or atom_idx in blocked: + continue + reachable.add(atom_idx) + stack.extend(_carbon_neighbor_indices(mol, atom_idx, blocked | reachable)) + return reachable + + +def _longest_carbon_path( + mol: Chem.Mol, start_idx: int, allowed: set[int] +) -> list[int]: + """Longest simple carbon path starting at ``start_idx`` within ``allowed``.""" + + def dfs(atom_idx: int, visited: set[int]) -> list[int]: + best_tail: list[int] = [] + for neighbor_idx in _carbon_neighbor_indices(mol, atom_idx, visited): + if neighbor_idx not in allowed: + continue + tail = dfs(neighbor_idx, visited | {neighbor_idx}) + if len(tail) > len(best_tail): + best_tail = tail + return [atom_idx] + best_tail + + return dfs(start_idx, {start_idx}) + + +def _label_steroid_atom( + atom_extensions: dict[str, list], + iupac_to_atom: dict[int, int], + position: int, + atom_idx: int, +) -> None: + atom_extensions.setdefault(f"steroid_{position}", []).append(atom_idx) + iupac_to_atom[position] = atom_idx + + +def _add_side_chain( + mol: Chem.Mol, + atom_extensions: dict[str, list], + iupac_to_atom: dict[int, int], +) -> None: + """Label the C17 side chain with IUPAC numbers C20, C21, … + + C20 is the unique non-core carbon neighbour of C17. From there the longest + carbon path is treated as the main chain; a methyl branch on C20 (if any) + is C21, then the main chain continues as C22, C23, … Terminal methyls on + the last main-chain carbon (e.g. C26/C27 of cholestane) get the next + numbers. Shorter chains (pregnane: C20–C21 only) and missing chains + (estrone) are handled naturally. + """ + c17_atom_idx = iupac_to_atom.get(17) + if c17_atom_idx is None: + return + + # Core C1–C17 plus any already labeled atoms (C18/C19, …). + labeled_atom_indices = set(iupac_to_atom.values()) + for atoms in atom_extensions.values(): + labeled_atom_indices.update(atoms) + + c20_candidates = _carbon_neighbor_indices(mol, c17_atom_idx, labeled_atom_indices) + if len(c20_candidates) != 1: + return + + c20_atom_idx = c20_candidates[0] + _label_steroid_atom(atom_extensions, iupac_to_atom, 20, c20_atom_idx) + + side_chain_atoms = _reachable_carbon_indices(mol, c20_atom_idx, labeled_atom_indices) + side_chain_beyond_c20 = side_chain_atoms - {c20_atom_idx} + next_position = 21 + c20_side_neighbors = _carbon_neighbor_indices( + mol, c20_atom_idx, labeled_atom_indices | {c17_atom_idx} + ) + + if not c20_side_neighbors: + return + + if len(c20_side_neighbors) == 1: + main_chain = _longest_carbon_path( + mol, c20_side_neighbors[0], side_chain_beyond_c20 + ) + else: + neighbor_paths = { + neighbor_idx: _longest_carbon_path(mol, neighbor_idx, side_chain_beyond_c20) + for neighbor_idx in c20_side_neighbors + } + main_chain_start = max( + neighbor_paths, key=lambda neighbor_idx: len(neighbor_paths[neighbor_idx]) + ) + main_chain = neighbor_paths[main_chain_start] + for branch_idx in sorted( + neighbor_idx + for neighbor_idx in c20_side_neighbors + if neighbor_idx != main_chain_start + ): + _label_steroid_atom(atom_extensions, iupac_to_atom, next_position, branch_idx) + next_position += 1 + + # Cholestane terminus: longest path may absorb one of C26/C27; stop at C25. + while len(main_chain) >= 2: + last_atom_idx = main_chain[-1] + previous_atom_idx = main_chain[-2] + leaf_methyls_on_previous = [ + neighbor_idx + for neighbor_idx in _carbon_neighbor_indices( + mol, previous_atom_idx, {c20_atom_idx, c17_atom_idx} + ) + if neighbor_idx in side_chain_atoms + and neighbor_idx not in main_chain[:-1] + and len(_carbon_neighbor_indices(mol, neighbor_idx, set())) == 1 + ] + if last_atom_idx in leaf_methyls_on_previous and len(leaf_methyls_on_previous) >= 2: + main_chain = main_chain[:-1] + else: + break + + for atom_idx in main_chain: + _label_steroid_atom(atom_extensions, iupac_to_atom, next_position, atom_idx) + next_position += 1 + + if not main_chain: + return + + labeled_now = set(iupac_to_atom.values()) + for atoms in atom_extensions.values(): + labeled_now.update(atoms) + terminal_methyls = sorted( + neighbor_idx + for neighbor_idx in _carbon_neighbor_indices(mol, main_chain[-1], labeled_now) + if neighbor_idx in side_chain_atoms + ) + for atom_idx in terminal_methyls: + _label_steroid_atom(atom_extensions, iupac_to_atom, next_position, atom_idx) + next_position += 1 + + +def _c5_steroid_parity(mol: Chem.Mol, iupac_to_atom: dict[int, int]) -> int | None: + """Local C5 parity relative to steroid numbering C4, C6, C10. + + Remaps RDKit CW/CCW to a fixed neighbour order. Returns ``+1`` / ``-1``, + or ``None`` if C5 is missing, not tetrahedral, or has unexpected neighbours. + On the current labelled set: ``+1`` ↔ 5β, ``-1`` ↔ 5α. + """ + try: + c4_idx = iupac_to_atom[4] + c5_idx = iupac_to_atom[5] + c6_idx = iupac_to_atom[6] + c10_idx = iupac_to_atom[10] + except KeyError: + return None + + c5_atom = mol.GetAtomWithIdx(c5_idx) + chiral_tag = c5_atom.GetChiralTag() + if chiral_tag == Chem.ChiralType.CHI_UNSPECIFIED: + return None + + neighbor_indices = [neighbor.GetIdx() for neighbor in c5_atom.GetNeighbors()] + heavy_neighbor_indices = [ + atom_idx + for atom_idx in neighbor_indices + if mol.GetAtomWithIdx(atom_idx).GetAtomicNum() != 1 + ] + hydrogen_neighbor_indices = [ + atom_idx + for atom_idx in neighbor_indices + if mol.GetAtomWithIdx(atom_idx).GetAtomicNum() == 1 + ] + wanted_heavy_order = [c4_idx, c6_idx, c10_idx] + if set(heavy_neighbor_indices) != set(wanted_heavy_order): + return None + if len(heavy_neighbor_indices) != 3: + return None + + # RDKit CW/CCW follows GetNeighbors() order. With implicit H that is the + # three heavy atoms; with explicit H the H is included in the order. + if not hydrogen_neighbor_indices: + observed_order = heavy_neighbor_indices + wanted_order = wanted_heavy_order + elif len(hydrogen_neighbor_indices) == 1: + observed_order = neighbor_indices + wanted_order = wanted_heavy_order + hydrogen_neighbor_indices + else: + return None + + permutation = [wanted_order.index(atom_idx) for atom_idx in observed_order] + inversions = sum( + 1 + for left in range(len(permutation)) + for right in range(left + 1, len(permutation)) + if permutation[left] > permutation[right] + ) + permutation_sign = -1 if inversions % 2 else 1 + + if chiral_tag == Chem.ChiralType.CHI_TETRAHEDRAL_CW: + rdkit_sign = 1 + elif chiral_tag == Chem.ChiralType.CHI_TETRAHEDRAL_CCW: + rdkit_sign = -1 + else: + return None + + return rdkit_sign * permutation_sign + + +def _add_c5_ab_configuration( + mol: Chem.Mol, + atom_extensions: dict[str, list], + iupac_to_atom: dict[int, int], +) -> None: + """Label C5 as ``steroid_5_alpha`` or ``steroid_5_beta`` when tetrahedral. + + Mapping from local parity (C4/C6/C10 order): ``+1`` → beta, ``-1`` → alpha. + """ + c5_atom_idx = iupac_to_atom.get(5) + if c5_atom_idx is None: + return + parity = _c5_steroid_parity(mol, iupac_to_atom) + if parity == 1: + atom_extensions.setdefault("steroid_5_beta", []).append(c5_atom_idx) + elif parity == -1: + atom_extensions.setdefault("steroid_5_alpha", []).append(c5_atom_idx) + + def get_steroid_positions(mol: Chem.Mol) -> dict[str, list]: """Extract steroid-nucleus position predicates. Matches the molecule against the gonane core and, on a match, labels the ring atoms with their IUPAC steroid position as predicates ``steroid_1`` … - ``steroid_17``. Molecules without a gonane core yield no predicates. + ``steroid_17``. When present, angular methyls are added as ``steroid_18`` + (on C13) and ``steroid_19`` (on C10), the C17 side chain as + ``steroid_20``, ``steroid_21``, … (cholestane-style), and tetrahedral C5 + configuration as ``steroid_5_alpha`` / ``steroid_5_beta``. Molecules without + a gonane core yield no predicates. Parameters ---------- @@ -253,16 +522,25 @@ def get_steroid_positions(mol: Chem.Mol) -> dict[str, list]: Returns ------- dict[str, list] - Mapping from ``steroid_{position}`` to the ``list[int]`` of matched atom - indices. Empty when the molecule has no steroid nucleus. + Mapping from ``steroid_{position}`` (and optional ``steroid_5_alpha`` / + ``steroid_5_beta``) to the ``list[int]`` of matched atom indices. Empty + when the molecule has no steroid nucleus. """ atom_extensions: dict[str, list] = {} steroid_match = mol.GetSubstructMatch(_GONANE_PATTERN, useChirality=False) - if steroid_match: - for pat_idx, atom_idx in enumerate(steroid_match): - iupac = _GONANE_IDX_TO_IUPAC.get(pat_idx) - if iupac is not None: - atom_extensions.setdefault(f"steroid_{iupac}", []).append(atom_idx) + if not steroid_match: + return atom_extensions + + iupac_to_atom: dict[int, int] = {} + for pat_idx, atom_idx in enumerate(steroid_match): + iupac = _GONANE_IDX_TO_IUPAC.get(pat_idx) + if iupac is not None: + atom_extensions.setdefault(f"steroid_{iupac}", []).append(atom_idx) + iupac_to_atom[iupac] = atom_idx + + _add_angular_methyls(mol, atom_extensions, iupac_to_atom) + _add_side_chain(mol, atom_extensions, iupac_to_atom) # NEW: C20+ + _add_c5_ab_configuration(mol, atom_extensions, iupac_to_atom) return atom_extensions @@ -289,3 +567,38 @@ def get_numerical_facts(mol: Chem.Mol) -> dict[str, list]: for ring in mol.GetRingInfo().AtomRings(): atom_extensions.setdefault("ring_size", []).append(len(ring)) return atom_extensions + + +"""Manual check for steroid numbering (not part of the library API). + +Run: python -m chebi_utils.extract_properties +Expect: cholesterol -> steroid_1..27 (no 5α/5β: Δ5), estrone -> steroid_1..18, +pregnenolone -> steroid_1..21 (with 18/19), 5β-androstane -> … + steroid_5_beta, +5α-pregnene -> … + steroid_5_alpha, benzene -> [] +""" +if __name__ == "__main__": + from chebi_utils.read_molecule import smiles_or_inchi_to_mol + + for name, smiles in { + "cholesterol": ( + "C[C@H](CCCC(C)C)[C@H]1CC[C@@H]2[C@@]1(CC[C@H]3[C@H]2CC=C4[C@@]3(CC[C@@H](C4)O)C)C" + ), + "estrone": "C[C@]12CC[C@H]3[C@H]([C@@H]1CCC2=O)CCc4c3ccc(O)c4", + "pregnenolone": ( + "CC(=O)[C@H]1CC[C@@H]2[C@@]1(CC[C@H]3[C@H]2CC=C4[C@@]3(CC[C@@H](C4)O)C)C" + ), + "5beta-androstane-3beta,17alpha-diol": ( + "[H][C@]12CC[C@]3([H])[C@]([H])(CC[C@]4(C)[C@H](O)CC[C@@]34[H])" + "[C@@]1(C)CC[C@H](O)C2" + ), + "5alpha-Pregn-2-en-20-one": ( + "[H][C@]12CC=CC[C@]1(C)[C@@]1([H])CC[C@]3(C)[C@@H](C(C)=O)" + "CC[C@@]3([H])[C@]1([H])CC2" + ), + "benzene": "c1ccccc1", + }.items(): + keys = sorted( + get_steroid_positions(smiles_or_inchi_to_mol(smiles)), + key=lambda k: (int(k.split("_")[1]), k), + ) + print(name, keys) diff --git a/chebi_utils/steroid_cip_ab.py b/chebi_utils/steroid_cip_ab.py new file mode 100644 index 0000000..2cefbf2 --- /dev/null +++ b/chebi_utils/steroid_cip_ab.py @@ -0,0 +1,56 @@ +"""Demo: C5 alpha/beta from get_steroid_positions. + +Run: python -m chebi_utils.steroid_cip_ab +""" + +from chebi_utils.extract_properties import get_steroid_positions +from chebi_utils.read_molecule import smiles_or_inchi_to_mol + + +def _c5_ab_label(positions: dict[str, list]) -> str: + if "steroid_5_alpha" in positions: + return "alpha" + if "steroid_5_beta" in positions: + return "beta" + return "-" + + +# Paste substances here: (name, SMILES_or_InChI, "alpha"|"beta") +EXAMPLES: list[tuple[str, str, str]] = [ + ("3-oxo-5β-steroid", "*C1CCC2C3CC[C@]4([H])CC(=O)CCC4(C)C3CCC12C", "beta"), + ("5beta-androstane-3beta,17alpha-diol", "[H][C@]12CC[C@]3([H])[C@]([H])(CC[C@]4(C)[C@H](O)CC[C@@]34[H])[C@@]1(C)CC[C@H](O)C2", "beta"), + ("5beta-tetrahydro-11-deoxycorticosterone", "[H][C@]12CC[C@]3([H])[C@]([H])(CC[C@@]4(C)[C@@]3([H])CC[C@]4([H])C(=O)CO)[C@@]1(C)CC[C@@H](O)C2", "beta"), + ("5β-chola-3,8(14),11-trien-24-oic acid", "[H][C@@]12C=CCC[C@]1(C)[C@@]1([H])C=C[C@@]3(C)C(=C1CC2)CC[C@]3([H])[C@H](C)CCC(=O)O", "beta"), + ("5β-chola-3,8(14),11-trien-24-oic acid", "[H][C@@]12C=CCC[C@]1(C)[C@@]1([H])C=C[C@@]3(C)C(=C1CC2)CC[C@]3([H])[C@H](C)CCC(=O)O", "beta"), + ("andrastin A", "[H][C@@]12C=C(C)[C@@]3(C)C(=O)C(C)=C(O)[C@@]3(C(=O)OC)[C@@]1(C)CC[C@]1([H])C(C)(C)[C@@H](OC(C)=O)CC[C@@]21C=O", "beta"), + ("andrastin B", "[H][C@@]12C=C(C)[C@@]3(C)C(=O)C(C)=C(O)[C@@]3(C(=O)OC)[C@@]1(C)CC[C@]1([H])C(C)(C)[C@@H](OC(C)=O)CC[C@@]21CO", "beta"), + ("andrastin C", "[H][C@@]12C=C(C)[C@@]3(C)C(=O)C(C)=C(O)[C@@]3(C(=O)OC)[C@@]1(C)CC[C@]1([H])C(C)(C)[C@@H](OC(C)=O)CC[C@@]21C", "beta"), + ("andrastin D", "[H][C@@]12C=C(C)[C@@]3(C)C(=O)C(C)=C(O)[C@@]3(C(=O)OC)[C@@]1(C)CC[C@]1([H])C(C)(C)C(=O)CC[C@@]21C", "beta"), + ("andrastin E", "[H][C@@]12C=C(C)[C@@]3(C)C(=O)C(C)=C(O)[C@@]3(C(=O)OC)[C@@]1(C)CC[C@]1([H])C(C)(C)[C@H](O)CC[C@@]21C", "beta"), + ("andrastin F", "[H][C@@]12C=C(C)[C@@]3(C)C(=O)C(C)=C(O)[C@@]3(C(=O)OC)[C@@]1(C)CC[C@]1([H])C(C)(C)[C@@H](O)CC[C@@]21C", "beta"), + ("Tetrahydro-11-dehydrocorticosterone", "[H][C@]12CC[C@@]3([H])[C@]4([H])CC[C@]([H])(C(=O)CO)[C@@]4(C)CC(=O)[C@]3([H])[C@@]1(C)CC[C@@H](O)C2", "beta"), + ("25S)-5alpha-cholestan..", "[H][C@@]12[C@@H](O)[C@@H](O)[C@]3(O)[C@]([H])(CC[C@]4(C)[C@@]([H])([C@H](C)CCC[C@H](C)CO)[C@@H](O)[C@@H](O)[C@]43[H])[C@@]1(C)CC[C@H](O)[C@@H]2O", "alpha"), + ("(25S)-5alpha-cholestan-3beta ...", "[H][C@@]12[C@@H](O)C[C@]3(O)[C@]([H])(CC[C@]4(C)[C@@]([H])([C@H](C)CCC[C@H](C)CO)[C@@H](O)[C@H](O)[C@]43[H])[C@@]1(C)CC[C@H](O)[C@@H]2O", "alpha"), + ("(25S)-5alpha-cholestan-3beta,6alpha,7beta,8", "[H][C@]1([C@H](C)CCC[C@H](C)CO)[C@@H](O)[C@H](O)[C@]2([H])[C@]1(C)CC[C@]1([H])[C@@]3(C)CC[C@H](O)C[C@]3([H])[C@@H](O)[C@H](O)[C@]12O", "alpha"), + ("(25S)-5alpha-cholestan-3beta,6al...", "[H][C@]1([C@H](C)CCC[C@H](C)CO)[C@@H](O)[C@@H](O)[C@]2([H])[C@]1(C)CC[C@]1([H])[C@@]3(C)CC[C@H](O)C[C@]3([H])[C@@H](O)C[C@]12O", "alpha"), + ("(25S)-5alpha-cholestan-3beta,6bet..", "[H][C@]1([C@H](C)CCC[C@H](C)CO)[C@@H](O)[C@H](O)[C@@]2([H])[C@]3([H])C[C@@H](O)[C@@]4([H])C[C@@H](O)CC[C@]4(C)[C@@]3([H])CC[C@]12C", "alpha"), + ("5alpha-Dermophol", "[H][C@]12C[C@H](O)CC[C@]1(C)[C@@]1([H])C[C@H](O)[C@@]3(C)[C@@]([H])(CC[C@]3([H])[C@H](C)CCCC(O)(CO)CO)[C@]1([H])[C@H](O)C2", "alpha"), + ("5alpha-Pregn-2-en-20-one", "[H][C@]12CC=CC[C@]1(C)[C@@]1([H])CC[C@]3(C)[C@@H](C(C)=O)CC[C@@]3([H])[C@]1([H])CC2", "alpha"), +] + + +if __name__ == "__main__": + rows: list[tuple[str, str, str]] = [] + for name, smiles_or_inchi, expected in EXAMPLES: + mol = smiles_or_inchi_to_mol(smiles_or_inchi) + if mol is None: + rows.append((name, expected, "parse_failed")) + continue + rows.append((name, expected, _c5_ab_label(get_steroid_positions(mol)))) + + name_width = max(len(name) for name, _, _ in rows) + col_width = max(len("expected"), max(len(value) for _, expected, predicted in rows for value in (expected, predicted))) + print(f"{'molecule':<{name_width}} | {'expected':<{col_width}} | predicted") + print(f"{'-' * name_width}--+--{'-' * col_width}--+----------") + for name, expected, predicted in rows: + print(f"{name:<{name_width}} | {expected:<{col_width}} | {predicted}") diff --git a/tests/test_extract_properties.py b/tests/test_extract_properties.py new file mode 100644 index 0000000..d5dac22 --- /dev/null +++ b/tests/test_extract_properties.py @@ -0,0 +1,107 @@ +"""Tests for chebi_utils.extract_properties steroid numbering.""" + +from __future__ import annotations + +import pytest + +from chebi_utils.extract_properties import get_steroid_positions, mol_to_fol_atoms +from chebi_utils.read_molecule import smiles_or_inchi_to_mol + +# Cholesterol: gonane core + C18/C19 + cholestane side chain C20–C27 +CHOLESTEROL_SMILES = ( + "C[C@H](CCCC(C)C)[C@H]1CC[C@@H]2[C@@]1(CC[C@H]3[C@H]2CC=C4[C@@]3(CC[C@@H](C4)O)C)C" +) +# Estrone: aromatic ring A, C18 present, C19 and side chain absent +ESTRONE_SMILES = "C[C@]12CC[C@H]3[C@H]([C@@H]1CCC2=O)CCc4c3ccc(O)c4" +# Pregnenolone: short acetyl side chain C20–C21 +PREGNENOLONE_SMILES = ( + "CC(=O)[C@H]1CC[C@@H]2[C@@]1(CC[C@H]3[C@H]2CC=C4[C@@]3(CC[C@@H](C4)O)C)C" +) + +CORE_POSITIONS = {f"steroid_{n}" for n in range(1, 18)} +CHOLESTEROL_SIDE_CHAIN = {f"steroid_{n}" for n in range(20, 28)} + + +def _steroid_keys(smiles: str) -> set[str]: + mol = smiles_or_inchi_to_mol(smiles) + assert mol is not None + return set(get_steroid_positions(mol)) + + +class TestSteroidPositions: + def test_non_steroid_has_no_positions(self): + assert get_steroid_positions(smiles_or_inchi_to_mol("c1ccccc1")) == {} + + def test_cholesterol_has_core_angular_methyls_and_side_chain(self): + keys = _steroid_keys(CHOLESTEROL_SMILES) + assert CORE_POSITIONS <= keys + assert "steroid_18" in keys + assert "steroid_19" in keys + assert CHOLESTEROL_SIDE_CHAIN <= keys + + def test_estrone_has_c18_but_not_c19_or_side_chain(self): + keys = _steroid_keys(ESTRONE_SMILES) + assert CORE_POSITIONS <= keys + assert "steroid_18" in keys + assert "steroid_19" not in keys + assert not any(int(key.split("_")[1]) >= 20 for key in keys) + + def test_pregnenolone_has_c20_and_c21_only(self): + keys = _steroid_keys(PREGNENOLONE_SMILES) + assert CORE_POSITIONS <= keys + assert "steroid_20" in keys + assert "steroid_21" in keys + assert "steroid_22" not in keys + + @pytest.mark.parametrize( + "smiles,attachment,methyl", + [ + (CHOLESTEROL_SMILES, "steroid_13", "steroid_18"), + (CHOLESTEROL_SMILES, "steroid_10", "steroid_19"), + (ESTRONE_SMILES, "steroid_13", "steroid_18"), + ], + ) + def test_angular_methyl_is_bonded_to_attachment_atom(self, smiles, attachment, methyl): + mol = smiles_or_inchi_to_mol(smiles) + positions = get_steroid_positions(mol) + attachment_idx = positions[attachment][0] + methyl_idx = positions[methyl][0] + attachment_atom = mol.GetAtomWithIdx(attachment_idx) + neighbor_indices = {neighbor.GetIdx() for neighbor in attachment_atom.GetNeighbors()} + assert methyl_idx in neighbor_indices + assert mol.GetAtomWithIdx(methyl_idx).GetAtomicNum() == 6 + + @pytest.mark.parametrize( + "smiles,parent,child", + [ + (CHOLESTEROL_SMILES, "steroid_17", "steroid_20"), + (CHOLESTEROL_SMILES, "steroid_20", "steroid_21"), + (CHOLESTEROL_SMILES, "steroid_20", "steroid_22"), + (CHOLESTEROL_SMILES, "steroid_25", "steroid_26"), + (CHOLESTEROL_SMILES, "steroid_25", "steroid_27"), + (PREGNENOLONE_SMILES, "steroid_17", "steroid_20"), + (PREGNENOLONE_SMILES, "steroid_20", "steroid_21"), + ], + ) + def test_side_chain_parent_child_bond(self, smiles, parent, child): + mol = smiles_or_inchi_to_mol(smiles) + positions = get_steroid_positions(mol) + parent_idx = positions[parent][0] + child_idx = positions[child][0] + parent_atom = mol.GetAtomWithIdx(parent_idx) + neighbor_indices = {neighbor.GetIdx() for neighbor in parent_atom.GetNeighbors()} + assert child_idx in neighbor_indices + + def test_mol_to_fol_atoms_includes_steroid_predicates(self): + mol = smiles_or_inchi_to_mol(CHOLESTEROL_SMILES) + atom_facts, _mol_facts = mol_to_fol_atoms(mol, with_steroids=True) + assert "steroid_1" in atom_facts + assert "steroid_18" in atom_facts + assert "steroid_19" in atom_facts + assert "steroid_20" in atom_facts + assert "steroid_27" in atom_facts + + def test_mol_to_fol_atoms_can_skip_steroids(self): + mol = smiles_or_inchi_to_mol(CHOLESTEROL_SMILES) + atom_facts, _mol_facts = mol_to_fol_atoms(mol, with_steroids=False) + assert not any(key.startswith("steroid_") for key in atom_facts)