Improve small molecule dictionary file readers

pdbeccdutils/core/ccd_reader.py

@@ -28,7 +28,7 @@
 import logging
 import os
 from datetime import date
-from typing import Dict, List, NamedTuple
+from typing import Collection, Dict, List, NamedTuple, Optional
 
 import rdkit
 from pdbeccdutils.core.component import Component
@@ -101,7 +101,7 @@ def read_pdb_cif_file(path_to_cif: str, sanitize: bool = True) -> CCDReaderResul
 
 
 def read_pdb_components_file(
-    path_to_cif: str, sanitize: bool = True, include: list[str] = []
+    path_to_cif: str, sanitize: bool = True, include: Optional[Collection[str]] = None
 ) -> Dict[str, CCDReaderResult]:
     """
     Process multiple compounds stored in the wwPDB CCD
@@ -112,8 +112,10 @@ def read_pdb_components_file(
             multiple ligands in it.
         sanitize (bool): Whether or not the components should be sanitized
             Defaults to True.
-        include (list[str]): List of CCDs to be parsed. By default it is empty and parse
-        all the CCDs. If a list of CCDs provided, will only parse them
+        include (Optional[Collection[str]]): List of CCDs to be parsed. By default it
+            is None and parses all the CCDs. If a collection of CCDs is provided, only
+            those CCDs are parsed. Empty collections parse all CCDs, matching the
+            previous behavior.
 
     Raises:
         ValueError: if the file does not exist.
@@ -126,16 +128,27 @@ def read_pdb_components_file(
         raise ValueError("File '{}' does not exists".format(path_to_cif))
 
     result_bag = {}
-    for block in cif.read(path_to_cif):
-        if (block.name in include) or (len(include) == 0):
-            try:
-                result_bag[block.name] = _parse_pdb_mmcif(block, sanitize)
-            except CCDUtilsError as e:
-                logging.error(
-                    f"ERROR: Data block {block.name} not processed. Reason: ({str(e)})."
-                )
-    return result_bag
+    doc = cif.read(path_to_cif)
 
+    if not include:
+        blocks = doc
+    else:
+        blocks = []
+        for block_name in dict.fromkeys(include):
+            block = doc.find_block(block_name)
+            if block is None:
+                logging.warning(f"Data block {block_name} not found in {path_to_cif}.")
+                continue
+            blocks.append(block)
+
+    for block in blocks:
+        try:
+            result_bag[block.name] = _parse_pdb_mmcif(block, sanitize)
+        except CCDUtilsError as e:
+            logging.error(
+                f"ERROR: Data block {block.name} not processed. Reason: ({str(e)})."
+            )
+    return result_bag
 
 # region parse mmcif
 
@@ -316,19 +329,31 @@ def _parse_pdb_bonds(mol, cif_block, errors):
     bonds = cif_block.find(
         "_chem_comp_bond.", ["atom_id_1", "atom_id_2", "value_order"]
     )
-    atoms_ids = list(atoms.find_column("_chem_comp_atom.atom_id"))
-    for row in bonds:
+    atom_indices = {
+        atom_id: atom_index
+        for atom_index, _atom_id in enumerate(
+            atoms.find_column("_chem_comp_atom.atom_id")
+        ) if (atom_id := cif.as_string(_atom_id))
+    }
+    for index, row in enumerate(bonds):
+        atom_1 = cif.as_string(row["_chem_comp_bond.atom_id_1"])
+        atom_2 = cif.as_string(row["_chem_comp_bond.atom_id_2"])
+
         try:
-            atom_1 = row["_chem_comp_bond.atom_id_1"]
-            atom_1_id = atoms_ids.index(atom_1)
-            atom_2 = row["_chem_comp_bond.atom_id_2"]
-            atom_2_id = atoms_ids.index(atom_2)
-            bond_order = helper.bond_pdb_order(row["_chem_comp_bond.value_order"])
+            atom_1_id = atom_indices[atom_1]
+            atom_2_id = atom_indices[atom_2]
+            bond_order = helper.bond_pdb_order(cif.as_string(row["_chem_comp_bond.value_order"]))
+            if bond_order is None:
+                errors.append(
+                    f"""Unknown bond order {cif.as_string(row['_chem_comp_bond.value_order'])} in 
+                    {index} entry in _chem_comp_bond"""
+                )
+                continue
 
             mol.AddBond(atom_1_id, atom_2_id, bond_order)
-        except ValueError:
+        except KeyError:
             errors.append(
-                f"Error perceiving {atom_1} - {atom_2} bond in _chem_comp_bond"
+                f"Missing atom in {index} entry in _chem_comp_bond"
             )
         except RuntimeError:
             errors.append(f"Duplicit bond {atom_1} - {atom_2}")

pdbeccdutils/core/prd_reader.py

@@ -18,7 +18,7 @@
 
 import logging
 import os
-from typing import Dict
+from typing import Collection, Dict, Optional
 import rdkit
 from pdbeccdutils.core.component import Component
 from pdbeccdutils.core import ccd_reader
@@ -69,31 +69,47 @@ def read_pdb_cif_file(
 
 
 def read_pdb_components_file(
-    path_to_cif: str, sanitize: bool = True
+    path_to_cif: str, sanitize: bool = True, include: Optional[Collection[str]] = None
 ) -> Dict[str, ccd_reader.CCDReaderResult]:
     """
-    Process multiple compounds stored in the wwPDB CCD
-    `components.cif` file.
+    Process multiple compounds stored in the wwPDB PRD
+    `prdcc-all.cif` file.
 
     Args:
         path_to_cif (str): Path to the `prdcc-all.cif` file with
             multiple ligands in it.
         sanitize (bool): Whether or not the components should be sanitized
             Defaults to True.
