Catch errors in sanitize (#59)

* Add test for pred_mol input that dose not sanitize and raises error
* Catch errors in sanitize

Author: maabuu

posebusters/__init__.py

@@ -32,4 +32,4 @@
     "check_volume_overlap",
 ]
 
-__version__ = "0.3.5"
+__version__ = "0.3.6"

posebusters/modules/energy_ratio.py

@@ -11,15 +11,11 @@
 from rdkit.Chem.inchi import InchiReadWriteError, MolFromInchi
 from rdkit.Chem.rdchem import Mol
 from rdkit.Chem.rdDistGeom import EmbedMultipleConfs, ETKDGv3
-from rdkit.Chem.rdForceFieldHelpers import (
-    UFFGetMoleculeForceField,
-    UFFOptimizeMoleculeConfs,
-)
-from rdkit.Chem.rdmolops import AddHs, AssignStereochemistryFrom3D
+from rdkit.Chem.rdForceFieldHelpers import UFFGetMoleculeForceField, UFFOptimizeMoleculeConfs
+from rdkit.Chem.rdmolops import AddHs, AssignStereochemistryFrom3D, SanitizeMol
 
 from ..tools.inchi import get_inchi
 from ..tools.logging import CaptureLogger
-from ..tools.molecules import assert_sanity
 
 logger = logging.getLogger(__name__)
 
