diff --git a/.gitignore b/.gitignore index c2d747d..666b124 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ +# Local Claude Code notes (personal, not committed) +CLAUDE.local.md + # Byte-compiled / optimized / DLL files __pycache__/ *.py[codz] diff --git a/README.md b/README.md index e79034d..53828f0 100644 --- a/README.md +++ b/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 diff --git a/chebi_utils/read_molecule.py b/chebi_utils/read_molecule.py new file mode 100644 index 0000000..67885ae --- /dev/null +++ b/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) diff --git a/chebi_utils/sdf_extractor.py b/chebi_utils/sdf_extractor.py index 22aae48..2843d10 100644 --- a/chebi_utils/sdf_extractor.py +++ b/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()] diff --git a/pyproject.toml b/pyproject.toml index 7a6c9d6..ed59c4b 100644 --- a/pyproject.toml +++ b/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" } diff --git a/tests/test_read_molecule.py b/tests/test_read_molecule.py new file mode 100644 index 0000000..cd50ef3 --- /dev/null +++ b/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