fetching residue SMILES

Author: stefdoerr

moleculekit/rcsb.py

@@ -4,6 +4,8 @@
 # No redistribution in whole or part
 #
 import urllib.request
+import urllib.error
+import json
 import logging
 
 logger = logging.getLogger(__name__)
@@ -123,3 +125,136 @@ def rcsbFindLigands(pdbid: str) -> list:
                 name = td[0].find_all("a")[0].text.strip()
                 ligands.append(name)
     return ligands
+
+
+def _getRCSBjson(url, attempts=3):
+    import time
+
+    last_err = None
+    for _ in range(attempts):
+        try:
+            response = urllib.request.urlopen(url)
+            return json.loads(response.read())
+        except urllib.error.HTTPError as err:
+            # 404 etc. — the code is wrong, do not retry
+            raise RuntimeError(f"RCSB request failed for {url}: {err}") from err
+        except Exception as err:
+            last_err = err
+            logger.warning(
+                f"Failed to connect to URL {url} with error {err}. Sleeping 5s and retrying."
+            )
+            time.sleep(5)
+    raise RuntimeError(f"Failed to connect to URL {url}: {last_err}")
+
+
+def rcsbFetchLigandInfo(comp_id: str) -> dict:
+    """Fetch the full RCSB Chemical Component Dictionary record for a ligand.
+
+    Queries the RCSB data API for a 3-letter chemical component (CCD) code and
+    returns the complete record, including identifiers, formula, weight and all
+    descriptor variants (InChI plus SMILES from RCSB, CACTVS, OpenEye and ACDLabs).
+
+    Parameters
+    ----------
+    comp_id : str
+        The chemical component (CCD) 3-letter code, e.g. ``"BEN"``. Case-insensitive.
+
+    Returns
+    -------
+    info : dict
+        The parsed JSON record. The curated descriptors live under
+        ``info["rcsb_chem_comp_descriptor"]`` (``SMILES``, ``SMILES_stereo``,
+        ``InChI``, ``InChIKey``); per-program variants live under
+        ``info["pdbx_chem_comp_descriptor"]``.
+
+    Examples
+    --------
+    >>> info = rcsbFetchLigandInfo('BEN')
+    >>> info['rcsb_chem_comp_descriptor']['comp_id']
+    'BEN'
+    """
+    comp_id = comp_id.strip().upper()
+    url = f"https://data.rcsb.org/rest/v1/core/chemcomp/{comp_id}"
+    return _getRCSBjson(url)
+
+
+def rcsbFetchLigandSmiles(
+    comp_id: str, stereo: bool = True, program: str = "OpenEye"
+) -> str:
+    """Fetch a SMILES string for a ligand by its RCSB CCD code.
+
+    Thin wrapper over :func:`rcsbFetchLigandInfo`. RCSB stores SMILES computed by
+    several toolkits (OpenEye, CACTVS, ACDLabs). By default this returns the
+    OpenEye descriptor, which RCSB also curates into its top-level
+    ``rcsb_chem_comp_descriptor`` block. Pass ``program`` to pick a different
+    toolkit; for full control read ``pdbx_chem_comp_descriptor`` off
+    :func:`rcsbFetchLigandInfo` directly.
+
+    Parameters
+    ----------
+    comp_id : str
+        The chemical component (CCD) 3-letter code, e.g. ``"BEN"``. Case-insensitive.
+    stereo : bool
+        If True (default) return the isomeric SMILES (stereochemistry included);
+        if False return the plain SMILES. RCSB labels the isomeric variant
+        ``SMILES_CANONICAL``. Falls back to the other variant when the preferred
+        one is absent for the chosen program.
+    program : str
+        Which toolkit's descriptor to return. ``"OpenEye"`` (default) uses RCSB's
+        curated descriptor. Other typical values are ``"CACTVS"`` and
+        ``"ACDLabs"``. Matched case-insensitively as a substring of the program
+        name reported by RCSB; raises if the component has no SMILES from a
+        matching program.
+
+    Returns
+    -------
+    smiles : str
+        The SMILES string.
+
+    Examples
+    --------
+    >>> rcsbFetchLigandSmiles('BEN', stereo=False)
+    '[H]N=C(c1ccccc1)N'
+    >>> rcsbFetchLigandSmiles('BEN', program='CACTVS')
+    'NC(=N)c1ccccc1'
+    """
+    info = rcsbFetchLigandInfo(comp_id)
+    code = comp_id.strip().upper()
+    want = program.strip().lower()
+
+    # OpenEye is the default, and is exactly what RCSB curates into the top-level
+    # ``rcsb_chem_comp_descriptor`` block (always present) — use it directly.
+    if want in ("openeye", "openeye oetoolkits", "oe"):
+        desc = info.get("rcsb_chem_comp_descriptor", {})
+        primary, secondary = (
+            ("SMILES_stereo", "SMILES") if stereo else ("SMILES", "SMILES_stereo")
+        )
+        smiles = desc.get(primary) or desc.get(secondary)
+        if smiles:
+            return smiles
+        # else fall through to the per-program rows below
+
+    # Per-program descriptors (CACTVS, ACDLabs, or an OpenEye fallback).
+    rows = [
+        r
+        for r in info.get("pdbx_chem_comp_descriptor", [])
+        if "SMILES" in (r.get("type") or "")
+    ]
+    available = sorted({r.get("program") for r in rows if r.get("program")})
+    matches = [r for r in rows if want and want in (r.get("program") or "").lower()]
+    if not matches:
+        raise RuntimeError(
+            f"RCSB has no SMILES for component '{code}' from program '{program}'. "
+            f"Available programs: {available}"
+        )
+    # ``SMILES_CANONICAL`` is the isomeric (stereo-bearing) variant.
+    primary, secondary = (
+        ("SMILES_CANONICAL", "SMILES") if stereo else ("SMILES", "SMILES_CANONICAL")
+    )
+    by_type = {r.get("type"): r.get("descriptor") for r in matches}
+    smiles = by_type.get(primary) or by_type.get(secondary)
+    if not smiles:
+        raise RuntimeError(
+            f"RCSB returned no SMILES descriptor for component '{code}' from program '{program}'"
+        )
+    return smiles

