Update for RDKit 2024.9 (#51)

Dependencies
- Add support for Python 3.13
- Drop support for Python 3.8
- Require RDKit version 2024.9 as a minimum

Energy ratio module
- Update energy ratio tests to be compatible with InChI 1.07

Development
- Add GitHub stale action

Author: maabuu

.github/workflows/inactivity.yml

@@ -0,0 +1,45 @@
+name: Inactivity
+
+on:
+    workflow_dispatch:
+    # push:
+    #   branches: [main]
+    # pull_request:
+    schedule:
+      #        ┌───────────── minute (0 - 59)
+      #        │  ┌───────────── hour (0 - 23)
+      #        │  │ ┌───────────── day of the month (1 - 31)
+      #        │  │ │ ┌───────────── month (1 - 12 or JAN-DEC)
+      #        │  │ │ │ ┌───────────── day of the week (0 - 6 or SUN-SAT)
+      #        │  │ │ │ │
+      #        │  │ │ │ │
+      #        │  │ │ │ │
+      #        *  * * * *
+      - cron: '47 2 3 * *'
+
+permissions:
+  issues: write
+  pull-requests: write
+
+jobs:
+  close-issues:
+    runs-on: ubuntu-latest
+
+    steps:
+
+    - name: Close inactive issues, PRs, and stale issues
+      uses: actions/stale@v3
+      with:
+        repo-token: "${{ secrets.GITHUB_TOKEN }}"
+        days-before-stale: 60
+        days-before-close: 7
+        # Issues
+        stale-issue-message: "This issue has been automatically closed due to inactivity. Please feel free to reopen it if you still have this problem."
+        close-issue-message: "This issue was closed because it has been inactive for 14 days since being marked as stale."
+        stale-issue-label: "stale"
+        exempt-issue-labels: "in-progress"
+        # PRs
+        stale-pr-label: "stale"
+        exempt-all-pr-assignees: true
+        stale-pr-message: "This PR was marked as stale because it has been open for 60 days with no activity."
+        close-pr-message: "This PR was closed because it has been inactive for 14 days since being marked as stale."

.github/workflows/testing.yml

@@ -14,7 +14,7 @@ jobs:
       fail-fast: true
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        python: ["3.9", "3.10", "3.11", "3.12", "3.13"]
 
     steps:
 

posebusters/cli.py

@@ -1,11 +1,13 @@
 """Command line interface for PoseBusters."""
+
 from __future__ import annotations
 
 import argparse
 import logging
 import sys
+from collections.abc import Iterable
 from pathlib import Path
-from typing import Any, Iterable, TextIO
+from typing import Any, TextIO
 
 import pandas as pd
 from rdkit.Chem.rdchem import Mol

posebusters/modules/distance_geometry.py

@@ -1,9 +1,11 @@
 """Module to check bond lengths, bond angles, and internal clash of ligand conformations."""
+
 from __future__ import annotations
 
+from collections.abc import Iterable
 from copy import deepcopy
 from logging import getLogger
-from typing import Any, Iterable
+from typing import Any
 
 import numpy as np
 import pandas as pd

posebusters/modules/energy_ratio.py

@@ -4,9 +4,9 @@
 
 import logging
 from copy import deepcopy
-from functools import lru_cache
+from functools import cache
+from math import isfinite
 
-import numpy as np
 from rdkit import ForceField  # noqa: F401
 from rdkit.Chem.inchi import InchiReadWriteError, MolFromInchi
 from rdkit.Chem.rdchem import Mol
@@ -23,13 +23,14 @@
 
 logger = logging.getLogger(__name__)
 
+
 _warning_prefix = "WARNING: Energy ratio module "
 _empty_results = {
     "results": {
-        "ensemble_avg_energy": np.nan,
-        "mol_pred_energy": np.nan,
-        "energy_ratio": np.nan,
-        "energy_ratio_passes": np.nan,
+        "ensemble_avg_energy": float("nan"),
+        "mol_pred_energy": float("nan"),
+        "energy_ratio": float("nan"),
+        "energy_ratio_passes": float("nan"),
     }
 }
 
@@ -39,7 +40,8 @@ def check_energy_ratio(
     threshold_energy_ratio: float = 7.0,
     ensemble_number_conformations: int = 100,
     inchi_strict: bool = False,
-):
+    epsilon=1e-10,
+) -> dict[str, dict[str, float | bool]]:
     """Check whether the energy of the docked ligand is within user defined range.
 
     Args:
@@ -72,40 +74,65 @@ def check_energy_ratio(
         return _empty_results
 
     try:
-        conf_energy = get_conf_energy(mol_pred)
+        observed_energy = get_conf_energy(mol_pred)
     except Exception as e:
         logger.warning(_warning_prefix + "failed to calculate conformation energy for %s: %s", inchi, e)
-        conf_energy = np.nan
+        observed_energy = float("nan")
 
     try:
-        avg_energy = float(get_average_energy(inchi, ensemble_number_conformations))
+        energies = get_energies(inchi, ensemble_number_conformations)
+        mean_energy = sum(energies) / len(energies)
+        std_energy = sum((energy - mean_energy) ** 2 for energy in energies) / (len(energies) - 1)
+        std_energy = max(epsilon, std_energy)  # clipping
     except Exception as e:
         logger.warning(_warning_prefix + "failed to calculate ensemble conformation energy for %s: %s", inchi, e)
-        avg_energy = np.nan
+        mean_energy = float("nan")
+        std_energy = float("nan")
 
-    if avg_energy == 0:
+    if mean_energy == 0:
         logger.warning(_warning_prefix + "calculated average energy of molecule 0 for %s", inchi)
-        avg_energy = np.nan
+        mean_energy = epsilon  # clipping
+
+    # simple ratio
+    ratio = observed_energy / mean_energy
+    ratio_passes = ratio <= threshold_energy_ratio if isfinite(ratio) else float("nan")
 
-    pred_factor = conf_energy / avg_energy
-    ratio_passes = pred_factor <= threshold_energy_ratio
+    # ratio after subtracting mean
+    deviation = observed_energy - mean_energy
+    relative_deviation = deviation / mean_energy
+    relative_deviation_passes = (
+        relative_deviation <= threshold_energy_ratio if isfinite(relative_deviation) else float("nan")
+    )
+
+    # standard score (ratio after subtracting by population mean and dividing by population std)
+    z_value = (observed_energy - mean_energy) / std_energy
+    z_value_passes = z_value <= threshold_energy_ratio if isfinite(z_value) else float("nan")
 
     results = {
-        "ensemble_avg_energy": avg_energy,
-        "mol_pred_energy": conf_energy,
-        "energy_ratio": pred_factor,
+        "ensemble_avg_energy": mean_energy,
+        "mol_pred_energy": observed_energy,
+        "energy_ratio": ratio,
+        "relative_deviation": relative_deviation,
+        "z_value": z_value,
         "energy_ratio_passes": ratio_passes,
+        "relative_deviation_passes": relative_deviation_passes,
+        "z_value_passes": z_value_passes,
     }
     return {"results": results}
 
 
-@lru_cache(maxsize=None)
 def get_average_energy(inchi: str, n_confs: int = 50, num_threads: int = 0) -> Mol:
     """Get average energy of an ensemble of molecule conformations."""
+    energies = get_energies(inchi, n_confs, num_threads)
+    return sum(energies) / len(energies)
+
+
+@cache
+def get_energies(inchi: str, n_confs: int = 50, num_threads: int = 0) -> list[float]:
+    """Get energies of an ensemble of molecule conformations."""
     with CaptureLogger():
         mol = MolFromInchi(inchi)
-    energies = new_conformation(mol, n_confs, num_threads)["energies"]
-    return sum(energies) / len(energies)
+    return new_conformation(mol, n_confs, num_threads)["energies"]
 
 
 def new_conformation(mol: Mol, n_confs: int = 1, num_threads: int = 0, energy_minimization=True) -> Mol:

posebusters/modules/flatness.py

@@ -1,8 +1,10 @@
 """Module to check flatness of ligand substructures."""
+
 from __future__ import annotations
 
+from collections.abc import Iterable
 from copy import deepcopy
-from typing import Any, Iterable
+from typing import Any
 
 import numpy as np
 from numpy import ndarray as Array

posebusters/modules/loading.py

@@ -1,4 +1,5 @@
 """Module to check loading of protein and ligand files."""
+
 from __future__ import annotations
 
 from typing import Any

posebusters/modules/rmsd.py

@@ -1,4 +1,5 @@
 """Module to check RMSD between docked and crystal ligand."""
+
 from __future__ import annotations
 
 import logging

posebusters/modules/volume_overlap.py

@@ -1,4 +1,5 @@
 """Module to check volume overlap between docked ligand and protein."""
+
 from __future__ import annotations
 
 import logging

posebusters/posebusters.py

@@ -5,9 +5,10 @@
 import inspect
 import logging
 from collections import defaultdict
+from collections.abc import Generator, Iterable
 from functools import partial
 from pathlib import Path
-from typing import Any, Callable, Generator, Iterable
+from typing import Any, Callable
 
 import pandas as pd
 from rdkit.Chem.rdchem import Mol

posebusters/tools/loading.py

@@ -3,8 +3,8 @@
 from __future__ import annotations
 
 import logging
+from collections.abc import Generator
 from pathlib import Path
-from typing import Generator
 
 from rdkit.Chem.AllChem import AssignBondOrdersFromTemplate
 from rdkit.Chem.rdchem import Mol

posebusters/tools/molecules.py

@@ -2,9 +2,9 @@
 
 from __future__ import annotations
 
+from collections.abc import Iterable
 from copy import deepcopy
 from logging import getLogger
-from typing import Iterable
 
 import numpy as np
 from rdkit import RDLogger

posebusters/tools/parallel.py

@@ -2,8 +2,9 @@
 
 from __future__ import annotations
 
+from collections.abc import Iterable
 from multiprocessing import Pool
-from typing import Callable, Iterable
+from typing import Callable
 
 from tqdm import tqdm
 

posebusters/tools/protein.py

@@ -1,8 +1,9 @@
 """Protein related functions."""
+
 from __future__ import annotations
 
 import logging
-from typing import Iterable
+from collections.abc import Iterable
 
 from rdkit.Chem.rdchem import Atom, Mol
 

pyproject.toml

@@ -16,16 +16,16 @@ classifiers = [
     "Intended Audience :: Science/Research",
     "License :: OSI Approved :: BSD License",
     "Programming Language :: Python :: 3 :: Only",
-    "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
     "Operating System :: OS Independent",
 ]
-requires-python = "~=3.8"
+requires-python = "~=3.9"
 dynamic = ["version", "description"]
-dependencies = ['rdkit >= 2020.09', 'pandas', 'numpy', 'pyyaml']
+dependencies = ['rdkit >= 2024.9', 'pandas', 'numpy', 'pyyaml']
 
 [project.optional-dependencies]
 dev = ["ruff"]

tests/test_modules/test_energy_ratio.py

@@ -1,42 +1,63 @@
 from __future__ import annotations
 
+import math
+
 import numpy as np
 
 from posebusters.modules.energy_ratio import check_energy_ratio
 
 
 def test_check_energy_ratio(mol_pm2):
+    # molecule has valid conformation
     out = check_energy_ratio(mol_pm2)
+    assert math.isfinite(out["results"]["mol_pred_energy"])
+    assert math.isfinite(out["results"]["ensemble_avg_energy"])
+    assert math.isfinite(out["results"]["energy_ratio"])
     assert out["results"]["energy_ratio_passes"] is True
 
 
 def test_check_energy_ratio_14gs_0(mol_pred_14gs_gen0):
-    # nan because molecule cannot be converted to valid InChI
+    # molecule is chemically insane, but conformation looks alright
     out = check_energy_ratio(mol_pred_14gs_gen0)
+    assert math.isfinite(out["results"]["mol_pred_energy"])
+    assert math.isfinite(out["results"]["ensemble_avg_energy"])
+    assert math.isfinite(out["results"]["energy_ratio"])
     assert out["results"]["energy_ratio_passes"] is True
 
 
 def test_check_energy_ratio_1afs_87(mol_pred_1afs_gen87):
     # nan because molecule cannot be converted to valid InChI
     out = check_energy_ratio(mol_pred_1afs_gen87)
-    assert np.isnan(out["results"]["energy_ratio_passes"])
+    assert math.isnan(out["results"]["mol_pred_energy"])
+    assert math.isnan(out["results"]["ensemble_avg_energy"])
+    assert math.isnan(out["results"]["energy_ratio"])
+    assert math.isnan(out["results"]["energy_ratio_passes"])
 
 
 def test_check_energy_ratio_1afs_94(mol_pred_1afs_gen94):
-    # nan because molecule cannot be converted to valid InChI
+    # molecule is chemically insane, but conformation looks alright
     out = check_energy_ratio(mol_pred_1afs_gen94)
-    assert out["results"]["energy_ratio_passes"] is True
+    assert math.isfinite(out["results"]["mol_pred_energy"])
+    assert math.isnan(out["results"]["ensemble_avg_energy"])
+    assert math.isnan(out["results"]["energy_ratio"])
+    assert math.isnan(out["results"]["energy_ratio_passes"])
 
 
 def test_check_energy_ratio_1jn2_3(mol_pred_1jn2_gen3):
     # nan because molecule cannot be converted to valid InChI
     out = check_energy_ratio(mol_pred_1jn2_gen3)
-    assert np.isnan(out["results"]["energy_ratio_passes"])
+    assert math.isnan(out["results"]["mol_pred_energy"])
+    assert math.isnan(out["results"]["ensemble_avg_energy"])
+    assert math.isnan(out["results"]["energy_ratio"])
+    assert math.isnan(out["results"]["energy_ratio_passes"])
 
 
 def test_check_energy_ratio_1jn2_62(mol_pred_1jn2_gen62):
-    # nan because molecule cannot be converted to valid InChI
+    # molecule is chemically insane, but conformation looks alright
     out = check_energy_ratio(mol_pred_1jn2_gen62)
+    assert math.isfinite(out["results"]["mol_pred_energy"])
+    assert math.isfinite(out["results"]["ensemble_avg_energy"])
+    assert math.isfinite(out["results"]["energy_ratio"])
     assert out["results"]["energy_ratio_passes"] is True
 
 
@@ -49,5 +70,9 @@ def test_check_energy_ratio_approximate_consistency(mol_1a30_clash_2, mol_1a30_c
 
 
 def test_check_energy_ratio_disconnected_atoms(mol_disconnnected_atoms):
+    # no bonds, just disconnected atoms -> energy is 0
     out = check_energy_ratio(mol_disconnnected_atoms)
-    assert out["results"]["energy_ratio_passes"] is False
+    assert math.isfinite(out["results"]["mol_pred_energy"])
+    assert math.isfinite(out["results"]["ensemble_avg_energy"])
+    assert math.isfinite(out["results"]["energy_ratio"])
+    assert out["results"]["energy_ratio_passes"] is True