Feature/steroid positions

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)