From f61a80c4e520193794a539e30ea51163b0caf034 Mon Sep 17 00:00:00 2001 From: kzubko Date: Mon, 17 Aug 2026 17:30:31 +0200 Subject: [PATCH 1/5] Add C18/C19 steroid labels, tests, and a small demo --- chebi_utils/extract_properties.py | 69 +++++++++++++++++++++++++++--- tests/test_extract_properties.py | 70 +++++++++++++++++++++++++++++++ 2 files changed, 133 insertions(+), 6 deletions(-) create mode 100644 tests/test_extract_properties.py diff --git a/chebi_utils/extract_properties.py b/chebi_utils/extract_properties.py index c96f351..25418ae 100644 --- a/chebi_utils/extract_properties.py +++ b/chebi_utils/extract_properties.py @@ -238,12 +238,41 @@ 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] + ) + + 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). Molecules without a gonane core yield + no predicates. Parameters ---------- @@ -258,11 +287,17 @@ def get_steroid_positions(mol: Chem.Mol) -> dict[str, list]: """ 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) return atom_extensions @@ -289,3 +324,25 @@ 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..19, estrone -> steroid_1..18, 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", + "benzene": "c1ccccc1", + }.items(): + keys = sorted( + get_steroid_positions(smiles_or_inchi_to_mol(smiles)), + key=lambda k: int(k.split("_")[1]), + ) + print(name, keys) diff --git a/tests/test_extract_properties.py b/tests/test_extract_properties.py new file mode 100644 index 0000000..6c22080 --- /dev/null +++ b/tests/test_extract_properties.py @@ -0,0 +1,70 @@ +"""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 and C19 angular methyls +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 absent +ESTRONE_SMILES = "C[C@]12CC[C@H]3[C@H]([C@@H]1CCC2=O)CCc4c3ccc(O)c4" + +CORE_POSITIONS = {f"steroid_{n}" for n in range(1, 18)} + + +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_and_angular_methyls(self): + keys = _steroid_keys(CHOLESTEROL_SMILES) + assert CORE_POSITIONS <= keys + assert "steroid_18" in keys + assert "steroid_19" in keys + + def test_estrone_has_c18_but_not_c19(self): + keys = _steroid_keys(ESTRONE_SMILES) + assert CORE_POSITIONS <= keys + assert "steroid_18" in keys + assert "steroid_19" 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 + + 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 + + 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) From 4f4a8ecd3b9f48c39cd37181a16aa98ddcfdaba0 Mon Sep 17 00:00:00 2001 From: kzubko Date: Mon, 24 Aug 2026 15:44:38 +0200 Subject: [PATCH 2/5] Add C20+ steroid side-chain labels, tests, and demo --- chebi_utils/extract_properties.py | 164 +++++++++++++++++++++++++++++- tests/test_extract_properties.py | 45 +++++++- 2 files changed, 202 insertions(+), 7 deletions(-) diff --git a/chebi_utils/extract_properties.py b/chebi_utils/extract_properties.py index 25418ae..cc01647 100644 --- a/chebi_utils/extract_properties.py +++ b/chebi_utils/extract_properties.py @@ -265,14 +265,167 @@ def _add_angular_methyls( ) + +# --- 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 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``. When present, angular methyls are added as ``steroid_18`` - (on C13) and ``steroid_19`` (on C10). Molecules without a gonane core yield - no predicates. + (on C13) and ``steroid_19`` (on C10), and the C17 side chain as + ``steroid_20``, ``steroid_21``, … (cholestane-style). Molecules without a + gonane core yield no predicates. Parameters ---------- @@ -298,6 +451,7 @@ def get_steroid_positions(mol: Chem.Mol) -> dict[str, list]: 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+ return atom_extensions @@ -329,7 +483,8 @@ def get_numerical_facts(mol: Chem.Mol) -> dict[str, list]: """Manual check for steroid numbering (not part of the library API). Run: python -m chebi_utils.extract_properties -Expect: cholesterol -> steroid_1..19, estrone -> steroid_1..18, benzene -> [] +Expect: cholesterol -> steroid_1..27, estrone -> steroid_1..18, +pregnenolone -> steroid_1..21 (with 18/19), benzene -> [] """ if __name__ == "__main__": from chebi_utils.read_molecule import smiles_or_inchi_to_mol @@ -339,6 +494,9 @@ def get_numerical_facts(mol: Chem.Mol) -> dict[str, list]: "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" + ), "benzene": "c1ccccc1", }.items(): keys = sorted( diff --git a/tests/test_extract_properties.py b/tests/test_extract_properties.py index 6c22080..d5dac22 100644 --- a/tests/test_extract_properties.py +++ b/tests/test_extract_properties.py @@ -7,14 +7,19 @@ 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 and C19 angular methyls +# 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 absent +# 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]: @@ -27,17 +32,26 @@ 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_and_angular_methyls(self): + 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(self): + 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", @@ -57,12 +71,35 @@ def test_angular_methyl_is_bonded_to_attachment_atom(self, smiles, attachment, m 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) From ea6e576403a709c8b528fb79da925a91f58bc98e Mon Sep 17 00:00:00 2001 From: kzubko Date: Mon, 24 Aug 2026 16:17:22 +0200 Subject: [PATCH 3/5] Add helper to print steroid core CIP R/S labels Co-authored-by: Cursor --- chebi_utils/steroid_cip_ab.py | 63 +++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 chebi_utils/steroid_cip_ab.py diff --git a/chebi_utils/steroid_cip_ab.py b/chebi_utils/steroid_cip_ab.py new file mode 100644 index 0000000..00d8496 --- /dev/null +++ b/chebi_utils/steroid_cip_ab.py @@ -0,0 +1,63 @@ +"""Print CIP R/S at steroid core positions for a list of molecules. + +Run: python -m chebi_utils.steroid_cip_ab +""" + +from rdkit import Chem + +from chebi_utils.extract_properties import get_steroid_positions +from chebi_utils.read_molecule import smiles_or_inchi_to_mol + +STEREO_CORE_POSITIONS = (5, 8, 9, 10, 13, 14, 17) + + +def get_steroid_cip_fingerprint(mol: Chem.Mol) -> dict[int, str | None]: + """CIP labels at C5/C8/C9/C10/C13/C14/C17 (None if missing/achiral).""" + positions = get_steroid_positions(mol) + try: + Chem.rdCIPLabeler.AssignCIPLabels(mol) + except Exception: + pass + + fingerprint: dict[int, str | None] = {} + for iupac_position in STEREO_CORE_POSITIONS: + atom_indices = positions.get(f"steroid_{iupac_position}") + if not atom_indices: + fingerprint[iupac_position] = None + continue + atom = mol.GetAtomWithIdx(atom_indices[0]) + fingerprint[iupac_position] = ( + atom.GetProp("_CIPCode") if atom.HasProp("_CIPCode") else None + ) + return fingerprint + + +# 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"), +] + + +if __name__ == "__main__": + for name, smiles_or_inchi, label in EXAMPLES: + mol = smiles_or_inchi_to_mol(smiles_or_inchi) + if mol is None: + print(f"{name}: parse_failed") + continue + fingerprint = get_steroid_cip_fingerprint(mol) + cip_text = " ".join( + f"C{position}={fingerprint.get(position) or '-'}" + for position in STEREO_CORE_POSITIONS + ) + print(f"{name} [{label}]: {cip_text}") From a3f2031f30cfa69bb6d503df1e37b7398e1f5ec5 Mon Sep 17 00:00:00 2001 From: kzubko Date: Sun, 30 Aug 2026 20:12:24 +0200 Subject: [PATCH 4/5] Add C5 stroid parity (+1/-1) demo --- chebi_utils/steroid_cip_ab.py | 123 ++++++++++++++++++++++++++-------- 1 file changed, 95 insertions(+), 28 deletions(-) diff --git a/chebi_utils/steroid_cip_ab.py b/chebi_utils/steroid_cip_ab.py index 00d8496..a7c73c7 100644 --- a/chebi_utils/steroid_cip_ab.py +++ b/chebi_utils/steroid_cip_ab.py @@ -1,4 +1,7 @@ -"""Print CIP R/S at steroid core positions for a list of molecules. +"""Explore C5 steroid parity (+1/-1) vs known alpha/beta labels. + +Uses RDKit tetrahedral CW/CCW remapped to a fixed steroid neighbor order +(C4, C6, C10), not CIP. Run: python -m chebi_utils.steroid_cip_ab """ @@ -8,28 +11,80 @@ from chebi_utils.extract_properties import get_steroid_positions from chebi_utils.read_molecule import smiles_or_inchi_to_mol -STEREO_CORE_POSITIONS = (5, 8, 9, 10, 13, 14, 17) - -def get_steroid_cip_fingerprint(mol: Chem.Mol) -> dict[int, str | None]: - """CIP labels at C5/C8/C9/C10/C13/C14/C17 (None if missing/achiral).""" +def _iupac_to_atom(mol: Chem.Mol) -> dict[int, int]: positions = get_steroid_positions(mol) + return { + int(key.split("_")[1]): atom_indices[0] + for key, atom_indices in positions.items() + if key.startswith("steroid_") + } + + +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. + + Returns ``+1`` / ``-1``, or ``None`` if C5 is missing / not tetrahedral / + neighbours are unexpected. Does not name alpha/beta yet. + """ try: - Chem.rdCIPLabeler.AssignCIPLabels(mol) - except Exception: - pass - - fingerprint: dict[int, str | None] = {} - for iupac_position in STEREO_CORE_POSITIONS: - atom_indices = positions.get(f"steroid_{iupac_position}") - if not atom_indices: - fingerprint[iupac_position] = None - continue - atom = mol.GetAtomWithIdx(atom_indices[0]) - fingerprint[iupac_position] = ( - atom.GetProp("_CIPCode") if atom.HasProp("_CIPCode") else None - ) - return fingerprint + 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 # Paste substances here: (name, SMILES_or_InChI, "alpha"|"beta") @@ -37,7 +92,7 @@ def get_steroid_cip_fingerprint(mol: Chem.Mol) -> dict[int, str | None]: ("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"), ("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"), @@ -46,18 +101,30 @@ def get_steroid_cip_fingerprint(mol: Chem.Mol) -> dict[int, str | None]: ("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, label in EXAMPLES: mol = smiles_or_inchi_to_mol(smiles_or_inchi) if mol is None: - print(f"{name}: parse_failed") + rows.append((name, label, "parse_failed")) continue - fingerprint = get_steroid_cip_fingerprint(mol) - cip_text = " ".join( - f"C{position}={fingerprint.get(position) or '-'}" - for position in STEREO_CORE_POSITIONS - ) - print(f"{name} [{label}]: {cip_text}") + parity = c5_steroid_parity(mol, _iupac_to_atom(mol)) + parity_text = "-" if parity is None else f"{parity:+d}" + rows.append((name, label, parity_text)) + + name_width = max(len(name) for name, _, _ in rows) + label_width = max(len("expected"), max(len(label) for _, label, _ in rows)) + print(f"{'molecule':<{name_width}} | {'expected':<{label_width}} | parity") + print(f"{'-' * name_width}--+--{'-' * label_width}--+--------") + for name, label, parity_text in rows: + print(f"{name:<{name_width}} | {label:<{label_width}} | {parity_text}") From f764d0d1995d583da5f741253fbeffb2aa86b2b3 Mon Sep 17 00:00:00 2001 From: kzubko Date: Sun, 30 Aug 2026 20:55:30 +0200 Subject: [PATCH 5/5] Add stroid c5 alpha/beta labels --- chebi_utils/extract_properties.py | 114 +++++++++++++++++++++++++++--- chebi_utils/steroid_cip_ab.py | 104 ++++----------------------- 2 files changed, 121 insertions(+), 97 deletions(-) diff --git a/chebi_utils/extract_properties.py b/chebi_utils/extract_properties.py index cc01647..a012219 100644 --- a/chebi_utils/extract_properties.py +++ b/chebi_utils/extract_properties.py @@ -417,15 +417,102 @@ def _add_side_chain( 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``. When present, angular methyls are added as ``steroid_18`` - (on C13) and ``steroid_19`` (on C10), and the C17 side chain as - ``steroid_20``, ``steroid_21``, … (cholestane-style). Molecules without a - gonane core yield no predicates. + (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 ---------- @@ -435,8 +522,9 @@ 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) @@ -452,6 +540,7 @@ def get_steroid_positions(mol: Chem.Mol) -> dict[str, list]: _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 @@ -483,8 +572,9 @@ def get_numerical_facts(mol: Chem.Mol) -> dict[str, list]: """Manual check for steroid numbering (not part of the library API). Run: python -m chebi_utils.extract_properties -Expect: cholesterol -> steroid_1..27, estrone -> steroid_1..18, -pregnenolone -> steroid_1..21 (with 18/19), benzene -> [] +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 @@ -497,10 +587,18 @@ def get_numerical_facts(mol: Chem.Mol) -> dict[str, list]: "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]), + 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 index a7c73c7..2cefbf2 100644 --- a/chebi_utils/steroid_cip_ab.py +++ b/chebi_utils/steroid_cip_ab.py @@ -1,90 +1,18 @@ -"""Explore C5 steroid parity (+1/-1) vs known alpha/beta labels. - -Uses RDKit tetrahedral CW/CCW remapped to a fixed steroid neighbor order -(C4, C6, C10), not CIP. +"""Demo: C5 alpha/beta from get_steroid_positions. Run: python -m chebi_utils.steroid_cip_ab """ -from rdkit import Chem - from chebi_utils.extract_properties import get_steroid_positions from chebi_utils.read_molecule import smiles_or_inchi_to_mol -def _iupac_to_atom(mol: Chem.Mol) -> dict[int, int]: - positions = get_steroid_positions(mol) - return { - int(key.split("_")[1]): atom_indices[0] - for key, atom_indices in positions.items() - if key.startswith("steroid_") - } - - -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. - - Returns ``+1`` / ``-1``, or ``None`` if C5 is missing / not tetrahedral / - neighbours are unexpected. Does not name alpha/beta yet. - """ - 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 _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") @@ -113,18 +41,16 @@ def c5_steroid_parity(mol: Chem.Mol, iupac_to_atom: dict[int, int]) -> int | Non if __name__ == "__main__": rows: list[tuple[str, str, str]] = [] - for name, smiles_or_inchi, label in EXAMPLES: + for name, smiles_or_inchi, expected in EXAMPLES: mol = smiles_or_inchi_to_mol(smiles_or_inchi) if mol is None: - rows.append((name, label, "parse_failed")) + rows.append((name, expected, "parse_failed")) continue - parity = c5_steroid_parity(mol, _iupac_to_atom(mol)) - parity_text = "-" if parity is None else f"{parity:+d}" - rows.append((name, label, parity_text)) + rows.append((name, expected, _c5_ab_label(get_steroid_positions(mol)))) name_width = max(len(name) for name, _, _ in rows) - label_width = max(len("expected"), max(len(label) for _, label, _ in rows)) - print(f"{'molecule':<{name_width}} | {'expected':<{label_width}} | parity") - print(f"{'-' * name_width}--+--{'-' * label_width}--+--------") - for name, label, parity_text in rows: - print(f"{name:<{name_width}} | {label:<{label_width}} | {parity_text}") + 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}")