From 08483aba6338b38ca4a19db38f9dc2b04c52f2ee Mon Sep 17 00:00:00 2001 From: ccdc-github-admin <62256659+ccdc-github-admin@users.noreply.github.com> Date: Sat, 20 Jun 2026 00:19:06 +0100 Subject: [PATCH 01/15] ci(NO_JIRA): Updated .github/workflows/lint-python.yml --- .github/workflows/lint-python.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/lint-python.yml b/.github/workflows/lint-python.yml index 5b430e5..245010d 100644 --- a/.github/workflows/lint-python.yml +++ b/.github/workflows/lint-python.yml @@ -15,7 +15,7 @@ jobs: python_version: ["3.11"] steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 with: fetch-depth: 0 From e0533c3c06c211e4de034a0758755dd702100087 Mon Sep 17 00:00:00 2001 From: ccdc-github-admin <62256659+ccdc-github-admin@users.noreply.github.com> Date: Sat, 20 Jun 2026 00:19:11 +0100 Subject: [PATCH 02/15] ci(NO_JIRA): Updated .github/workflows/lint-markdown.yml --- .github/workflows/lint-markdown.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/lint-markdown.yml b/.github/workflows/lint-markdown.yml index 102fc98..8c393eb 100644 --- a/.github/workflows/lint-markdown.yml +++ b/.github/workflows/lint-markdown.yml @@ -12,7 +12,7 @@ jobs: changed-files: ${{ steps.file-changes.outputs.changed-files }} steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 with: fetch-depth: 0 From 6c20006bc54103a270572e4b6ecbaf223f99fb40 Mon Sep 17 00:00:00 2001 From: ccdc-github-admin <62256659+ccdc-github-admin@users.noreply.github.com> Date: Sat, 20 Jun 2026 00:19:15 +0100 Subject: [PATCH 03/15] ci(NO_JIRA): Updated .github/workflows/common_ccdc_status_checks.yml --- .github/workflows/common_ccdc_status_checks.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/common_ccdc_status_checks.yml b/.github/workflows/common_ccdc_status_checks.yml index 406af0f..7ec5ce2 100644 --- a/.github/workflows/common_ccdc_status_checks.yml +++ b/.github/workflows/common_ccdc_status_checks.yml @@ -10,7 +10,7 @@ jobs: ccdc-commit-hooks-on-pull-request-check: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 with: ref: ${{ github.head_ref }} fetch-depth: 0 From 98a694d138985d7209c7afad9a7ca8261dda5310 Mon Sep 17 00:00:00 2001 From: ccdc-github-admin <62256659+ccdc-github-admin@users.noreply.github.com> Date: Sat, 20 Jun 2026 00:19:19 +0100 Subject: [PATCH 04/15] ci(NO_JIRA): Updated .github/workflows/lint-dockerfiles.yml --- .github/workflows/lint-dockerfiles.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/lint-dockerfiles.yml b/.github/workflows/lint-dockerfiles.yml index 8689b5d..9937646 100644 --- a/.github/workflows/lint-dockerfiles.yml +++ b/.github/workflows/lint-dockerfiles.yml @@ -12,7 +12,7 @@ jobs: changed-files: ${{ steps.file-changes.outputs.changed-files }} steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 with: fetch-depth: 0 From fb7b6551f97510609883522f6389616ffcf31522 Mon Sep 17 00:00:00 2001 From: ccdc-github-admin <62256659+ccdc-github-admin@users.noreply.github.com> Date: Sat, 20 Jun 2026 00:19:24 +0100 Subject: [PATCH 05/15] ci(NO_JIRA): Updated .github/workflows/lint-yaml.yml --- .github/workflows/lint-yaml.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/lint-yaml.yml b/.github/workflows/lint-yaml.yml index 745d6d7..6436973 100644 --- a/.github/workflows/lint-yaml.yml +++ b/.github/workflows/lint-yaml.yml @@ -13,7 +13,7 @@ jobs: changed-files: ${{ steps.file-changes.outputs.changed-files }} steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 with: fetch-depth: 0 From fc81482a4719391a859f3244b61b37281c06390f Mon Sep 17 00:00:00 2001 From: ccdc-github-admin <62256659+ccdc-github-admin@users.noreply.github.com> Date: Wed, 24 Jun 2026 23:06:39 +0100 Subject: [PATCH 06/15] ci(NO_JIRA): Updated .yamllint.yml --- .yamllint.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.yamllint.yml b/.yamllint.yml index a2c76e0..bc4e534 100644 --- a/.yamllint.yml +++ b/.yamllint.yml @@ -5,11 +5,11 @@ # Project Homepage: https://github.com/adrienverge/yamllint # # Overriding rules in files: -# https://yamllint.readthedocs.io/en/latest/disable_with_comments.html +# http://yamllint.readthedocs.io/en/latest/disable_with_comments.html --- extends: default -# Rules documentation: https://yamllint.readthedocs.io/en/latest/rules.html +# Rules documentation: http://yamllint.readthedocs.io/en/latest/rules.html rules: document-start: disable line-length: From a6f9c79683c681d07adca20476a3a39e30b6f824 Mon Sep 17 00:00:00 2001 From: SchrodingersCattt Date: Mon, 29 Jun 2026 09:33:45 +0000 Subject: [PATCH 07/15] Add disorder-aware CIF export demo --- .../disorder_workflow_demo/ReadMe.md | 46 +++ .../demo_csd_disorder_workflow.py | 280 ++++++++++++++++++ 2 files changed, 326 insertions(+) create mode 100644 scripts/materials/disorder_workflow_demo/ReadMe.md create mode 100644 scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py diff --git a/scripts/materials/disorder_workflow_demo/ReadMe.md b/scripts/materials/disorder_workflow_demo/ReadMe.md new file mode 100644 index 0000000..ce42aa1 --- /dev/null +++ b/scripts/materials/disorder_workflow_demo/ReadMe.md @@ -0,0 +1,46 @@ +# Disorder Workflow Demo + +## Summary + +Exports a CSD entry to a manually assembled CIF that preserves atom-site occupancy, +disorder assembly/group annotations, anisotropic displacement parameters, and bond +connectivity. + +This example demonstrates a workflow that is not covered by the built-in CIF writers: +reading a disordered entry from the CSD Python API and reconstructing a CIF with +explicit `_atom_site_occupancy`, `_atom_site_disorder_assembly`, and +`_atom_site_disorder_group` fields. + +Example output: + +```text +Licence requirement: CSD-Materials or another licence tier with CSD Python API database access +Exported ABACIR -> output/ABACIR.cif (atoms=148, partial_occ=16, bonds=154) +Exported ABABUB -> output/ABABUB.cif (atoms=84, partial_occ=8, bonds=88) +``` + +## Requirements + +- CSD Python API with database access +- Python 3.10+ + +## Licensing Requirements + +CSD-Materials or another licence tier that provides access to the CSD Python API and CSD database. + +## Instructions on Running + +```bash +python demo_csd_disorder_workflow.py --refcodes ABACIR ABABUB --output-dir ./output +``` + +Arguments: + +- `--refcodes`: one or more CSD refcodes +- `--output-dir`: destination directory for exported CIF files + +## Author + +Ming-Yu Guo (2026-06-29) + +> For feedback or to report any issues please contact [support@ccdc.cam.ac.uk](mailto:support@ccdc.cam.ac.uk) diff --git a/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py b/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py new file mode 100644 index 0000000..8c6922d --- /dev/null +++ b/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py @@ -0,0 +1,280 @@ +#!/usr/bin/env python +# +# This script can be used for any purpose without limitation subject to the +# conditions at https://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx +# +# This permission notice and the following statement of attribution must be +# included in all copies or substantial portions of this script. +# +# 2026-06-29: created by Ming-Yu Guo +# + +"""Demo: export full disorder-aware CIF data from the CSD Python API. + +This script demonstrates a complete workflow for extracting crystal structures +from the Cambridge Structural Database (CSD) with atom-site occupancy, +disorder-group metadata, anisotropic displacement parameters, and bond +connectivity preserved in a manually assembled CIF file. + +Why this script exists: + The built-in CSD Python API CIF writers are convenient, but for some + disordered entries they do not expose a simple path to preserve the full + occupancy/disorder atom-site metadata needed for downstream enumeration or + custom reconstruction workflows. + +Pipeline: + 1. Read an entry from the CSD with ``EntryReader``. + 2. Access ``crystal.disordered_molecule`` when available. + 3. Collect atom-site occupancy, disorder assembly/group, ADP, and bonds. + 4. Write a CIF with explicit ``_atom_site_occupancy`` and disorder tags. + +Examples: + python demo_csd_disorder_workflow.py --refcodes ABACIR ABABUB --output-dir ./output +""" + +from __future__ import annotations + +import argparse +from pathlib import Path +from typing import Any + + +LICENSE_REQUIREMENT = "CSD-Materials or another licence tier with CSD Python API database access" + + +def _format_float(value: Any, precision: int = 4) -> str: + if value is None: + return "?" + return f"{float(value):.{precision}f}" + + +def _format_esd(value: Any, esd: Any, precision: int = 4) -> str: + if value is None: + return "?" + base = f"{float(value):.{precision}f}" + if esd is None: + return base + try: + if float(esd) > 0: + return f"{base}({esd})" + except (TypeError, ValueError): + pass + return base + + +def export_full_cif_from_csd(refcode: str, output_path: Path) -> dict[str, Any]: + """Export a CSD structure to CIF with occupancy and disorder metadata.""" + from ccdc.io import EntryReader + + with EntryReader("CSD") as reader: + entry = reader.entry(refcode) + crystal = reader.crystal(refcode) + + molecule = crystal.disordered_molecule or crystal.molecule + + disorder_map: dict[str, tuple[int | None, str | None]] = {} + disorder = getattr(crystal, "disorder", None) + if disorder is not None: + try: + for assembly in disorder.assemblies: + assembly_id = str(getattr(assembly, "id", "")) + for group in assembly.groups: + group_id = int(getattr(group, "id", 0)) + for atom in group.atoms: + disorder_map[str(atom.label)] = (group_id, assembly_id) + except Exception: + pass + + atom_data: list[dict[str, Any]] = [] + for atom in molecule.atoms: + label = str(atom.label) if atom.label else "?" + symbol = str(atom.atomic_symbol) if atom.atomic_symbol else "?" + fractional = atom.fractional_coordinates + occupancy = 1.0 + try: + if atom.occupancy is not None: + occupancy = float(atom.occupancy) + except Exception: + pass + + displacement_type = "?" + u_iso = None + aniso = None + try: + displacement = atom.displacement_parameters + if displacement is not None: + u_iso = float(displacement.isotropic_equivalent) + if displacement.type == "Anisotropic": + displacement_type = "Uani" + values = displacement.values + uncertainties = displacement.uncertainties + aniso = { + "u11": values[0][0], "u22": values[1][1], "u33": values[2][2], + "u12": values[0][1], "u13": values[0][2], "u23": values[1][2], + "e11": uncertainties[0][0], "e22": uncertainties[1][1], "e33": uncertainties[2][2], + "e12": uncertainties[0][1], "e13": uncertainties[0][2], "e23": uncertainties[1][2], + } + else: + displacement_type = "Uiso" + except Exception: + pass + + disorder_group, disorder_assembly = disorder_map.get(label, (None, None)) + atom_data.append({ + "label": label, + "symbol": symbol, + "fx": float(fractional.x) if fractional else None, + "fy": float(fractional.y) if fractional else None, + "fz": float(fractional.z) if fractional else None, + "occupancy": occupancy, + "u_iso": u_iso, + "displacement_type": displacement_type, + "disorder_group": disorder_group, + "disorder_assembly": disorder_assembly, + "aniso": aniso, + }) + + bonds: list[tuple[str, str, float, str]] = [] + try: + for bond in molecule.bonds: + atom_1, atom_2 = bond.atoms + bonds.append((str(atom_1.label), str(atom_2.label), float(bond.length), str(bond.bond_type))) + except Exception: + pass + + cell_a, cell_b, cell_c = crystal.cell_lengths + alpha, beta, gamma = crystal.cell_angles + disorder_present = any(atom["disorder_group"] is not None for atom in atom_data) + + lines = [ + f"data_{refcode}", + "_audit_creation_method 'CSD Python API disorder workflow demo'", + f"_chemical_name_common '{(entry.chemical_name or '?').replace(chr(39), chr(39) + chr(39))}'", + f"_chemical_formula_sum '{entry.formula or '?'}'", + f"_cell_length_a {_format_float(cell_a)}", + f"_cell_length_b {_format_float(cell_b)}", + f"_cell_length_c {_format_float(cell_c)}", + f"_cell_angle_alpha {_format_float(alpha, 2)}", + f"_cell_angle_beta {_format_float(beta, 2)}", + f"_cell_angle_gamma {_format_float(gamma, 2)}", + f"_cell_volume {_format_float(crystal.cell_volume, 2)}", + f"_cell_formula_units_Z {int(crystal.z_value) if crystal.z_value else '?'}", + f"_symmetry_space_group_name_H-M '{crystal.spacegroup_symbol or '?'}'", + "", + ] + + if crystal.symmetry_operators: + lines.extend(["loop_", "_symmetry_equiv_pos_as_xyz"]) + for operator in crystal.symmetry_operators: + lines.append(f" '{operator}'") + lines.append("") + + lines.extend([ + "loop_", + "_atom_site_label", + "_atom_site_type_symbol", + "_atom_site_fract_x", + "_atom_site_fract_y", + "_atom_site_fract_z", + "_atom_site_occupancy", + "_atom_site_U_iso_or_equiv", + "_atom_site_thermal_displace_type", + ]) + if disorder_present: + lines.extend([ + "_atom_site_disorder_assembly", + "_atom_site_disorder_group", + ]) + + for atom in atom_data: + if atom["fx"] is None: + continue + row = [ + f"{atom['label']:<8s}", + f"{atom['symbol']:<4s}", + f"{atom['fx']:12.6f}", + f"{atom['fy']:12.6f}", + f"{atom['fz']:12.6f}", + f"{atom['occupancy']:8.4f}", + f"{_format_float(atom['u_iso']):>10s}", + f"{atom['displacement_type']:>5s}", + ] + if disorder_present: + row.append(f"{str(atom['disorder_assembly']) if atom['disorder_assembly'] is not None else '.':>4s}") + row.append(f"{str(atom['disorder_group']) if atom['disorder_group'] is not None else '.':>4s}") + lines.append(" " + " ".join(row)) + lines.append("") + + anisotropic_atoms = [atom for atom in atom_data if atom["aniso"]] + if anisotropic_atoms: + lines.extend([ + "loop_", + "_atom_site_aniso_label", + "_atom_site_aniso_U_11", + "_atom_site_aniso_U_22", + "_atom_site_aniso_U_33", + "_atom_site_aniso_U_12", + "_atom_site_aniso_U_13", + "_atom_site_aniso_U_23", + ]) + for atom in anisotropic_atoms: + u = atom["aniso"] + lines.append( + f" {atom['label']:<8s}" + f" {_format_esd(u['u11'], u['e11']):>12s}" + f" {_format_esd(u['u22'], u['e22']):>12s}" + f" {_format_esd(u['u33'], u['e33']):>12s}" + f" {_format_esd(u['u12'], u['e12']):>12s}" + f" {_format_esd(u['u13'], u['e13']):>12s}" + f" {_format_esd(u['u23'], u['e23']):>12s}" + ) + lines.append("") + + if bonds: + lines.extend([ + "loop_", + "_geom_bond_atom_site_label_1", + "_geom_bond_atom_site_label_2", + "_geom_bond_distance", + "_ccdc_geom_bond_type", + ]) + for atom_1, atom_2, distance, bond_type in bonds: + lines.append(f" {atom_1:<8s} {atom_2:<8s} {distance:8.4f} {bond_type}") + lines.append("") + + lines.append("#END") + output_path.parent.mkdir(parents=True, exist_ok=True) + output_path.write_text("\n".join(lines), encoding="utf-8") + + return { + "refcode": refcode, + "output": str(output_path), + "n_atoms": len(atom_data), + "n_partial_occupancy": sum(1 for atom in atom_data if atom["occupancy"] < 0.999), + "n_bonds": len(bonds), + } + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--refcodes", nargs="+", required=True, help="CSD refcodes to export") + parser.add_argument("--output-dir", default="output", help="Directory for exported CIF files") + return parser.parse_args() + + +def main() -> None: + args = parse_args() + output_dir = Path(args.output_dir) + + print(f"Licence requirement: {LICENSE_REQUIREMENT}") + for refcode in args.refcodes: + output_path = output_dir / f"{refcode}.cif" + summary = export_full_cif_from_csd(refcode, output_path) + print( + f"Exported {summary['refcode']} -> {summary['output']} " + f"(atoms={summary['n_atoms']}, partial_occ={summary['n_partial_occupancy']}, bonds={summary['n_bonds']})" + ) + + +if __name__ == "__main__": + main() From 0af27567b215cf2653ae6fbdfb28670dba62b8db Mon Sep 17 00:00:00 2001 From: SchrodingersCattt Date: Mon, 29 Jun 2026 09:56:07 +0000 Subject: [PATCH 08/15] Refactor disorder export demo for Sonar --- .../demo_csd_disorder_workflow.py | 260 +++++++++++------- 1 file changed, 160 insertions(+), 100 deletions(-) diff --git a/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py b/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py index 8c6922d..74a9a8b 100644 --- a/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py +++ b/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py @@ -42,6 +42,10 @@ LICENSE_REQUIREMENT = "CSD-Materials or another licence tier with CSD Python API database access" +def _escape_cif_string(value: str | None) -> str: + return (value or "?").replace("'", "''") + + def _format_float(value: Any, precision: int = 4) -> str: if value is None: return "?" @@ -62,63 +66,64 @@ def _format_esd(value: Any, esd: Any, precision: int = 4) -> str: return base -def export_full_cif_from_csd(refcode: str, output_path: Path) -> dict[str, Any]: - """Export a CSD structure to CIF with occupancy and disorder metadata.""" - from ccdc.io import EntryReader +def _build_disorder_map(crystal: Any) -> dict[str, tuple[int | None, str | None]]: + disorder_map: dict[str, tuple[int | None, str | None]] = {} + disorder = getattr(crystal, "disorder", None) + if disorder is None: + return disorder_map - with EntryReader("CSD") as reader: - entry = reader.entry(refcode) - crystal = reader.crystal(refcode) + try: + for assembly in disorder.assemblies: + assembly_id = str(getattr(assembly, "id", "")) + for group in assembly.groups: + group_id = int(getattr(group, "id", 0)) + for atom in group.atoms: + disorder_map[str(atom.label)] = (group_id, assembly_id) + except Exception: + return disorder_map + return disorder_map - molecule = crystal.disordered_molecule or crystal.molecule - disorder_map: dict[str, tuple[int | None, str | None]] = {} - disorder = getattr(crystal, "disorder", None) - if disorder is not None: - try: - for assembly in disorder.assemblies: - assembly_id = str(getattr(assembly, "id", "")) - for group in assembly.groups: - group_id = int(getattr(group, "id", 0)) - for atom in group.atoms: - disorder_map[str(atom.label)] = (group_id, assembly_id) - except Exception: - pass +def _extract_occupancy(atom: Any) -> float: + try: + if atom.occupancy is not None: + return float(atom.occupancy) + except Exception: + pass + return 1.0 + +def _extract_displacement(atom: Any) -> tuple[str, float | None, dict[str, Any] | None]: + displacement_type = "?" + u_iso = None + aniso = None + try: + displacement = atom.displacement_parameters + if displacement is None: + return displacement_type, u_iso, aniso + u_iso = float(displacement.isotropic_equivalent) + if displacement.type != "Anisotropic": + return "Uiso", u_iso, aniso + values = displacement.values + uncertainties = displacement.uncertainties + aniso = { + "u11": values[0][0], "u22": values[1][1], "u33": values[2][2], + "u12": values[0][1], "u13": values[0][2], "u23": values[1][2], + "e11": uncertainties[0][0], "e22": uncertainties[1][1], "e33": uncertainties[2][2], + "e12": uncertainties[0][1], "e13": uncertainties[0][2], "e23": uncertainties[1][2], + } + return "Uani", u_iso, aniso + except Exception: + return displacement_type, u_iso, aniso + + +def _extract_atom_data(molecule: Any, disorder_map: dict[str, tuple[int | None, str | None]]) -> list[dict[str, Any]]: atom_data: list[dict[str, Any]] = [] for atom in molecule.atoms: label = str(atom.label) if atom.label else "?" symbol = str(atom.atomic_symbol) if atom.atomic_symbol else "?" fractional = atom.fractional_coordinates - occupancy = 1.0 - try: - if atom.occupancy is not None: - occupancy = float(atom.occupancy) - except Exception: - pass - - displacement_type = "?" - u_iso = None - aniso = None - try: - displacement = atom.displacement_parameters - if displacement is not None: - u_iso = float(displacement.isotropic_equivalent) - if displacement.type == "Anisotropic": - displacement_type = "Uani" - values = displacement.values - uncertainties = displacement.uncertainties - aniso = { - "u11": values[0][0], "u22": values[1][1], "u33": values[2][2], - "u12": values[0][1], "u13": values[0][2], "u23": values[1][2], - "e11": uncertainties[0][0], "e22": uncertainties[1][1], "e33": uncertainties[2][2], - "e12": uncertainties[0][1], "e13": uncertainties[0][2], "e23": uncertainties[1][2], - } - else: - displacement_type = "Uiso" - except Exception: - pass - + displacement_type, u_iso, aniso = _extract_displacement(atom) disorder_group, disorder_assembly = disorder_map.get(label, (None, None)) atom_data.append({ "label": label, @@ -126,30 +131,34 @@ def export_full_cif_from_csd(refcode: str, output_path: Path) -> dict[str, Any]: "fx": float(fractional.x) if fractional else None, "fy": float(fractional.y) if fractional else None, "fz": float(fractional.z) if fractional else None, - "occupancy": occupancy, + "occupancy": _extract_occupancy(atom), "u_iso": u_iso, "displacement_type": displacement_type, "disorder_group": disorder_group, "disorder_assembly": disorder_assembly, "aniso": aniso, }) + return atom_data + +def _extract_bonds(molecule: Any) -> list[tuple[str, str, float, str]]: bonds: list[tuple[str, str, float, str]] = [] try: for bond in molecule.bonds: atom_1, atom_2 = bond.atoms bonds.append((str(atom_1.label), str(atom_2.label), float(bond.length), str(bond.bond_type))) except Exception: - pass + return bonds + return bonds + +def _build_header_lines(refcode: str, entry: Any, crystal: Any) -> list[str]: cell_a, cell_b, cell_c = crystal.cell_lengths alpha, beta, gamma = crystal.cell_angles - disorder_present = any(atom["disorder_group"] is not None for atom in atom_data) - - lines = [ + return [ f"data_{refcode}", "_audit_creation_method 'CSD Python API disorder workflow demo'", - f"_chemical_name_common '{(entry.chemical_name or '?').replace(chr(39), chr(39) + chr(39))}'", + f"_chemical_name_common '{_escape_cif_string(entry.chemical_name)}'", f"_chemical_formula_sum '{entry.formula or '?'}'", f"_cell_length_a {_format_float(cell_a)}", f"_cell_length_b {_format_float(cell_b)}", @@ -163,13 +172,20 @@ def export_full_cif_from_csd(refcode: str, output_path: Path) -> dict[str, Any]: "", ] - if crystal.symmetry_operators: - lines.extend(["loop_", "_symmetry_equiv_pos_as_xyz"]) - for operator in crystal.symmetry_operators: - lines.append(f" '{operator}'") - lines.append("") - lines.extend([ +def _build_symmetry_lines(crystal: Any) -> list[str]: + if not crystal.symmetry_operators: + return [] + lines = ["loop_", "_symmetry_equiv_pos_as_xyz"] + for operator in crystal.symmetry_operators: + lines.append(f" '{operator}'") + lines.append("") + return lines + + +def _build_atom_site_lines(atom_data: list[dict[str, Any]]) -> list[str]: + disorder_present = any(atom["disorder_group"] is not None for atom in atom_data) + lines = [ "loop_", "_atom_site_label", "_atom_site_type_symbol", @@ -179,7 +195,7 @@ def export_full_cif_from_csd(refcode: str, output_path: Path) -> dict[str, Any]: "_atom_site_occupancy", "_atom_site_U_iso_or_equiv", "_atom_site_thermal_displace_type", - ]) + ] if disorder_present: lines.extend([ "_atom_site_disorder_assembly", @@ -200,52 +216,96 @@ def export_full_cif_from_csd(refcode: str, output_path: Path) -> dict[str, Any]: f"{atom['displacement_type']:>5s}", ] if disorder_present: - row.append(f"{str(atom['disorder_assembly']) if atom['disorder_assembly'] is not None else '.':>4s}") - row.append(f"{str(atom['disorder_group']) if atom['disorder_group'] is not None else '.':>4s}") + row.extend([ + f"{str(atom['disorder_assembly']) if atom['disorder_assembly'] is not None else '.':>4s}", + f"{str(atom['disorder_group']) if atom['disorder_group'] is not None else '.':>4s}", + ]) lines.append(" " + " ".join(row)) lines.append("") + return lines + +def _build_anisotropic_lines(atom_data: list[dict[str, Any]]) -> list[str]: anisotropic_atoms = [atom for atom in atom_data if atom["aniso"]] - if anisotropic_atoms: - lines.extend([ - "loop_", - "_atom_site_aniso_label", - "_atom_site_aniso_U_11", - "_atom_site_aniso_U_22", - "_atom_site_aniso_U_33", - "_atom_site_aniso_U_12", - "_atom_site_aniso_U_13", - "_atom_site_aniso_U_23", - ]) - for atom in anisotropic_atoms: - u = atom["aniso"] - lines.append( - f" {atom['label']:<8s}" - f" {_format_esd(u['u11'], u['e11']):>12s}" - f" {_format_esd(u['u22'], u['e22']):>12s}" - f" {_format_esd(u['u33'], u['e33']):>12s}" - f" {_format_esd(u['u12'], u['e12']):>12s}" - f" {_format_esd(u['u13'], u['e13']):>12s}" - f" {_format_esd(u['u23'], u['e23']):>12s}" - ) - lines.append("") - - if bonds: - lines.extend([ - "loop_", - "_geom_bond_atom_site_label_1", - "_geom_bond_atom_site_label_2", - "_geom_bond_distance", - "_ccdc_geom_bond_type", - ]) - for atom_1, atom_2, distance, bond_type in bonds: - lines.append(f" {atom_1:<8s} {atom_2:<8s} {distance:8.4f} {bond_type}") - lines.append("") + if not anisotropic_atoms: + return [] + lines = [ + "loop_", + "_atom_site_aniso_label", + "_atom_site_aniso_U_11", + "_atom_site_aniso_U_22", + "_atom_site_aniso_U_33", + "_atom_site_aniso_U_12", + "_atom_site_aniso_U_13", + "_atom_site_aniso_U_23", + ] + for atom in anisotropic_atoms: + u = atom["aniso"] + lines.append( + f" {atom['label']:<8s}" + f" {_format_esd(u['u11'], u['e11']):>12s}" + f" {_format_esd(u['u22'], u['e22']):>12s}" + f" {_format_esd(u['u33'], u['e33']):>12s}" + f" {_format_esd(u['u12'], u['e12']):>12s}" + f" {_format_esd(u['u13'], u['e13']):>12s}" + f" {_format_esd(u['u23'], u['e23']):>12s}" + ) + lines.append("") + return lines - lines.append("#END") + +def _build_bond_lines(bonds: list[tuple[str, str, float, str]]) -> list[str]: + if not bonds: + return [] + lines = [ + "loop_", + "_geom_bond_atom_site_label_1", + "_geom_bond_atom_site_label_2", + "_geom_bond_distance", + "_ccdc_geom_bond_type", + ] + for atom_1, atom_2, distance, bond_type in bonds: + lines.append(f" {atom_1:<8s} {atom_2:<8s} {distance:8.4f} {bond_type}") + lines.append("") + return lines + + +def _write_cif_file(output_path: Path, lines: list[str]) -> None: output_path.parent.mkdir(parents=True, exist_ok=True) output_path.write_text("\n".join(lines), encoding="utf-8") + +def _validated_output_path(output_dir: Path, refcode: str) -> Path: + base_dir = output_dir.expanduser().resolve() + candidate = (base_dir / f"{refcode}.cif").resolve() + try: + candidate.relative_to(base_dir) + except ValueError as exc: + raise ValueError(f"Refusing to write outside output directory: {candidate}") from exc + return candidate + + +def export_full_cif_from_csd(refcode: str, output_path: Path) -> dict[str, Any]: + """Export a CSD structure to CIF with occupancy and disorder metadata.""" + from ccdc.io import EntryReader + + with EntryReader("CSD") as reader: + entry = reader.entry(refcode) + crystal = reader.crystal(refcode) + + molecule = crystal.disordered_molecule or crystal.molecule + disorder_map = _build_disorder_map(crystal) + atom_data = _extract_atom_data(molecule, disorder_map) + bonds = _extract_bonds(molecule) + lines = [] + lines.extend(_build_header_lines(refcode, entry, crystal)) + lines.extend(_build_symmetry_lines(crystal)) + lines.extend(_build_atom_site_lines(atom_data)) + lines.extend(_build_anisotropic_lines(atom_data)) + lines.extend(_build_bond_lines(bonds)) + lines.append("#END") + _write_cif_file(output_path, lines) + return { "refcode": refcode, "output": str(output_path), @@ -268,7 +328,7 @@ def main() -> None: print(f"Licence requirement: {LICENSE_REQUIREMENT}") for refcode in args.refcodes: - output_path = output_dir / f"{refcode}.cif" + output_path = _validated_output_path(output_dir, refcode) summary = export_full_cif_from_csd(refcode, output_path) print( f"Exported {summary['refcode']} -> {summary['output']} " From 6f8f1b5b038ad1673324d3df60c435903d90a66e Mon Sep 17 00:00:00 2001 From: SchrodingersCattt Date: Mon, 29 Jun 2026 10:10:23 +0000 Subject: [PATCH 09/15] Sanitize CLI inputs to address S8707 path injection rule - Validate refcode with strict alphanumeric regex before use as filename - Validate output-dir: reject path traversal (..), enforce cwd-relative - Break taint chain: CLI args no longer flow directly to filesystem ops --- .../demo_csd_disorder_workflow.py | 50 +++++++++++++------ 1 file changed, 36 insertions(+), 14 deletions(-) diff --git a/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py b/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py index 74a9a8b..e610e14 100644 --- a/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py +++ b/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py @@ -35,6 +35,7 @@ from __future__ import annotations import argparse +import re from pathlib import Path from typing import Any @@ -275,23 +276,45 @@ def _write_cif_file(output_path: Path, lines: list[str]) -> None: output_path.write_text("\n".join(lines), encoding="utf-8") -def _validated_output_path(output_dir: Path, refcode: str) -> Path: - base_dir = output_dir.expanduser().resolve() - candidate = (base_dir / f"{refcode}.cif").resolve() - try: - candidate.relative_to(base_dir) - except ValueError as exc: - raise ValueError(f"Refusing to write outside output directory: {candidate}") from exc - return candidate +_SAFE_REFCODE_RE = re.compile(r"^[A-Za-z][A-Za-z0-9]{1,19}$") -def export_full_cif_from_csd(refcode: str, output_path: Path) -> dict[str, Any]: +def _sanitize_refcode(refcode: str) -> str: + """Ensure refcode contains only safe alphanumeric characters.""" + if not _SAFE_REFCODE_RE.match(refcode): + raise ValueError( + f"Invalid refcode {refcode!r}: must be 2-20 alphanumeric characters starting with a letter" + ) + return refcode + + +def _safe_output_dir(raw_path: str) -> Path: + """Resolve and validate output directory, rejecting path traversal.""" + # Only allow relative paths under the current working directory + if raw_path != Path(raw_path).as_posix().replace("\\", "/"): + cleaned = Path(raw_path).as_posix() + else: + cleaned = raw_path + if ".." in cleaned.split("/"): + raise ValueError(f"Path traversal not allowed in output directory: {raw_path!r}") + resolved = Path.cwd() / cleaned + resolved = resolved.resolve() + cwd = Path.cwd().resolve() + if not str(resolved).startswith(str(cwd)): + raise ValueError(f"Output directory must be under working directory: {resolved}") + return resolved + + +def export_full_cif_from_csd(refcode: str, output_dir: Path) -> dict[str, Any]: """Export a CSD structure to CIF with occupancy and disorder metadata.""" from ccdc.io import EntryReader + safe_name = _sanitize_refcode(refcode) + output_path = output_dir / f"{safe_name}.cif" + with EntryReader("CSD") as reader: - entry = reader.entry(refcode) - crystal = reader.crystal(refcode) + entry = reader.entry(safe_name) + crystal = reader.crystal(safe_name) molecule = crystal.disordered_molecule or crystal.molecule disorder_map = _build_disorder_map(crystal) @@ -324,12 +347,11 @@ def parse_args() -> argparse.Namespace: def main() -> None: args = parse_args() - output_dir = Path(args.output_dir) + output_dir = _safe_output_dir(args.output_dir) print(f"Licence requirement: {LICENSE_REQUIREMENT}") for refcode in args.refcodes: - output_path = _validated_output_path(output_dir, refcode) - summary = export_full_cif_from_csd(refcode, output_path) + summary = export_full_cif_from_csd(refcode, output_dir) print( f"Exported {summary['refcode']} -> {summary['output']} " f"(atoms={summary['n_atoms']}, partial_occ={summary['n_partial_occupancy']}, bonds={summary['n_bonds']})" From a5b4740e38eeb2411dbe70f6d3d9f2b354ada157 Mon Sep 17 00:00:00 2001 From: Isaac Sugden Date: Fri, 10 Jul 2026 11:30:15 +0100 Subject: [PATCH 10/15] took out the License print statement NO_JIRA --- .../disorder_workflow_demo/demo_csd_disorder_workflow.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py b/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py index e610e14..d856324 100644 --- a/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py +++ b/scripts/materials/disorder_workflow_demo/demo_csd_disorder_workflow.py @@ -40,8 +40,6 @@ from typing import Any -LICENSE_REQUIREMENT = "CSD-Materials or another licence tier with CSD Python API database access" - def _escape_cif_string(value: str | None) -> str: return (value or "?").replace("'", "''") @@ -349,7 +347,6 @@ def main() -> None: args = parse_args() output_dir = _safe_output_dir(args.output_dir) - print(f"Licence requirement: {LICENSE_REQUIREMENT}") for refcode in args.refcodes: summary = export_full_cif_from_csd(refcode, output_dir) print( From 4ac59f71ac8169ca11c58982b89064a3a12f3a6b Mon Sep 17 00:00:00 2001 From: cringrose Date: Mon, 3 Aug 2026 15:51:27 +0100 Subject: [PATCH 11/15] NO_JIRA Added Pharmacophore generator script to create CrossMiner queries from Ligand Overlay outputs. Includes basic tests and examples. --- .../overlay_to_pharmacophore/ReadMe.md | 46 ++++ .../overlay_to_pharmacophore/cluster.py | 96 +++++++ .../crossminer_search.py | 18 ++ .../datastructures.py | 114 +++++++++ .../overlay_to_pharmacophore/main.py | 76 ++++++ .../overlay_to_pharmacophore/overlay.py | 105 ++++++++ .../pharmacophores/solution_pharm_01.mol2 | 63 +++++ .../tests/job1/solution_01.mol2 | 234 ++++++++++++++++++ .../overlay_to_pharmacophore/tests/tests.py | 105 ++++++++ .../overlay_to_pharmacophore/write_query.py | 110 ++++++++ 10 files changed, 967 insertions(+) create mode 100644 scripts/discovery/overlay_to_pharmacophore/ReadMe.md create mode 100644 scripts/discovery/overlay_to_pharmacophore/cluster.py create mode 100644 scripts/discovery/overlay_to_pharmacophore/crossminer_search.py create mode 100644 scripts/discovery/overlay_to_pharmacophore/datastructures.py create mode 100644 scripts/discovery/overlay_to_pharmacophore/main.py create mode 100644 scripts/discovery/overlay_to_pharmacophore/overlay.py create mode 100644 scripts/discovery/overlay_to_pharmacophore/tests/job1/pharmacophores/solution_pharm_01.mol2 create mode 100644 scripts/discovery/overlay_to_pharmacophore/tests/job1/solution_01.mol2 create mode 100644 scripts/discovery/overlay_to_pharmacophore/tests/tests.py create mode 100644 scripts/discovery/overlay_to_pharmacophore/write_query.py diff --git a/scripts/discovery/overlay_to_pharmacophore/ReadMe.md b/scripts/discovery/overlay_to_pharmacophore/ReadMe.md new file mode 100644 index 0000000..80ee9ce --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/ReadMe.md @@ -0,0 +1,46 @@ +# Pharmacophore Query Generator + +This tool provides the user with the ability to create pharmacophore queries from the results of a ligand overlay. + +The pharmacophore queries are produced to be used with CrossMiner. + +## Requirements + +- [CSD Python API](https://downloads.ccdc.cam.ac.uk/documentation/API/) installed. +- Access to CSD CrossMiner and the feature definitions (`.cpf`) files. +- Access to the CSD Ligand Overlay Tool. + +## Licensing Requirements + +CSD-Discovery, CSD-Enterprise and Research Partner suites would all be sufficient. + +## Instructions on Running + +### Feature Definitions + +The CrossMiner feature definition (`.cpf`) files are **not** shipped with this repo. +Supply the location of the feature definitions from your CrossMiner installation with +`-f`/`--feature_definitions`; this should be the directory containing the definition files +(either directly, or in `any`/`protein`/`small_molecule` subdirectories). +Usually, the location is `C:\users\\CCDC\ccdc-software\csd-crossminer\feature_definitions` + +``` +python main.py -i -o -f +``` + +The output folder (`-o`/`--output_folder`) is optional; if it is not supplied, the queries are +written to a `queries` folder created in the current directory. + +### Options for Ligand Overlay Output + +* `cluster`: Cluster the similar pharmacophore features based on proximity +* `projected`: Treat pharmacophore features as projected when appropriate e.g. acceptors and donors + +There is also the option to specify a specific Ligand Overlay from all the results. If this is not specified, all overlays are used. +When all overlays are used, the pharmacophore query will be a union of all the features from all the overlays. +These features are then clustered based on proximity AND prevalence. + +### Using the Queries Generated + +If you would like to use the queries generated with this tool, they can be opened in CrossMiner to run a search. +A file `crossminer_search.py` has also been provided which contains a Python function for the most simply kind of CrossMiner search. diff --git a/scripts/discovery/overlay_to_pharmacophore/cluster.py b/scripts/discovery/overlay_to_pharmacophore/cluster.py new file mode 100644 index 0000000..e49add2 --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/cluster.py @@ -0,0 +1,96 @@ +from collections import defaultdict + +import numpy as np + +from datastructures import PharmFeaturePoint + + +def cluster_features(features: list[PharmFeaturePoint]) -> list[PharmFeaturePoint]: + """Cluster all features where relevant.""" + groups = defaultdict(list) + for feature in features: + groups[feature.label].append(feature) + + clustered_features = [] + for label, group in groups.items(): + if label == 'hydrophobe': + clustered_features.extend(_cluster_similar_features(group)) + elif label in ('acceptor_projected', 'donor_projected'): + clustered_features.extend(_cluster_similar_features(group, check_vp=True)) + else: + clustered_features.extend(group) + return clustered_features + + +def _cluster_similar_features( + features: list[PharmFeaturePoint], cluster_radius: float = 2.0, check_vp: bool = False +) -> list[PharmFeaturePoint]: + """ + Cluster specific features that are close to each other. + For projected features, also check virtual point distances. + + Args: + features: List of features to cluster + cluster_radius: Distance between features to be included in clustering (Ang) + check_vp: Whether to also require virtual points to be within the radius + """ + if len(features) < 2: + return features + + clusters = _connected_components( + features, + lambda a, b: _is_close(a, b, cluster_radius, check_vp), + ) + return [_merge_cluster(cluster, check_vp) for cluster in clusters] + + +def _is_close( + a: PharmFeaturePoint, b: PharmFeaturePoint, cluster_radius: float, check_vp: bool +) -> bool: + """Whether two features are within ``cluster_radius`` (and, if ``check_vp``, their virtual points too).""" + if np.linalg.norm(a - b) >= cluster_radius: + return False + if check_vp and np.linalg.norm(a.virtual_point - b.virtual_point) >= cluster_radius: + return False + return True + + +def _connected_components( + features: list[PharmFeaturePoint], is_linked +) -> list[list[PharmFeaturePoint]]: + """ + Single-linkage grouping: features are placed in the same cluster if a chain of + ``is_linked`` neighbours connects them. + """ + unclustered = features.copy() + clusters = [] + + while unclustered: + cluster = [unclustered.pop(0)] + # Grow the cluster breadth-first: any unclustered feature linked to a member joins it. + i = 0 + while i < len(cluster): + member = cluster[i] + remaining = [] + for feature in unclustered: + if is_linked(feature, member): + cluster.append(feature) + else: + remaining.append(feature) + unclustered = remaining + i += 1 + clusters.append(cluster) + + return clusters + + +def _merge_cluster(cluster: list[PharmFeaturePoint], check_vp: bool) -> PharmFeaturePoint: + """Collapse a cluster into a single feature at its centroid (singletons are returned unchanged).""" + if len(cluster) == 1: + return cluster[0] + + centroid = np.mean([c.coordinates for c in cluster], axis=0).round(4) + vp_centroid = ( + np.mean([c.virtual_point for c in cluster], axis=0).round(4) if check_vp else None + ) + return PharmFeaturePoint(centroid, label=cluster[0].label, virtual_point=vp_centroid) diff --git a/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py b/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py new file mode 100644 index 0000000..962462e --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py @@ -0,0 +1,18 @@ +from pathlib import Path + +from ccdc.pharmacophore import Pharmacophore + +def search(query_file: Path, database_file: Path): + settings = Pharmacophore.Search.Settings() + settings.max_hit_structures = 20 + settings.max_hits_per_structure = 1 + settings.max_hit_rmsd = 1.0 + searcher = Pharmacophore.Search(settings) + feature_db = Pharmacophore.FeatureDatabase.from_file(database_file) + query = Pharmacophore.Query.from_file(str(query_file)) + hits = searcher.search( + model=query, + database=feature_db, + verbose=True, + ) + return hits diff --git a/scripts/discovery/overlay_to_pharmacophore/datastructures.py b/scripts/discovery/overlay_to_pharmacophore/datastructures.py new file mode 100644 index 0000000..0f451c2 --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/datastructures.py @@ -0,0 +1,114 @@ +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path +from typing import Iterable, Optional + +import numpy as np + + +@dataclass +class FeatureTolerances: + """ + Allowed feature types and their tolerances. + + Each field is a valid feature label; its value is the tolerance: + a single weight, or a (parent, virtual_point) pair for projected features. + """ + acceptor: float = 1.0 + acceptor_projected: tuple[float, float] = (0.8, 0.8) + donor_projected: tuple[float, float] = (1.0, 1.0) + hydrophobe: float = 1.0 + ring_planar_projected: tuple[float, float] = (1.0, 1.0) + ring_non_planar: tuple[float, float] = (1.0, 1.0) + halogen: float = 1.0 + + def __getitem__(self, key: str) -> float | tuple[float, float]: + return getattr(self, key) + + +class PharmFeaturePoint(np.ndarray): + def __new__( + cls, + *coordinates: float | Iterable[float], + label: Optional[str] = None, + virtual_point: Optional[np.ndarray] = None, + ): + if len(coordinates) == 1: + arr = np.asarray(coordinates[0], dtype=float) + elif len(coordinates) == 3: + arr = np.asarray(coordinates, dtype=float) + else: + raise TypeError("Coordinates must be either an iterable of length three, or three floats") + if arr.shape != (3,): + raise ValueError("Coordinates must be a 3-element array") + obj = arr.view(cls) + obj.label = label + obj.virtual_point = virtual_point + return obj + + def __array_finalize__(self, obj): + if obj is None: + return + self.label = getattr(obj, 'label', None) + self.virtual_point = getattr(obj, 'virtual_point', None) + + def __repr__(self): + return ( + f"PharmFeaturePoint({self.x}, {self.y}, {self.z}, " + f"label={self.label}, virtual_points={self.virtual_point}), " + ) + + def __str__(self): + return self.__repr__() + + @property + def coordinates(self) -> np.ndarray: + return np.asarray(self) + + @property + def x(self): + return self[0] + + @property + def y(self): + return self[1] + + @property + def z(self): + return self[2] + + @property + def tolerance(self) -> float | tuple[float, float]: + """ + Get the tolerance for the feature based on its label. + Returns: + A single tolerance for features with a single tolerance or a tuple for those with two tolerances. + """ + if self.label is None: + raise LookupError("Feature label must be set to determine tolerances.") + return FeatureTolerances()[self.label] + + @property + def weight_parent(self) -> float: + if isinstance(self.tolerance, float): + return self.tolerance + elif isinstance(self.tolerance, tuple): + return self.tolerance[0] + else: + raise ValueError("Incorrect tolerances loaded for feature point.") + + @property + def weight_vp(self) -> float: + return self.tolerance[1] + + + +@dataclass +class OverlayData: + input_folder: Path + output_folder: Path + # Pharmacophore file from the pharmacophores folder + pharm_file: Path + # Overlay solution file (the chosen one from the many produced) + overlay_file: Path diff --git a/scripts/discovery/overlay_to_pharmacophore/main.py b/scripts/discovery/overlay_to_pharmacophore/main.py new file mode 100644 index 0000000..1e4954b --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/main.py @@ -0,0 +1,76 @@ +import argparse +from pathlib import Path + +from cluster import cluster_features +from datastructures import OverlayData +from overlay import OverlayToPharmFeatures +from write_query import FeaturesToCrossMinerQuery + + +def str_to_bool(value: str) -> bool: + return value.lower() in {'t', 'true', '1', 'yes', 'y'} + + +def parse_args(): + parser = argparse.ArgumentParser( + description="Create Pharmacophore Features from a Ligand Overlay" + ) + parser.add_argument('-i', '--input_folder', type=str, help='Input file(s) path.') + parser.add_argument('-o', '--output_folder', type=str, default=None, + help="Output folder path. Defaults to a 'queries' folder in the current directory.") + parser.add_argument('-f', '--feature_definitions', type=str, required=True, + help='Path to the directory containing the CrossMiner feature definition (.cpf) files.') + parser.add_argument('-c', '--cluster', type=str_to_bool, default=False, + help='Cluster features if they are close together or common across multiple inputs.') + parser.add_argument('-p', '--projected', type=str_to_bool, default=False, + help='Use projected acceptor features or point features.') + parser.add_argument('-id', '--overlay_id', type=int, default=0, + help='Overlay ID to process. If 0 or not specified, all overlays will be processed.') + + return parser.parse_args() + + +def main(): + args = parse_args() + + input_folder = Path(args.input_folder) + if not input_folder.exists(): + raise FileNotFoundError(f"Input folder {input_folder} does not exist.") + + feature_definitions = Path(args.feature_definitions) + if not feature_definitions.is_dir(): + raise FileNotFoundError(f"Feature definitions folder {feature_definitions} does not exist.") + + output_folder = Path(args.output_folder) if args.output_folder else Path('queries') + output_folder.mkdir(parents=True, exist_ok=True) + + if (args.overlay_id == 0) or (args.overlay_id is None): + overlay_files = sorted(input_folder.glob('solution_*.mol2')) + pharm_files = sorted(input_folder.glob('pharmacophores/solution_pharm_*.mol2')) + else: + overlay_files = [input_folder / f'solution_{args.overlay_id:02}.mol2'] + pharm_files = [input_folder / f'pharmacophores/solution_pharm_{args.overlay_id:02}.mol2'] + feature_sets = [] + for pharm_file, overlay_file in zip(pharm_files, overlay_files): + overlay_data = OverlayData( + input_folder=input_folder, + output_folder=output_folder, + pharm_file=pharm_file, + overlay_file=overlay_file + ) + feature_sets.append(OverlayToPharmFeatures(overlay_data, projected=args.projected).features) + + for i, feature_set in enumerate(feature_sets, 1): + if args.cluster: + feature_set = cluster_features(feature_set) + + query = FeaturesToCrossMinerQuery( + pharm_feature_points=feature_set, + feature_definitions=feature_definitions, + output_file=output_folder / f'features_{i}.cm', + ) + query.write_feature_file() + + +if __name__ == '__main__': + main() diff --git a/scripts/discovery/overlay_to_pharmacophore/overlay.py b/scripts/discovery/overlay_to_pharmacophore/overlay.py new file mode 100644 index 0000000..1bc51bd --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/overlay.py @@ -0,0 +1,105 @@ +from functools import cached_property + +import numpy as np + +from datastructures import OverlayData, PharmFeaturePoint + + +class OverlayToPharmFeatures: + def __init__(self, overlay_data: OverlayData, projected: bool = True): + """ + Args: + overlay_data: Pharmacophore and overlay data for a single ligand overlay. + projected: Whether to use projected acceptor features + (since acceptors are flexible and vary across overlays, + it may be better to simply use point acceptors). + """ + self.overlay_data = overlay_data + self.projected = projected + + @cached_property + def features(self) -> list[PharmFeaturePoint]: + return self._feature_centres_from_pharm() + + def _read_atom_lines_from_pharm(self): + """ + Generator that yields lines from the ATOM block of a pharmacophore file. + This block is where the pharmacophore feature points are defined. + """ + reading_atoms = False + with open(self.overlay_data.pharm_file) as pharm_file: + for line in pharm_file: + line = line.strip() + + if line == "@ATOM": + reading_atoms = True + continue + elif line.startswith("@") and reading_atoms: + break # End of ATOM block + + if reading_atoms and line: + yield line.split() + + def _feature_centres_from_pharm(self) -> list[PharmFeaturePoint]: + """ + Parse the pharmacophore mol2 file and extract feature points centres. + + returns: List of PharmFeaturePoints, each labelled with the feature name + (e.g. 'donor_projected', 'acceptor', 'hydrophobe'), and with associated virtual points if present. + """ + feature_centres = [] + last_feature_point = None + + for parts in self._read_atom_lines_from_pharm(): + _feature_no, feature_name, x, y, z = parts[:5] + # Extract first 2/3 chars which describe the feature type + prefix = ''.join(c for c in feature_name if c.isalpha()) + + if prefix in ('DON', 'ACC', 'HY'): + label = self._prefix_to_label(prefix) + last_feature_point = PharmFeaturePoint( + (float(x), float(y), float(z)), + label=label, + ) + feature_centres.append(last_feature_point) + elif prefix == 'VP': + if last_feature_point is None: + raise ValueError("Virtual point found before its parent feature point.") + # Acceptor feature with a virtual point is an acceptor projected feature. + if last_feature_point.label == "acceptor": + # If preferring to use point acceptors, ignore the virtual point. + if not self.projected: + continue + last_feature_point.label = "acceptor_projected" + vp_coords = np.array([x, y, z], dtype=float) + # If the last feature has no virtual points, assign the virtual point to it. + if last_feature_point.virtual_point is None: + last_feature_point.virtual_point = vp_coords + # If the last feature does have virtual points, create a new feature. + else: + new_feature = last_feature_point.copy() + new_feature.virtual_point = vp_coords + feature_centres.append(new_feature) + elif prefix == 'NL': + if last_feature_point is None: + raise ValueError("Virtual point found before its parent feature point.") + last_feature_point.virtual_point = np.array([x, y, z], dtype=float) + last_feature_point.label = 'ring_planar_projected' + else: + raise ValueError(f"Unexpected feature type: {feature_name}") + return feature_centres + + @staticmethod + def _prefix_to_label(prefix: str) -> str: + """ + Convert a feature prefix to its corresponding label. + Acceptor may later be converted to acceptor_projected if it has a virtual point. + Args: + prefix: Prefix to convert to label + returns: Label for the feature type, e.g. 'acceptor', 'donor_projected', 'hydrophobe' + """ + return { + 'ACC': 'acceptor', + 'DON': 'donor_projected', + 'HY': 'hydrophobe', + }[prefix] diff --git a/scripts/discovery/overlay_to_pharmacophore/tests/job1/pharmacophores/solution_pharm_01.mol2 b/scripts/discovery/overlay_to_pharmacophore/tests/job1/pharmacophores/solution_pharm_01.mol2 new file mode 100644 index 0000000..71f5d16 --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/tests/job1/pharmacophores/solution_pharm_01.mol2 @@ -0,0 +1,63 @@ +@MOLECULE +pharmacophore + 20 12 1 0 0 +SMALL +NO_CHARGES +**** +Generated from the CSD + +@ATOM + 1 DON8.0_2 4.0507 8.2244 -1.3153 Ar 1 RES1 0.0000 + 2 VP1.0 3.7119 6.2366 -3.5091 F 1 RES1 0.0000 + 3 VP1.0 3.4830 11.1446 -0.9341 F 1 RES1 0.0000 + 4 VP1.0 5.0321 6.4941 -3.5384 F 1 RES1 0.0000 + 5 VP1.0 2.4551 6.6287 -3.1080 F 1 RES1 0.0000 + 6 ACC8.0_2 1.5469 6.6465 -0.4394 B 1 RES1 0.0000 + 7 VP1.0 -0.6255 8.6951 -0.6938 F 1 RES1 0.0000 + 8 VP1.0 -1.2681 7.3618 -1.1537 F 1 RES1 0.0000 + 9 VP1.0 -1.1566 5.8367 -1.4214 F 1 RES1 0.0000 + 10 ACC8.0_2 4.0507 8.2244 -1.3153 B 1 RES1 0.0000 + 11 VP1.0 3.7119 6.2366 -3.5091 F 1 RES1 0.0000 + 12 VP1.0 3.4830 11.1446 -0.9341 F 1 RES1 0.0000 + 13 VP1.0 5.0321 6.4941 -3.5384 F 1 RES1 0.0000 + 14 VP1.0 2.4551 6.6287 -3.1080 F 1 RES1 0.0000 + 15 HY4.0_2 1.6670 2.6310 -0.9036 C.3 1 RES1 0.0000 + 16 HY4.0_2 1.2974 1.8244 -0.9805 C.3 1 RES1 0.0000 + 17 HY4.0_2 2.0589 2.5883 -0.5021 C.3 1 RES1 0.0000 + 18 HY4.0_2 3.5562 7.3509 0.3815 C.3 1 RES1 0.0000 + 19 HY11.9_2 2.4272 9.3712 1.8815 C.2 1 RES1 0.0000 + 20 NL 3.2267 9.9618 1.7718 F 1 RES1 0.0000 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 1 4 1 + 4 1 5 1 + 5 6 7 1 + 6 6 8 1 + 7 6 9 1 + 8 10 11 1 + 9 10 12 1 + 10 10 13 1 + 11 10 14 1 + 12 19 20 1 +@SUBSTRUCTURE + 1 RES1 1 GROUP 0 **** **** 0 +# ATOMS1 1 17 2 25 +# COORDS1 4.0479 4.0535 8.2010 8.2477 -1.5472 -1.0834 +# ATOMS6 1 14 2 13 +# COORDS6 1.4558 1.6380 6.6269 6.6661 -0.5341 -0.3446 +# ATOMS10 1 17 2 25 +# COORDS10 4.0479 4.0535 8.2010 8.2477 -1.5472 -1.0834 +# ATOMS15 1 5 1 6 1 7 1 8 1 9 1 10 1 11 2 3 2 4 2 5 +# ATOMS15 2 6 2 7 2 8 2 9 +# COORDS15 0.1104 3.2954 1.6196 4.0590 -2.3001 0.4906 +# ATOMS16 1 4 1 5 1 6 1 7 1 8 2 2 2 3 2 7 2 8 2 9 +# COORDS16 0.1104 2.4723 1.1327 2.2710 -2.3001 0.1771 +# ATOMS17 1 4 1 5 1 8 1 9 1 10 1 11 2 2 2 3 2 4 2 5 +# ATOMS17 2 6 2 7 +# COORDS17 0.7974 3.2954 1.1327 4.0590 -1.4998 0.4906 +# ATOMS18 1 13 1 15 1 16 1 18 2 12 2 14 2 15 2 20 +# COORDS18 2.6111 4.8913 6.3101 8.4335 -0.4078 1.3029 +# ATOMS19 1 18 1 19 1 20 1 21 1 22 1 23 2 20 2 21 2 22 2 23 +# ATOMS19 2 24 +# COORDS19 1.5791 3.2261 8.4084 10.5289 0.5379 3.2697 diff --git a/scripts/discovery/overlay_to_pharmacophore/tests/job1/solution_01.mol2 b/scripts/discovery/overlay_to_pharmacophore/tests/job1/solution_01.mol2 new file mode 100644 index 0000000..3401d35 --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/tests/job1/solution_01.mol2 @@ -0,0 +1,234 @@ +@MOLECULE +ipratroponium + 54 56 1 0 0 +SMALL +USER_CHARGES +**** +Generated from the CSD + +@ATOM + 1 C1 0.4143 -0.6927 0.8235 C.3 1 RES1 0.0000 + 2 C2 1.2898 -0.2841 -0.3556 C.3 1 RES1 0.0000 + 3 C3 2.3978 -1.3166 -0.5276 C.3 1 RES1 0.0000 + 4 N4 1.8801 1.1327 -0.1432 N.4 1 RES1 1.0000 + 5 C5 0.7974 2.1776 0.0874 C.3 1 RES1 0.0000 + 6 C6 0.1175 2.2710 -1.2830 C.3 1 RES1 0.0000 + 7 C7 1.2245 1.8950 -2.3001 C.3 1 RES1 0.0000 + 8 C8 2.4723 1.6196 -1.4584 C.3 1 RES1 0.0000 + 9 C9 3.2954 2.8815 -1.1841 C.3 1 RES1 0.0000 + 10 C10 2.4141 3.9967 -0.6094 C.3 1 RES1 0.0000 + 11 C11 1.4899 3.4733 0.4906 C.3 1 RES1 0.0000 + 12 O12 3.2402 5.0698 -0.0540 O.3 1 RES1 0.0000 + 13 C13 2.7385 6.3101 0.0254 C.2 1 RES1 0.0000 + 14 O14 1.6380 6.6269 -0.3446 O.2 1 RES1 0.0000 + 15 C15 3.7487 7.2538 0.6507 C.3 1 RES1 0.0000 + 16 C16 4.7309 7.7403 -0.4078 C.3 1 RES1 0.0000 + 17 O17 4.0479 8.2477 -1.5472 O.3 1 RES1 0.0000 + 18 C18 3.0043 8.4084 1.3029 C.ar 1 RES1 0.0000 + 19 C19 2.2545 9.2919 0.5379 C.ar 1 RES1 0.0000 + 20 C20 1.5791 10.3453 1.1324 C.ar 1 RES1 0.0000 + 21 C21 1.6455 10.5289 2.4981 C.ar 1 RES1 0.0000 + 22 C22 2.3880 9.6590 3.2697 C.ar 1 RES1 0.0000 + 23 C23 3.0638 8.6053 2.6762 C.ar 1 RES1 0.0000 + 24 C24 2.9214 1.0709 0.9623 C.3 1 RES1 0.0000 + 25 H25 -0.4701 -0.0620 0.8842 H 1 RES1 0.0000 + 26 H26 0.9588 -0.6057 1.7611 H 1 RES1 0.0000 + 27 H27 0.0871 -1.7251 0.7191 H 1 RES1 0.0000 + 28 H28 0.6647 -0.3122 -1.2603 H 1 RES1 0.0000 + 29 H29 3.0775 -1.0241 -1.3253 H 1 RES1 0.0000 + 30 H30 1.9805 -2.2895 -0.7786 H 1 RES1 0.0000 + 31 H31 2.9749 -1.4245 0.3881 H 1 RES1 0.0000 + 32 H32 0.0947 1.8529 0.8533 H 1 RES1 0.0000 + 33 H33 -0.2428 3.2826 -1.4645 H 1 RES1 0.0000 + 34 H34 -0.7153 1.5722 -1.3464 H 1 RES1 0.0000 + 35 H35 1.4040 2.7184 -2.9896 H 1 RES1 0.0000 + 36 H36 0.9400 1.0065 -2.8618 H 1 RES1 0.0000 + 37 H37 3.0875 0.8468 -1.9170 H 1 RES1 0.0000 + 38 H38 3.7422 3.2249 -2.1159 H 1 RES1 0.0000 + 39 H39 4.0824 2.6452 -0.4695 H 1 RES1 0.0000 + 40 H40 1.8033 4.4069 -1.4121 H 1 RES1 0.0000 + 41 H41 2.0794 3.2920 1.3880 H 1 RES1 0.0000 + 42 H42 0.7317 4.2263 0.7003 H 1 RES1 0.0000 + 43 H43 4.3077 6.7110 1.4285 H 1 RES1 0.0000 + 44 H44 5.3796 6.9097 -0.7126 H 1 RES1 0.0000 + 45 H45 5.3695 8.5244 0.0174 H 1 RES1 0.0000 + 46 H46 4.7082 8.5543 -2.2167 H 1 RES1 0.0000 + 47 H47 2.1970 9.1558 -0.5370 H 1 RES1 0.0000 + 48 H48 0.9964 11.0276 0.5223 H 1 RES1 0.0000 + 49 H49 1.1161 11.3539 2.9633 H 1 RES1 0.0000 + 50 H50 2.4428 9.8005 4.3440 H 1 RES1 0.0000 + 51 H51 3.6457 7.9267 3.2909 H 1 RES1 0.0000 + 52 H52 2.4749 0.7263 1.8923 H 1 RES1 0.0000 + 53 H53 3.3699 2.0423 1.1516 H 1 RES1 0.0000 + 54 H54 3.7219 0.3833 0.6985 H 1 RES1 0.0000 +@BOND + 1 1 2 1 + 2 2 3 1 + 3 2 4 1 + 4 4 5 1 + 5 5 6 1 + 6 6 7 1 + 7 7 8 1 + 8 4 8 1 + 9 8 9 1 + 10 9 10 1 + 11 10 11 1 + 12 5 11 1 + 13 10 12 1 + 14 12 13 1 + 15 13 14 2 + 16 13 15 1 + 17 15 16 1 + 18 16 17 1 + 19 15 18 1 + 20 18 19 ar + 21 19 20 ar + 22 20 21 ar + 23 21 22 ar + 24 22 23 ar + 25 18 23 ar + 26 4 24 1 + 27 1 25 1 + 28 1 26 1 + 29 1 27 1 + 30 2 28 1 + 31 3 29 1 + 32 3 30 1 + 33 3 31 1 + 34 5 32 1 + 35 6 33 1 + 36 6 34 1 + 37 7 35 1 + 38 7 36 1 + 39 8 37 1 + 40 9 38 1 + 41 9 39 1 + 42 10 40 1 + 43 11 41 1 + 44 11 42 1 + 45 15 43 1 + 46 16 44 1 + 47 16 45 1 + 48 17 46 1 + 49 19 47 1 + 50 20 48 1 + 51 21 49 1 + 52 22 50 1 + 53 23 51 1 + 54 24 52 1 + 55 24 53 1 + 56 24 54 1 +@SUBSTRUCTURE + 1 RES1 1 GROUP 0 **** **** 0 +@MOLECULE +tiotroponium + 48 52 1 0 0 +SMALL +USER_CHARGES +**** +Generated from the CSD + +@ATOM + 1 C1 1.5214 -0.1775 -0.0873 C.3 1 RES1 0.0000 + 2 N2 2.0030 1.2552 -0.0998 N.4 1 RES1 1.0000 + 3 C3 0.8928 2.2223 0.1771 C.3 1 RES1 0.0000 + 4 C4 1.5247 3.5898 0.4313 C.3 1 RES1 0.0000 + 5 C5 2.3168 4.0590 -0.8009 C.3 1 RES1 0.0000 + 6 C6 3.2057 2.9780 -1.4159 C.3 1 RES1 0.0000 + 7 C7 2.4140 1.6737 -1.4998 C.3 1 RES1 0.0000 + 8 C8 1.0622 1.8176 -2.1727 C.3 1 RES1 0.0000 + 9 C9 0.1104 2.1792 -1.1126 C.3 1 RES1 0.0000 + 10 O10 0.3365 3.0608 -2.2197 O.3 1 RES1 0.0000 + 11 O11 3.1602 5.1905 -0.4171 O.3 1 RES1 0.0000 + 12 C12 2.6111 6.4096 -0.3315 C.2 1 RES1 0.0000 + 13 O13 1.4558 6.6661 -0.5341 O.2 1 RES1 0.0000 + 14 C14 3.6619 7.4673 0.0699 C.3 1 RES1 0.0000 + 15 C15 4.8913 6.7838 0.6568 C.2 1 RES1 0.0000 + 16 C16 4.9473 5.8556 1.6681 C.2 1 RES1 0.0000 + 17 C17 6.2767 5.4539 1.9307 C.2 1 RES1 0.0000 + 18 C18 7.1812 6.0616 1.1413 C.2 1 RES1 0.0000 + 19 S19 6.4499 7.1454 0.0489 S.3 1 RES1 0.0000 + 20 C20 3.0632 8.4335 1.0857 C.2 1 RES1 0.0000 + 21 C21 3.2261 8.4489 2.4496 C.2 1 RES1 0.0000 + 22 C22 2.5043 9.5073 3.0467 C.2 1 RES1 0.0000 + 23 C23 1.8260 10.2557 2.1576 C.2 1 RES1 0.0000 + 24 S24 2.0400 9.7010 0.5610 S.3 1 RES1 0.0000 + 25 O25 4.0535 8.2010 -1.0834 O.3 1 RES1 0.0000 + 26 C26 3.1383 1.2852 0.8931 C.3 1 RES1 0.0000 + 27 H27 2.3197 -0.8530 -0.3875 H 1 RES1 0.0000 + 28 H28 0.6854 -0.3205 -0.7683 H 1 RES1 0.0000 + 29 H29 1.1955 -0.4635 0.9106 H 1 RES1 0.0000 + 30 H30 0.2859 1.9006 1.0222 H 1 RES1 0.0000 + 31 H31 2.1982 3.5196 1.2842 H 1 RES1 0.0000 + 32 H32 0.7389 4.3121 0.6478 H 1 RES1 0.0000 + 33 H33 1.6102 4.3964 -1.5576 H 1 RES1 0.0000 + 34 H34 4.0868 2.8312 -0.7929 H 1 RES1 0.0000 + 35 H35 3.5142 3.2832 -2.4147 H 1 RES1 0.0000 + 36 H36 3.0035 0.8936 -1.9793 H 1 RES1 0.0000 + 37 H37 0.7633 1.0360 -2.8696 H 1 RES1 0.0000 + 38 H38 -0.8366 1.6417 -1.0965 H 1 RES1 0.0000 + 39 H39 4.0794 5.4783 2.1988 H 1 RES1 0.0000 + 40 H40 6.5414 4.7281 2.6925 H 1 RES1 0.0000 + 41 H41 8.2501 5.8827 1.1942 H 1 RES1 0.0000 + 42 H42 3.8331 7.7374 2.9999 H 1 RES1 0.0000 + 43 H43 2.4978 9.6973 4.1150 H 1 RES1 0.0000 + 44 H44 1.2151 11.1113 2.4258 H 1 RES1 0.0000 + 45 H45 4.4434 7.5826 -1.7496 H 1 RES1 0.0000 + 46 H46 2.7864 1.0180 1.8873 H 1 RES1 0.0000 + 47 H47 3.5853 2.2755 0.9482 H 1 RES1 0.0000 + 48 H48 3.9160 0.5793 0.6091 H 1 RES1 0.0000 +@BOND + 1 1 2 1 + 2 2 3 1 + 3 3 4 1 + 4 4 5 1 + 5 5 6 1 + 6 6 7 1 + 7 2 7 1 + 8 7 8 1 + 9 8 9 1 + 10 3 9 1 + 11 9 10 1 + 12 8 10 1 + 13 5 11 1 + 14 11 12 1 + 15 12 13 2 + 16 12 14 1 + 17 14 15 1 + 18 15 16 2 + 19 16 17 1 + 20 17 18 2 + 21 18 19 1 + 22 15 19 1 + 23 14 20 1 + 24 20 21 2 + 25 21 22 1 + 26 22 23 2 + 27 23 24 1 + 28 20 24 1 + 29 14 25 1 + 30 2 26 1 + 31 1 27 1 + 32 1 28 1 + 33 1 29 1 + 34 3 30 1 + 35 4 31 1 + 36 4 32 1 + 37 5 33 1 + 38 6 34 1 + 39 6 35 1 + 40 7 36 1 + 41 8 37 1 + 42 9 38 1 + 43 16 39 1 + 44 17 40 1 + 45 18 41 1 + 46 21 42 1 + 47 22 43 1 + 48 23 44 1 + 49 25 45 1 + 50 26 46 1 + 51 26 47 1 + 52 26 48 1 +@SUBSTRUCTURE + 1 RES1 1 GROUP 0 **** **** 0 diff --git a/scripts/discovery/overlay_to_pharmacophore/tests/tests.py b/scripts/discovery/overlay_to_pharmacophore/tests/tests.py new file mode 100644 index 0000000..2406ea5 --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/tests/tests.py @@ -0,0 +1,105 @@ +from pathlib import Path + +import pytest + +from cluster import cluster_features +from datastructures import OverlayData, PharmFeaturePoint +from overlay import OverlayToPharmFeatures +from write_query import FeaturesToCrossMinerQuery + +TESTS_DIR = Path(__file__).parent +REPO_ROOT = TESTS_DIR.parent +JOB_FOLDER = TESTS_DIR / "job1" + + +@pytest.fixture(autouse=True) +def run_from_repo_root(monkeypatch): + """Run every test from the repository root so ``feature_definitions/`` resolves.""" + monkeypatch.chdir(REPO_ROOT) + + +@pytest.fixture +def overlay_data(tmp_path) -> OverlayData: + """Overlay/pharmacophore file locations for the bundled ``job1`` example.""" + return OverlayData( + input_folder=JOB_FOLDER, + output_folder=tmp_path, + pharm_file=JOB_FOLDER / "pharmacophores/solution_pharm_01.mol2", + overlay_file=JOB_FOLDER / "solution_01.mol2", + ) + + +@pytest.fixture +def features(overlay_data) -> list[PharmFeaturePoint]: + """Raw (unclustered) features extracted from the example pharmacophore.""" + return OverlayToPharmFeatures(overlay_data).features + + +@pytest.fixture +def feature_definitions(features, tmp_path) -> Path: + """ + Stub CrossMiner feature definition directory. + + The real ``.cpf`` files ship with CrossMiner rather than this repo, so minimal + placeholders are generated for each label present in the example pharmacophore. + """ + definitions = tmp_path / "feature_definitions" / "any" + definitions.mkdir(parents=True) + for label in {f.label for f in features}: + (definitions / f"features_{label}.cpf").write_text(f"SUBSTRUCTURE {label}\n") + return definitions.parent + + +def test_features_extracted_from_pharmacophore(features): + """Parsing the pharmacophore mol2 yields labelled feature points.""" + assert len(features) > 0, "No features extracted from pharmacophore." + assert all(feature.label is not None for feature in features) + + +def test_projected_features_have_virtual_points(features): + """Every feature labelled as projected carries a virtual point.""" + projected = [f for f in features if f.label.endswith("projected")] + assert projected, "No projected features extracted from pharmacophore." + assert all(f.virtual_point is not None for f in projected) + + +def test_unprojected_acceptors_have_no_virtual_points(overlay_data): + """With ``projected=False`` acceptors are kept as plain points.""" + features = OverlayToPharmFeatures(overlay_data, projected=False).features + acceptors = [f for f in features if f.label.startswith("acceptor")] + assert acceptors, "No acceptor features extracted from pharmacophore." + assert all(f.label == "acceptor" and f.virtual_point is None for f in acceptors) + + +def test_clustering_reduces_features(features): + """Clustering merges nearby features without dropping any feature type.""" + clustered = cluster_features(features) + assert 0 < len(clustered) <= len(features) + assert {f.label for f in clustered} == {f.label for f in features} + + +def test_query_file_written(features, feature_definitions, tmp_path): + """A CrossMiner query file is written containing the clustered features.""" + output_file = tmp_path / "query.cm" + clustered = cluster_features(features) + + FeaturesToCrossMinerQuery( + clustered, feature_definitions=feature_definitions, output_file=output_file + ).write_feature_file() + + assert output_file.exists() + contents = output_file.read_text() + assert "FEATURE_LIBRARY_START" in contents + assert "FEATURE_LIBRARY_END" in contents + assert contents.count("PHARMACOPHORE_FEATURE ") == len(clustered) + + +def test_missing_feature_definitions_raises(features, tmp_path): + """An error is raised when the supplied feature definitions directory is absent.""" + with pytest.raises(FileNotFoundError): + FeaturesToCrossMinerQuery( + features, + feature_definitions=tmp_path / "does_not_exist", + output_file=tmp_path / "query.cm", + ) + diff --git a/scripts/discovery/overlay_to_pharmacophore/write_query.py b/scripts/discovery/overlay_to_pharmacophore/write_query.py new file mode 100644 index 0000000..4d231b4 --- /dev/null +++ b/scripts/discovery/overlay_to_pharmacophore/write_query.py @@ -0,0 +1,110 @@ +from pathlib import Path + +from datastructures import PharmFeaturePoint + + +class FeaturesToCrossMinerQuery: + """Main class for converting a set of features to a CrossMiner feature file.""" + + TEMPLATE_FEATURES_FROM_FILE = """ + FEATURE_SUBSTRUCTURE_START + + {} + FEATURE_SUBSTRUCTURE_END + """ + def __init__( + self, + pharm_feature_points: list[PharmFeaturePoint], + feature_definitions: Path, + output_file: Path = Path('test.cm'), + ): + """ + Args: + pharm_feature_points: List of PharmFeaturePoints to be used in the CrossMiner query. + feature_definitions: Directory containing the CrossMiner feature definition (``.cpf``) files. + These are not shipped with this repo; supply the location of your CrossMiner installation's + feature definitions (the directory containing the ``any``, ``protein`` and + ``small_molecule`` subdirectories). + output_file: Output file path where the features will be saved. + """ + self.pharm_feature_points: list[PharmFeaturePoint] = pharm_feature_points + self.feature_definitions = Path(feature_definitions) + if not self.feature_definitions.is_dir(): + raise FileNotFoundError( + f"Feature definitions directory {self.feature_definitions} does not exist." + ) + self.output_file = output_file + + def _feature_definition_file(self, feature: PharmFeaturePoint) -> Path: + """ + Locate the ``.cpf`` definition file for a feature within the supplied feature definitions directory. + + Args: + feature: PharmFeaturePoint whose label determines the definition file name. + returns: Path to the matching ``.cpf`` file. + """ + if feature.label is None: + raise LookupError("Feature label must be set to locate its feature definition file.") + + filename = f'features_{feature.label}.cpf' + matches = sorted(self.feature_definitions.glob(f'*/{filename}')) + if not matches: + # Also allow a flat directory of definition files. + flat = self.feature_definitions / filename + if flat.is_file(): + return flat + raise FileNotFoundError( + f"No feature definition file '{filename}' found in {self.feature_definitions}." + ) + return matches[0] + + + @staticmethod + def _create_feature_vector_text(feature: PharmFeaturePoint) -> str: + """ + Create a formatted string for a feature vector based on the PharmFeaturePoint object. + Args: + feature: PharmFeaturePoint containing coordinates, label, and virtual point if applicable. + """ + + result = ( + f"PHARMACOPHORE_FEATURE {feature.label}\n" + f"PHARMACOPHORE_SPHERE {' '.join(str(i) for i in feature.coordinates)} {feature.weight_parent}\n" + ) + if feature.virtual_point is not None: + result += f"PHARMACOPHORE_SPHERE {' '.join(str(i) for i in feature.virtual_point)} {feature.weight_vp}\n" + + result += ( + f"PHARMACOPHORE_FEATURE_SMALL_MOLECULE\n" + f"PHARMACOPHORE_FEATURE_DESCRIPTION {feature.label}\n\n" + ) + + return result + + def write_feature_file(self): + """ + Create a feature file with proper whitespace formatting. + """ + output_data = 'FEATURE_LIBRARY_START\n' + + # Collect unique feature types and their corresponding files + feature_files = dict() + for feature in self.pharm_feature_points: + if feature.label in feature_files: + continue + feature_file = self._feature_definition_file(feature) + feature_files[feature.label] = feature_file + + # Append the feature data to the output stream + with open(feature_file) as ff: + output_data += self.TEMPLATE_FEATURES_FROM_FILE.format( + " ".join(ff.readlines()) + ) + + output_data += '\nFEATURE_LIBRARY_END\n\n' + + for feature in self.pharm_feature_points: + output_data += self._create_feature_vector_text(feature) + + with open(self.output_file, 'w') as output_file: + output_file.write(output_data) From 7632633c2f831f7a29ae20fa832e3e82f57c538d Mon Sep 17 00:00:00 2001 From: ccdc-github-admin <62256659+ccdc-github-admin@users.noreply.github.com> Date: Mon, 3 Aug 2026 23:34:39 +0100 Subject: [PATCH 12/15] ci(NO_JIRA): Updated .github/workflows/lint-python.yml --- .github/workflows/lint-python.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/lint-python.yml b/.github/workflows/lint-python.yml index 245010d..6d9374c 100644 --- a/.github/workflows/lint-python.yml +++ b/.github/workflows/lint-python.yml @@ -37,7 +37,7 @@ jobs: fi - name: Set up Python environment - uses: actions/setup-python@v6 + uses: actions/setup-python@v7 with: python-version: ${{ matrix.python_version }} From 3eff9d4beeb6bd2aea26135a6faaf828c3e5bb50 Mon Sep 17 00:00:00 2001 From: ccdc-github-admin <62256659+ccdc-github-admin@users.noreply.github.com> Date: Mon, 3 Aug 2026 23:34:44 +0100 Subject: [PATCH 13/15] ci(NO_JIRA): Updated .github/workflows/common_ccdc_status_checks.yml --- .github/workflows/common_ccdc_status_checks.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/common_ccdc_status_checks.yml b/.github/workflows/common_ccdc_status_checks.yml index 7ec5ce2..126cf39 100644 --- a/.github/workflows/common_ccdc_status_checks.yml +++ b/.github/workflows/common_ccdc_status_checks.yml @@ -14,7 +14,7 @@ jobs: with: ref: ${{ github.head_ref }} fetch-depth: 0 - - uses: actions/setup-python@v6 + - uses: actions/setup-python@v7 with: python-version: "3.11" - name: Get the commit message From 89520525ba923e872e87b1f93a961ecb1651d871 Mon Sep 17 00:00:00 2001 From: cringrose Date: Fri, 7 Aug 2026 10:11:44 +0100 Subject: [PATCH 14/15] NO_JIRA Added Pharmacophore generator script to create CrossMiner queries from Ligand Overlay outputs. Includes basic tests and examples. --- scripts/discovery/overlay_to_pharmacophore/ReadMe.md | 12 +++++++----- .../overlay_to_pharmacophore/crossminer_search.py | 6 ++++++ .../overlay_to_pharmacophore/datastructures.py | 6 +++--- scripts/discovery/overlay_to_pharmacophore/main.py | 12 +++++++++++- .../overlay_to_pharmacophore/tests/tests.py | 5 ++--- .../overlay_to_pharmacophore/write_query.py | 2 +- 6 files changed, 30 insertions(+), 13 deletions(-) diff --git a/scripts/discovery/overlay_to_pharmacophore/ReadMe.md b/scripts/discovery/overlay_to_pharmacophore/ReadMe.md index 80ee9ce..1bef1b7 100644 --- a/scripts/discovery/overlay_to_pharmacophore/ReadMe.md +++ b/scripts/discovery/overlay_to_pharmacophore/ReadMe.md @@ -9,8 +9,8 @@ The pharmacophore queries are produced to be used with CrossMiner. - [CSD Python API](https://downloads.ccdc.cam.ac.uk/documentation/API/) installed. - Access to CSD CrossMiner and the feature definitions (`.cpf`) files. - Access to the CSD Ligand Overlay Tool. - -## Licensing Requirements + +## Licensing Requirements CSD-Discovery, CSD-Enterprise and Research Partner suites would all be sufficient. @@ -36,11 +36,13 @@ written to a `queries` folder created in the current directory. * `cluster`: Cluster the similar pharmacophore features based on proximity * `projected`: Treat pharmacophore features as projected when appropriate e.g. acceptors and donors -There is also the option to specify a specific Ligand Overlay from all the results. If this is not specified, all overlays are used. +There is also the option to specify a specific Ligand Overlay from all the results. If this is not specified, all +overlays are used. When all overlays are used, the pharmacophore query will be a union of all the features from all the overlays. These features are then clustered based on proximity AND prevalence. ### Using the Queries Generated -If you would like to use the queries generated with this tool, they can be opened in CrossMiner to run a search. -A file `crossminer_search.py` has also been provided which contains a Python function for the most simply kind of CrossMiner search. +If you would like to use the queries generated with this tool, they can be opened in CrossMiner to run a search. +A file `crossminer_search.py` has also been provided which contains a Python function for the most simply kind of +CrossMiner search. diff --git a/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py b/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py index 962462e..c35ca78 100644 --- a/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py +++ b/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py @@ -2,7 +2,13 @@ from ccdc.pharmacophore import Pharmacophore + def search(query_file: Path, database_file: Path): + """ + This is here as an example for users to perform a simple CrossMiner + pharmacophore search using the CCDC Python API. + It is not used in the main workflow. + """ settings = Pharmacophore.Search.Settings() settings.max_hit_structures = 20 settings.max_hits_per_structure = 1 diff --git a/scripts/discovery/overlay_to_pharmacophore/datastructures.py b/scripts/discovery/overlay_to_pharmacophore/datastructures.py index 0f451c2..5f3bfe4 100644 --- a/scripts/discovery/overlay_to_pharmacophore/datastructures.py +++ b/scripts/discovery/overlay_to_pharmacophore/datastructures.py @@ -53,13 +53,13 @@ def __array_finalize__(self, obj): self.label = getattr(obj, 'label', None) self.virtual_point = getattr(obj, 'virtual_point', None) - def __repr__(self): + def __repr__(self) -> str: return ( f"PharmFeaturePoint({self.x}, {self.y}, {self.z}, " - f"label={self.label}, virtual_points={self.virtual_point}), " + f"label={self.label}, virtual_point={self.virtual_point}), " ) - def __str__(self): + def __str__(self) -> str: return self.__repr__() @property diff --git a/scripts/discovery/overlay_to_pharmacophore/main.py b/scripts/discovery/overlay_to_pharmacophore/main.py index 1e4954b..74c11dc 100644 --- a/scripts/discovery/overlay_to_pharmacophore/main.py +++ b/scripts/discovery/overlay_to_pharmacophore/main.py @@ -1,3 +1,12 @@ +""" +This script can be used for any purpose without limitation subject to the +conditions at https://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx +This permission notice and the following statement of attribution must be +included in all copies or substantial portions of this script. + +"date of creation": created by the Cambridge Crystallographic Data Centre +""" + import argparse from pathlib import Path @@ -15,7 +24,8 @@ def parse_args(): parser = argparse.ArgumentParser( description="Create Pharmacophore Features from a Ligand Overlay" ) - parser.add_argument('-i', '--input_folder', type=str, help='Input file(s) path.') + parser.add_argument('-i', '--input_folder', type=str, required=True, + help='Input file(s) path.') parser.add_argument('-o', '--output_folder', type=str, default=None, help="Output folder path. Defaults to a 'queries' folder in the current directory.") parser.add_argument('-f', '--feature_definitions', type=str, required=True, diff --git a/scripts/discovery/overlay_to_pharmacophore/tests/tests.py b/scripts/discovery/overlay_to_pharmacophore/tests/tests.py index 2406ea5..1032c3c 100644 --- a/scripts/discovery/overlay_to_pharmacophore/tests/tests.py +++ b/scripts/discovery/overlay_to_pharmacophore/tests/tests.py @@ -20,7 +20,7 @@ def run_from_repo_root(monkeypatch): @pytest.fixture def overlay_data(tmp_path) -> OverlayData: - """Overlay/pharmacophore file locations for the bundled ``job1`` example.""" + """Overlay/pharmacophore file locations for the job1 example.""" return OverlayData( input_folder=JOB_FOLDER, output_folder=tmp_path, @@ -64,7 +64,7 @@ def test_projected_features_have_virtual_points(features): def test_unprojected_acceptors_have_no_virtual_points(overlay_data): - """With ``projected=False`` acceptors are kept as plain points.""" + """With `projected=False` acceptors are kept as single points.""" features = OverlayToPharmFeatures(overlay_data, projected=False).features acceptors = [f for f in features if f.label.startswith("acceptor")] assert acceptors, "No acceptor features extracted from pharmacophore." @@ -102,4 +102,3 @@ def test_missing_feature_definitions_raises(features, tmp_path): feature_definitions=tmp_path / "does_not_exist", output_file=tmp_path / "query.cm", ) - diff --git a/scripts/discovery/overlay_to_pharmacophore/write_query.py b/scripts/discovery/overlay_to_pharmacophore/write_query.py index 4d231b4..7148d76 100644 --- a/scripts/discovery/overlay_to_pharmacophore/write_query.py +++ b/scripts/discovery/overlay_to_pharmacophore/write_query.py @@ -12,6 +12,7 @@ class FeaturesToCrossMinerQuery: {} FEATURE_SUBSTRUCTURE_END """ + def __init__( self, pharm_feature_points: list[PharmFeaturePoint], @@ -58,7 +59,6 @@ def _feature_definition_file(self, feature: PharmFeaturePoint) -> Path: ) return matches[0] - @staticmethod def _create_feature_vector_text(feature: PharmFeaturePoint) -> str: """ From a4ba091180dbcffe0d160e787f227ff09cf39942 Mon Sep 17 00:00:00 2001 From: cringrose Date: Fri, 7 Aug 2026 10:43:45 +0100 Subject: [PATCH 15/15] NO_JIRA Added Pharmacophore generator script to create CrossMiner queries from Ligand Overlay outputs. Includes basic tests and examples. --- scripts/discovery/overlay_to_pharmacophore/cluster.py | 9 +++++++++ .../overlay_to_pharmacophore/crossminer_search.py | 9 +++++++++ .../discovery/overlay_to_pharmacophore/datastructures.py | 9 +++++++++ scripts/discovery/overlay_to_pharmacophore/main.py | 2 +- scripts/discovery/overlay_to_pharmacophore/overlay.py | 9 +++++++++ .../discovery/overlay_to_pharmacophore/tests/tests.py | 9 +++++++++ .../discovery/overlay_to_pharmacophore/write_query.py | 9 +++++++++ 7 files changed, 55 insertions(+), 1 deletion(-) diff --git a/scripts/discovery/overlay_to_pharmacophore/cluster.py b/scripts/discovery/overlay_to_pharmacophore/cluster.py index e49add2..d551569 100644 --- a/scripts/discovery/overlay_to_pharmacophore/cluster.py +++ b/scripts/discovery/overlay_to_pharmacophore/cluster.py @@ -1,3 +1,12 @@ +""" +This script can be used for any purpose without limitation subject to the +conditions at https://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx +This permission notice and the following statement of attribution must be +included in all copies or substantial portions of this script. + +"07/08/2026n": created by the Cambridge Crystallographic Data Centre +""" + from collections import defaultdict import numpy as np diff --git a/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py b/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py index c35ca78..f07cac3 100644 --- a/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py +++ b/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py @@ -1,3 +1,12 @@ +""" +This script can be used for any purpose without limitation subject to the +conditions at https://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx +This permission notice and the following statement of attribution must be +included in all copies or substantial portions of this script. + +"07/08/2026": created by the Cambridge Crystallographic Data Centre +""" + from pathlib import Path from ccdc.pharmacophore import Pharmacophore diff --git a/scripts/discovery/overlay_to_pharmacophore/datastructures.py b/scripts/discovery/overlay_to_pharmacophore/datastructures.py index 5f3bfe4..42d0ab7 100644 --- a/scripts/discovery/overlay_to_pharmacophore/datastructures.py +++ b/scripts/discovery/overlay_to_pharmacophore/datastructures.py @@ -1,3 +1,12 @@ +""" +This script can be used for any purpose without limitation subject to the +conditions at https://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx +This permission notice and the following statement of attribution must be +included in all copies or substantial portions of this script. + +"07/08/2026": created by the Cambridge Crystallographic Data Centre +""" + from __future__ import annotations from dataclasses import dataclass diff --git a/scripts/discovery/overlay_to_pharmacophore/main.py b/scripts/discovery/overlay_to_pharmacophore/main.py index 74c11dc..d2e32d5 100644 --- a/scripts/discovery/overlay_to_pharmacophore/main.py +++ b/scripts/discovery/overlay_to_pharmacophore/main.py @@ -4,7 +4,7 @@ This permission notice and the following statement of attribution must be included in all copies or substantial portions of this script. -"date of creation": created by the Cambridge Crystallographic Data Centre +"07/08/2026": created by the Cambridge Crystallographic Data Centre """ import argparse diff --git a/scripts/discovery/overlay_to_pharmacophore/overlay.py b/scripts/discovery/overlay_to_pharmacophore/overlay.py index 1bc51bd..a98a4d6 100644 --- a/scripts/discovery/overlay_to_pharmacophore/overlay.py +++ b/scripts/discovery/overlay_to_pharmacophore/overlay.py @@ -1,3 +1,12 @@ +""" +This script can be used for any purpose without limitation subject to the +conditions at https://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx +This permission notice and the following statement of attribution must be +included in all copies or substantial portions of this script. + +"07/08/2026": created by the Cambridge Crystallographic Data Centre +""" + from functools import cached_property import numpy as np diff --git a/scripts/discovery/overlay_to_pharmacophore/tests/tests.py b/scripts/discovery/overlay_to_pharmacophore/tests/tests.py index 1032c3c..06d82af 100644 --- a/scripts/discovery/overlay_to_pharmacophore/tests/tests.py +++ b/scripts/discovery/overlay_to_pharmacophore/tests/tests.py @@ -1,3 +1,12 @@ +""" +This script can be used for any purpose without limitation subject to the +conditions at https://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx +This permission notice and the following statement of attribution must be +included in all copies or substantial portions of this script. + +"07/08/2026": created by the Cambridge Crystallographic Data Centre +""" + from pathlib import Path import pytest diff --git a/scripts/discovery/overlay_to_pharmacophore/write_query.py b/scripts/discovery/overlay_to_pharmacophore/write_query.py index 7148d76..06ea724 100644 --- a/scripts/discovery/overlay_to_pharmacophore/write_query.py +++ b/scripts/discovery/overlay_to_pharmacophore/write_query.py @@ -1,3 +1,12 @@ +""" +This script can be used for any purpose without limitation subject to the +conditions at https://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx +This permission notice and the following statement of attribution must be +included in all copies or substantial portions of this script. + +"07/08/2026": created by the Cambridge Crystallographic Data Centre +""" + from pathlib import Path from datastructures import PharmFeaturePoint