@@ -58,7 +54,7 @@ def check_energy_ratio(
 
     try:
         assert mol_pred.GetNumConformers() > 0, "Molecule does not have a conformer."
-        mol_pred = assert_sanity(mol_pred)
+        assert not SanitizeMol(mol_pred, catchErrors=True), "Molecule does not sanitize."
         AddHs(mol_pred, addCoords=True)
     except Exception as e:
         logger.warning(_warning_prefix + "failed because RDKit sanitization failed for molecule: %s", e)

posebusters/modules/rmsd.py

@@ -100,8 +100,8 @@ def robust_rmsd(  # noqa: PLR0913
         RemoveStereochemistry(mol_ref)
 
     if heavy_only:
-        mol_probe = RemoveHs(mol_probe)
-        mol_ref = RemoveHs(mol_ref)
+        mol_probe = RemoveHs(mol_probe, sanitize=False)
+        mol_ref = RemoveHs(mol_ref, sanitize=False)
 
     # combine parameters
     params = dict(symmetrizeConjugatedTerminalGroups=symmetrizeConjugatedTerminalGroups, kabsch=kabsch, **params)
@@ -186,11 +186,15 @@ def intercentroid(
     mol_probe: Mol, mol_ref: Mol, conf_id_probe: int = -1, conf_id_ref: int = -1, heavy_only: bool = True
 ) -> float:
     """Distance between centroids of two molecules."""
-    if heavy_only:
-        mol_probe = RemoveHs(mol_probe)
-        mol_ref = RemoveHs(mol_ref)
-
-    centroid_probe = mol_probe.GetConformer(conf_id_probe).GetPositions().mean(axis=0)
-    centroid_ref = mol_ref.GetConformer(conf_id_ref).GetPositions().mean(axis=0)
 
+    centroid_probe = get_centroid(mol_probe, heavy_only, conf_id_probe)
+    centroid_ref = get_centroid(mol_ref, heavy_only, conf_id_ref)
     return float(np.linalg.norm(centroid_probe - centroid_ref))
+
+
+def get_centroid(mol: Mol, heavy_only: bool = True, conf_id: int = -1) -> np.ndarray:
+    """Get centroid of molecule."""
+    pos = mol.GetConformer(conf_id).GetPositions()
+    if heavy_only:
+        pos = pos[[atom.GetAtomicNum() != 1 for atom in mol.GetAtoms()], :]
+    return pos.mean(axis=0)

posebusters/tools/inchi.py

@@ -11,7 +11,7 @@
 from rdkit.Chem.rdmolops import AssignStereochemistryFrom3D, RemoveHs, RemoveStereochemistry, SanitizeMol
 
 from .logging import CaptureLogger
-from .molecules import add_stereo_hydrogens, assert_sanity, neutralize_atoms, remove_isotopic_info
+from .molecules import add_stereo_hydrogens, neutralize_atoms, remove_isotopic_info
 
 logger = logging.getLogger(__name__)
 
@@ -33,8 +33,12 @@ def get_inchi(mol: Mol, inchi_strict: bool = False) -> str:
 
 def standardize_and_get_inchi(mol: Mol, options: str = "", log_level=None, warnings_as_errors=False) -> str:
     """Return InChI after standardising molecule and inferring stereo from coordinates."""
+
     mol = deepcopy(mol)
-    mol = assert_sanity(mol)
+
+    if flags := SanitizeMol(mol, catchErrors=True):
+        logger.debug("Cannot get InChI because molecule doesn't sanitize due to %s.", flags)
+        return ""
 
     # standardise molecule
     mol = remove_isotopic_info(mol)
@@ -65,8 +69,8 @@ def is_valid_inchi(inchi: str) -> bool:
     """Check that InChI can be parsed and sanitization does not fail."""
     try:
         mol = MolFromInchi(inchi)
-        assert_sanity(mol)
-        assert mol is not None
+        assert mol is not None, "Molecule is None."
+        assert not SanitizeMol(mol, catchErrors=True), "Molecule does not sanitize."
         return True
     except Exception:
         return False

posebusters/tools/molecules.py

@@ -94,13 +94,6 @@ def add_stereo_hydrogens(mol: Mol) -> Mol:
     return mol
 
 
-def assert_sanity(mol: Mol) -> Mol:
-    """Check that RDKit sanitization does not fail."""
-    flags = SanitizeMol(mol)
-    assert flags == 0, f"Sanitization failed with flags {flags}"
-    return mol
-
-
 def remove_isotopic_info(mol: Mol) -> Mol:
     """Remove isotopic information from molecule."""
     for atom in mol.GetAtoms():

tests/conftest/sanity_error/pred.sdf

@@ -0,0 +1,32 @@
+/pred.pdb
+ OpenBabel03212512233D
+
+ 13 12  0  0  1  0  0  0  0  0999 V2000
+   20.0370   46.5970   12.1140 N   0  0  0  0  0  0  0  0  0  0  0  0
+   18.9380   46.2330   16.2850 P   0  0  2  0  0  0  0  0  0  0  0  0
+   20.0790   47.1960   18.3360 C   0  0  0  0  0  0  0  0  0  0  0  0
+   19.6880   48.4400   19.0620 C   0  0  1  0  0  0  0  0  0  0  0  0
+   19.7350   48.1890   20.5450 C   0  0  0  0  0  0  0  0  0  0  0  0
+   18.8130   45.4360   13.8410 C   0  0  0  0  0  0  0  0  0  0  0  0
+   19.7710   47.3330   16.9860 O   0  0  0  0  0  0  0  0  0  0  0  0
+   19.6730   45.2860   12.6310 C   0  0  0  0  0  0  0  0  0  0  0  0
+   19.5130   46.1760   14.8290 O   0  0  0  0  0  0  0  0  0  0  0  0
+   17.4840   46.6470   16.2100 O   0  5  0  0  0  0  0  0  0  0  0  0
+   19.1520   44.8900   16.9230 O   0  0  0  0  0  0  0  0  0  0  0  0
+   20.5340   49.4930   18.6810 O   0  0  0  0  0  0  0  0  0  0  0  0
+   19.2020   49.3120   21.2340 O   0  0  0  0  0  0  0  0  0  0  0  0
+  1  8  1  0  0  0  0
+  2  7  1  0  0  0  0
+  2  9  1  0  0  0  0
+  2 10  2  6  0  0  0
+  2 11  1  0  0  0  0
+  3  4  1  0  0  0  0
+  3  7  1  0  0  0  0
+  4  5  1  0  0  0  0
+  4 12  1  1  0  0  0
+  5 13  1  0  0  0  0
+  6  8  1  0  0  0  0
+  6  9  1  0  0  0  0
+M  CHG  1  10  -1
+M  END
+$$$$