+        include (Optional[Collection[str]]): List of PRDs to be parsed. By default it
+            is None and parses all the PRDs. If a collection of PRDs is provided, only
+            those PRDs are parsed. Empty collections parse all PRDs, matching the
+            previous behavior.
 
     Raises:
         ValueError: if the file does not exist.
 
     Returns:
         dict[str, CCDReaderResult]: Internal representation of all
-        the components in the `components.cif` file.
+        the components in the `prdcc-all.cif` file.
     """
     if not os.path.isfile(path_to_cif):
         raise ValueError("File '{}' does not exists".format(path_to_cif))
 
     result_bag = {}
+    doc = cif.read(path_to_cif)
 
-    for block in cif.read(path_to_cif):
+    if not include:
+        blocks = doc
+    else:
+        blocks = []
+        for block_name in dict.fromkeys(include):
+            block = doc.find_block(block_name)
+            if block is None:
+                logging.warning(f"Data block {block_name} not found in {path_to_cif}.")
+                continue
+            blocks.append(block)
+
+    for block in blocks:
         try:
             result_bag[block.name] = _parse_pdb_mmcif(block, sanitize)
         except CCDUtilsError as e:

pdbeccdutils/tests/test_ccd_reader.py

@@ -0,0 +1,93 @@
+import logging
+from pathlib import Path
+
+from gemmi import cif
+from rdkit import Chem
+
+from pdbeccdutils.core import ccd_reader
+from pdbeccdutils.tests.tst_utilities import cif_filename
+
+
+def _components_cif(tmp_path):
+    path = tmp_path / "components.cif"
+    ccds = [Path(cif_filename("00O")), Path(cif_filename("007"))]
+    path.write_text("\n".join(ccd.read_text() for ccd in ccds))
+    return path
+
+
+def test_read_pdb_components_file_parses_all_for_empty_include(tmp_path):
+    path = _components_cif(tmp_path)
+
+    for include in (None, []):
+        result = ccd_reader.read_pdb_components_file(
+            str(path), sanitize=False, include=include
+        )
+
+        assert list(result) == ["00O", "007"]
+
+
+def test_read_pdb_components_file_uses_include_order_and_skips_missing(
+    tmp_path, caplog
+):
+    path = _components_cif(tmp_path)
+
+    with caplog.at_level(logging.WARNING):
+        result = ccd_reader.read_pdb_components_file(
+            str(path), sanitize=False, include=["007", "MISSING", "00O", "007"]
+        )
+
+    assert list(result) == ["007", "00O"]
+    assert "Data block MISSING not found" in caplog.text
+
+
+def test_parse_pdb_bonds_uses_atom_id_lookup():
+    block = cif.read_string(
+        """
+        data_TST
+        loop_
+        _chem_comp_atom.atom_id
+        A
+        B
+        loop_
+        _chem_comp_bond.atom_id_1
+        _chem_comp_bond.atom_id_2
+        _chem_comp_bond.value_order
+        A B SING
+        """
+    ).sole_block()
+    mol = Chem.RWMol()
+    mol.AddAtom(Chem.Atom("C"))
+    mol.AddAtom(Chem.Atom("O"))
+    errors = []
+
+    ccd_reader._parse_pdb_bonds(mol, block, errors)
+
+    assert errors == []
+    assert mol.GetNumBonds() == 1
+    assert mol.GetBondBetweenAtoms(0, 1) is not None
+
+
+def test_parse_pdb_bonds_missing_atom_1_does_not_mask_error():
+    block = cif.read_string(
+        """
+        data_TST
+        loop_
+        _chem_comp_atom.atom_id
+        A
+        loop_
+        _chem_comp_bond.atom_id_1
+        _chem_comp_bond.atom_id_2
+        _chem_comp_bond.value_order
+        ? A SING
+        """
+    ).sole_block()
+    mol = Chem.RWMol()
+    mol.AddAtom(Chem.Atom("C"))
+    errors = []
+
+    ccd_reader._parse_pdb_bonds(mol, block, errors)
+
+    assert errors == [
+        f"Missing atom in 0 entry in _chem_comp_bond"
+    ]
+    assert mol.GetNumBonds() == 0

pdbeccdutils/tests/test_prd_reader.py

@@ -0,0 +1,28 @@
+import logging
+
+from pdbeccdutils.core import prd_reader
+from pdbeccdutils.tests.tst_utilities import prd_cif_filename
+
+
+def test_read_prd_components_file_parses_all_for_empty_include():
+    path = prd_cif_filename("PRDCC_000103")
+
+    for include in (None, []):
+        result = prd_reader.read_pdb_components_file(
+            path, sanitize=False, include=include
+        )
+
+        assert list(result) == ["PRD_000103"]
+
+def test_read_prd_components_file_uses_include_order_and_skips_missing(caplog):
+    path = prd_cif_filename("PRDCC_000103")
+
+    with caplog.at_level(logging.WARNING):
+        result = prd_reader.read_pdb_components_file(
+            path,
+            sanitize=False,
+            include=["PRD_000103", "MISSING", "PRD_000103"],
+        )
+
+    assert list(result) == ["PRD_000103"]
+    assert "Data block MISSING not found" in caplog.text

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pdbeccdutils"
-version = "1.0.3"
+version = "1.0.4"
 description = "Toolkit to parse and process small molecules in wwPDB"
 authors = ["Protein Data Bank in Europe <pdbehelp@ebi.ac.uk>"]
 license = "Apache License 2.0."