tests/test_rcsb.py

@@ -0,0 +1,59 @@
+import pytest
+from moleculekit.rcsb import rcsbFetchLigandInfo, rcsbFetchLigandSmiles
+
+
+def test_fetch_ligand_info_returns_record():
+    info = rcsbFetchLigandInfo("BEN")
+    assert isinstance(info, dict)
+    assert "rcsb_chem_comp_descriptor" in info
+    assert info["rcsb_chem_comp_descriptor"]["comp_id"] == "BEN"
+    # the full record exposes per-program SMILES variants too
+    assert "pdbx_chem_comp_descriptor" in info
+
+
+def test_fetch_ligand_smiles_stereo_default():
+    # benzamidine, stereo SMILES carries the /N=C(\...)/ double-bond geometry
+    smi = rcsbFetchLigandSmiles("BEN")
+    assert isinstance(smi, str) and len(smi) > 0
+    assert "c1ccccc1" in smi
+
+
+def test_fetch_ligand_smiles_non_stereo():
+    smi = rcsbFetchLigandSmiles("BEN", stereo=False)
+    assert smi == "[H]N=C(c1ccccc1)N"
+
+
+def test_fetch_ligand_info_unknown_code_raises():
+    with pytest.raises(RuntimeError):
+        rcsbFetchLigandInfo("ZZZZ")
+
+
+def test_lowercase_code_is_accepted():
+    smi = rcsbFetchLigandSmiles("ben")
+    assert "c1ccccc1" in smi
+
+
+def test_default_program_is_openeye():
+    # explicit OpenEye matches the default (curated rcsb_chem_comp_descriptor)
+    assert rcsbFetchLigandSmiles("BEN", program="OpenEye") == rcsbFetchLigandSmiles("BEN")
+
+
+def test_program_cactvs():
+    assert rcsbFetchLigandSmiles("BEN", program="CACTVS") == "NC(=N)c1ccccc1"
+    # CACTVS differs from the OpenEye default for this ligand
+    assert rcsbFetchLigandSmiles("BEN", program="CACTVS") != rcsbFetchLigandSmiles("BEN")
+
+
+def test_program_match_is_case_insensitive():
+    assert rcsbFetchLigandSmiles("BEN", program="cactvs") == "NC(=N)c1ccccc1"
+
+
+def test_program_falls_back_to_other_type_when_canonical_absent():
+    # ACDLabs only provides a plain SMILES row for BEN (no SMILES_CANONICAL),
+    # so a stereo request falls back to it rather than raising.
+    assert rcsbFetchLigandSmiles("BEN", program="ACDLabs") == "[N@H]=C(N)c1ccccc1"
+
+
+def test_unknown_program_raises():
+    with pytest.raises(RuntimeError):
+        rcsbFetchLigandSmiles("BEN", program="Nonesuch")