NO_JIRA Added Pharmacophore generator script to create CrossMiner que…

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\<username>\CCDC\ccdc-software\csd-crossminer\feature_definitions`
+
+```
+python main.py -i <overlay_folder> -o <output_folder> -f <feature_definitions_folder>
+```
+
+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.

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)

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

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

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()