add molecule parsing functionality

.gitignore

@@ -1,3 +1,6 @@
+# Local Claude Code notes (personal, not committed)
+CLAUDE.local.md
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[codz]

README.md

@@ -27,8 +27,8 @@ pip install -e ".[dev]"
 ```python
 from chebi_utils import download_chebi_obo, download_chebi_sdf
 
-obo_path = download_chebi_obo(version=248, dest_dir="data/")   # downloads chebi.obo
-sdf_path = download_chebi_sdf(version=248, dest_dir="data/")   # downloads chebi.sdf.gz
+obo_path = download_chebi_obo(version=248, dest_dir="data/")  # downloads chebi.obo
+sdf_path = download_chebi_sdf(version=248, dest_dir="data/")  # downloads chebi.sdf.gz
 ```
 
 A specific ChEBI release `version` (e.g. `230`, `245`, `248`) must be provided.
@@ -92,8 +92,8 @@ from chebi_utils import create_multilabel_splits
 
 splits = create_multilabel_splits(dataset, train_ratio=0.8, val_ratio=0.1, test_ratio=0.1)
 train_df = splits["train"]
-val_df   = splits["validation"]   # renamed from "val" in v0.3
-test_df  = splits["test"]
+val_df = splits["validation"]  # renamed from "val" in v0.3
+test_df = splits["test"]
 ```
 
 Columns 0 and 1 (`chebi_id`, `mol`) are treated as metadata; all remaining

chebi_utils/read_molecule.py

@@ -0,0 +1,74 @@
+import warnings
+
+from chembl_structure_pipeline.standardizer import update_mol_valences
+from rdkit import Chem
+
+
+def _sanitize_molecule(mol: Chem.Mol) -> Chem.Mol:
+    """Sanitize molecule, falling back to the unsanitized molecule on failure"""
+    try:
+        mol = update_mol_valences(mol)
+        Chem.SanitizeMol(mol)
+    except Exception as e:
+        warnings.warn(f"Failed to sanitize molecule: {e}", stacklevel=2)
+        mol.UpdatePropertyCache(strict=False)
+        Chem.FastFindRings(mol)
+    return mol
+
+
+def parse_molblock(molblock: str, chebi_id: str | None = None) -> Chem.Mol | None:
+    """Parse a V2000/V3000 molblock into an RDKit Mol object.
+
+    Sanitize molecules with the ChEMBL structure pipeline for consistency with ChEBI.
+    If sanitization fails, the unsanitized molecule is returned.
+
+    Parameters
+    ----------
+    molblock : str
+        The molblock string (header + atom/bond table + ``M  END``).
+    chebi_id : str or None
+        Used only for the warning message when parsing fails.
+
+    Returns
+    -------
+    Chem.Mol or None
+        Parsed molecule, or ``None`` if parsing failed.
+    """
+    mol = Chem.MolFromMolBlock(molblock, sanitize=False, removeHs=False)
+    if mol is None:
+        warnings.warn(f"Failed to parse molblock for {chebi_id}", stacklevel=2)
+        return None
+
+    return _sanitize_molecule(mol)
+
+
+def smiles_or_inchi_to_mol(smiles_or_inchi: str) -> Chem.rdchem.Mol | None:
+    """Parse a SMILES or InChI string into an RDKit Mol object.
+
+    Sanitize molecules with the ChEMBL structure pipeline for consistency with ChEBI.
+    If sanitization fails, the unsanitized molecule is returned.
+
+    Parameters
+    ----------
+    smiles_or_inchi : str
+        The SMILES or InChI string to parse.
+
+    Returns
+    -------
+    Chem.Mol or None
+        Parsed molecule, or ``None`` if parsing failed.
+    """
+
+    if smiles_or_inchi.startswith("InChI="):
+        mol = Chem.MolFromInchi(smiles_or_inchi, sanitize=False, removeHs=False)
+    else:
+        params = Chem.SmilesParserParams()
+        params.removeHs = False
+        params.sanitize = False
+        mol = Chem.MolFromSmiles(smiles_or_inchi, params)
+
+    if mol is None:
+        warnings.warn(f"RDKit failed at parsing {smiles_or_inchi} (returned None)", stacklevel=2)
+        return None
+
+    return _sanitize_molecule(mol)

chebi_utils/sdf_extractor.py

@@ -3,56 +3,12 @@
 from __future__ import annotations
 
 import gzip
-import warnings
 from pathlib import Path
-from typing import Optional
 
 import pandas as pd
-from rdkit import Chem
 
 from chebi_utils.obo_extractor import _chebi_id_to_str
-
-
-def _sanitize_molecule(mol: Chem.Mol) -> Optional[Chem.Mol]:
-    """Sanitize molecule"""
-    from chembl_structure_pipeline.standardizer import update_mol_valences
-
-    mol = update_mol_valences(mol)
-    try:
-        Chem.SanitizeMol(mol)
-    except Exception as e:
-        warnings.warn(f"Failed to sanitize molecule: {e}", stacklevel=2)
-        mol = None
-    return mol
-
-
-def _parse_molblock(molblock: str, chebi_id: str | None = None) -> Chem.Mol | None:
-    """Parse a V2000/V3000 molblock into an RDKit Mol object.
-
-    Uses partial sanitisation to handle ChEBI molecules with unusual valences
-    or radicals.
-
-    Parameters
-    ----------
-    molblock : str
-        The molblock string (header + atom/bond table + ``M  END``).
-    chebi_id : str or None
-        Used only for the warning message when parsing fails.
-
-    Returns
-    -------
-    Chem.Mol or None
-        Parsed molecule, or ``None`` if parsing failed.
-    """
-    mol = Chem.MolFromMolBlock(molblock, sanitize=False, removeHs=False)
-    if mol is None:
-        warnings.warn(f"Failed to parse molblock for {chebi_id}", stacklevel=2)
-        return None
-    mol = _sanitize_molecule(mol)
-    if mol is None:
-        warnings.warn(f"Failed to sanitize molblock for {chebi_id}", stacklevel=2)
-
-    return mol
+from chebi_utils.read_molecule import parse_molblock
 
 
 def _iter_sdf_records(filepath: str | Path):
@@ -164,7 +120,7 @@ def extract_molecules(filepath: str | Path) -> pd.DataFrame:
     df = df.rename(columns={k: v for k, v in rename_map.items() if k in df.columns})
 
     chebi_ids = df["chebi_id"].tolist() if "chebi_id" in df.columns else [None] * len(df)
-    df["mol"] = [_parse_molblock(mb, cid) for mb, cid in zip(molblocks, chebi_ids, strict=False)]
+    df["mol"] = [parse_molblock(mb, cid) for mb, cid in zip(molblocks, chebi_ids, strict=False)]
     df["chebi_id"] = df["chebi_id"].apply(_chebi_id_to_str)
 
     df = df[df["mol"].notna()]

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "chebi-utils"
-version = "0.3"
+version = "0.4"
 description = "Common processing functionality for the ChEBI ontology"
 readme = "README.md"
 license = { file = "LICENSE" }

tests/test_read_molecule.py

@@ -0,0 +1,172 @@
+"""Tests for chebi_utils.read_molecule."""
+
+from __future__ import annotations
+
+import pytest
+from rdkit import Chem
+from rdkit.Chem import rdchem
+
+from chebi_utils.read_molecule import parse_molblock, smiles_or_inchi_to_mol
+
+ETHANOL_SMILES = "CCO"
+ETHANOL_INCHI = "InChI=1S/C2H6O/c1-2-3/h3H,2H2,1H3"
+BENZENE_INCHI = "InChI=1S/C6H6/c1-2-4-6-5-3-1/h1-6H"
+
+ETHANOL_MOLBLOCK = """ethanol
+     RDKit          2D
+
+  3  2  0  0  0  0  0  0  0  0999 V2000
+    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.2990    0.7500    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.5981   -0.0000    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+M  END
+"""
+
+
+class TestSmilesToMol:
+    @pytest.mark.parametrize(
+        "smiles,num_atoms",
+        [
+            ("C", 1),  # methane
+            (ETHANOL_SMILES, 3),
+            ("c1ccccc1", 6),  # benzene
+            ("[Na+].[Cl-]", 2),  # multi-component
+            ("O=C(C)Oc1ccccc1C(=O)O", 13),  # aspirin
+        ],
+    )
+    def test_valid_smiles_atom_counts(self, smiles, num_atoms):
+        mol = smiles_or_inchi_to_mol(smiles)
+        assert isinstance(mol, rdchem.Mol)
+        assert mol.GetNumAtoms() == num_atoms
+
+    def test_canonical_smiles_roundtrip(self):
+        mol = smiles_or_inchi_to_mol("OCC")
+        assert Chem.MolToSmiles(mol) == "CCO"
+
+    def test_aromaticity_perceived(self):
+        # Kekulized input must come back aromatic, i.e. sanitization ran
+        mol = smiles_or_inchi_to_mol("C1=CC=CC=C1")
+        assert all(atom.GetIsAromatic() for atom in mol.GetAtoms())
+
+    def test_ring_info_available(self):
+        mol = smiles_or_inchi_to_mol("C1CCCCC1")
+        assert mol.GetRingInfo().NumRings() == 1
+
+    def test_explicit_hydrogens_are_kept(self):
+        # removeHs=False, so the explicit H atoms stay in the graph
+        mol = smiles_or_inchi_to_mol("[H]C([H])([H])[H]")
+        assert mol.GetNumAtoms() == 5
+
+    def test_nitro_group_normalized_by_chembl_pipeline(self):
+        # The ChEMBL structure pipeline rewrites N(=O)=O as the charge-separated form
+        mol = smiles_or_inchi_to_mol("CN(=O)=O")
+        assert Chem.MolToSmiles(mol) == "C[N+](=O)[O-]"
+
+    def test_empty_string_gives_empty_mol(self):
+        mol = smiles_or_inchi_to_mol("")
+        assert mol is not None
+        assert mol.GetNumAtoms() == 0
+
+    @pytest.mark.parametrize(
+        "smiles",
+        [
+            "not_a_molecule",
+            "C(C",  # unclosed branch
+            "C1CC",  # unclosed ring
+            "[CH3",  # unclosed bracket atom
+            "%%%",
+        ],
+    )
+    def test_invalid_smiles_returns_none(self, smiles):
+        with pytest.warns(UserWarning, match="RDKit failed at parsing"):
+            assert smiles_or_inchi_to_mol(smiles) is None
+
+
+class TestInchiToMol:
+    @pytest.mark.parametrize(
+        "inchi,num_atoms",
+        [
+            ("InChI=1S/CH4/h1H4", 1),  # methane
+            ("InChI=1S/H2O/h1H2", 1),  # water
+            (ETHANOL_INCHI, 3),
+            (BENZENE_INCHI, 6),
+        ],
+    )
+    def test_valid_inchi_atom_counts(self, inchi, num_atoms):
+        mol = smiles_or_inchi_to_mol(inchi)
+        assert isinstance(mol, rdchem.Mol)
+        assert mol.GetNumAtoms() == num_atoms
+
+    def test_inchi_and_smiles_agree(self):
+        from_inchi = smiles_or_inchi_to_mol(ETHANOL_INCHI)
+        from_smiles = smiles_or_inchi_to_mol(ETHANOL_SMILES)
+        assert Chem.MolToSmiles(from_inchi) == Chem.MolToSmiles(from_smiles)
+
+    def test_inchi_aromaticity_perceived(self):
+        mol = smiles_or_inchi_to_mol(BENZENE_INCHI)
+        assert Chem.MolToSmiles(mol) == "c1ccccc1"
+
+    @pytest.mark.parametrize(
+        "inchi",
+        [
+            "InChI=1S/garbage",
+            "InChI=",
+            "InChI=1S/CH4N/c1-2(1)1/h1H3",  # syntactically well-formed but not a real layer
+        ],
+    )
+    def test_invalid_inchi_returns_none(self, inchi):
+        with pytest.warns(UserWarning, match="RDKit failed at parsing"):
+            assert smiles_or_inchi_to_mol(inchi) is None
+
+    def test_inchi_without_prefix_is_treated_as_smiles(self):
+        # Dispatch is based on the "InChI=" prefix, so a bare InChI body fails
+        with pytest.warns(UserWarning, match="RDKit failed at parsing"):
+            assert smiles_or_inchi_to_mol("1S/C2H6O/c1-2-3/h3H,2H2,1H3") is None
+
+
+class TestSanitizationFallback:
+    def test_unsanitizable_molecule_is_still_returned(self):
+        # Neutral 5-valent nitrogen: sanitization fails, the raw mol comes back
+        with pytest.warns(UserWarning, match="Failed to sanitize molecule"):
+            mol = smiles_or_inchi_to_mol("C[N](C)(C)C")
+        assert mol is not None
+        assert mol.GetNumAtoms() == 5
+
+    def test_ring_info_available_after_failed_sanitization(self):
+        # FastFindRings runs in the fallback path, so ring queries do not raise
+        with pytest.warns(UserWarning, match="Failed to sanitize molecule"):
+            mol = smiles_or_inchi_to_mol("C1CCCCC1[N](C)(C)C")
+        assert mol.GetRingInfo().NumRings() == 1
+
+    def test_valid_molecule_emits_no_warning(self, recwarn):
+        smiles_or_inchi_to_mol(ETHANOL_SMILES)
+        assert len(recwarn) == 0
+
+
+class TestParseMolblock:
+    def test_valid_molblock(self):
+        mol = parse_molblock(ETHANOL_MOLBLOCK, "CHEBI:16236")
+        assert isinstance(mol, rdchem.Mol)
+        assert mol.GetNumAtoms() == 3
+        assert Chem.MolToSmiles(mol) == "CCO"
+
+    def test_valid_molblock_without_chebi_id(self):
+        mol = parse_molblock(ETHANOL_MOLBLOCK)
+        assert mol is not None
+        assert mol.GetNumAtoms() == 3
+
+    @pytest.mark.parametrize("molblock", ["", "nonsense", "one\ntwo\nthree\n"])
+    def test_invalid_molblock_returns_none(self, molblock):
+        with pytest.warns(UserWarning, match="Failed to parse molblock"):
+            assert parse_molblock(molblock, "CHEBI:99") is None
+
+    def test_warning_mentions_chebi_id(self):
+        with pytest.warns(UserWarning, match="CHEBI:12345"):
+            parse_molblock("nonsense", "CHEBI:12345")
+
+    def test_atomless_molblock_returns_none(self):
+        molblock = "empty\n\n  0  0  0  0  0  0  0  0  0  0999 V2000\nM  END\n"
+        with pytest.warns(UserWarning, match="Failed to parse molblock"):
+            assert parse_molblock(molblock, "CHEBI:192499") is None