diff --git a/prody/proteins/_vorpy_aw.py b/prody/proteins/_vorpy_aw.py new file mode 100644 index 000000000..67c94fa21 --- /dev/null +++ b/prody/proteins/_vorpy_aw.py @@ -0,0 +1,644 @@ +# -*- coding: utf-8 -*- +"""Additively-weighted (Apollonius) Voronoi tessellation for :func:`.calcChannels` +``diagram="weighted"``, built on the third-party ``vorpy`` package. + +vorpy computes the *true* additively-weighted (a.k.a. Apollonius) Voronoi network, +which bakes van der Waals radii into the diagram exactly instead of approximating it +by homogenising atoms into uniform balls. This module: + +1. Ships a compiled (numba) drop-in for vorpy's ``calc_vert`` (the 4-sphere Apollonius + vertex solve). vorpy's own ``calc_vert`` solves a quadratic with ``numpy.roots`` + (a companion-matrix eigenvalue decomposition) and wraps every call in numpy + ``seterr`` context managers, costing ~116 us/call and ~80% of the total run time. + The compiled version uses the closed-form quadratic root and is bit-equivalent on + the common (non-degenerate) path, ~37x faster, deferring to the original only on + the rare singular (``F ~= 0``) branch. +2. Drives ``vorpy.Network`` in additively-weighted mode and adapts its vertex network + into the ``(simplices, neighbors, vertices)`` arrays the channel finder consumes, + plus an exact per-vertex clearance. + +The heavy imports (``vorpy``, ``numba``) happen at import time of *this* module, which +is imported lazily from the ``diagram="weighted"`` branch of :func:`.calcChannels`, so +they are not a hard dependency of ProDy or of the default channel path. +""" + +import os +import io +import re +import sys +import time +import hashlib +import contextlib +import tempfile +from math import sqrt + +import numpy as np + +__all__ = ['buildAwTessellation', 'accelerateVorpy', 'resolveCachePath'] + +# On-disk cache of the raw additively-weighted network. Building it with vorpy is the +# expensive step (minutes on a full structure); the downstream channel search is cheap +# and often re-run with different r2/bottleneck/sparsity/start_point, none of which +# change the diagram. Persisting the (simplices, neighbors, vertices, clearances) +# arrays lets those re-runs -- and repeated debugging passes -- skip vorpy entirely. +_CACHE_VERSION = 1 # bump whenever the cached array layout/semantics change + +_EPS_F = 1e-12 # |F| below this => degenerate spatial matrix, defer to vorpy +_TOL_IMAG = 1e-12 # matches vorpy _real_roots_quadratic imaginary tolerance + + +def _importNumbaKernel(): + """Build and return the numba-compiled AW vertex core. Kept in a function so the + numba import/compile only happens when the weighted path is actually used.""" + from numba import njit + + @njit(cache=True, fastmath=False) + def _awCore(b0x, b0y, b0z, b1x, b1y, b1z, b2x, b2y, b2z, b3x, b3y, b3z, + r0, r1, r2, r3): + # The four ball centres as 12 scalars followed by the four radii. Taking + # scalars (not (4,3)/(4,) arrays) lets calcVertFast skip building a numpy + # array per call -- np.ascontiguousarray was ~3 s of the accelerated run. + # Returns (degenerate, code, ax,ay,az,ar, bx,by,bz,br). + # degenerate=1 -> caller falls back to vorpy's original calc_vert (F ~= 0). + # code: 0 no vertex, 1 loc only, 2 loc+loc2 (ordered by |R| ascending). + l0x, l0y, l0z = b0x, b0y, b0z + l1x, l1y, l1z = b1x - b0x, b1y - b0y, b1z - b0z + l2x, l2y, l2z = b2x - b0x, b2y - b0y, b2z - b0z + l3x, l3y, l3z = b3x - b0x, b3y - b0y, b3z - b0z + r0_2 = r0 * r0 + + a1, b1, c1, d1 = 2*l1x, 2*l1y, 2*l1z, 2*(r1 - r0) + f1 = r0_2 - r1*r1 + l1x*l1x + l1y*l1y + l1z*l1z + a2, b2, c2, d2 = 2*l2x, 2*l2y, 2*l2z, 2*(r2 - r0) + f2 = r0_2 - r2*r2 + l2x*l2x + l2y*l2y + l2z*l2z + a3, b3, c3, d3 = 2*l3x, 2*l3y, 2*l3z, 2*(r3 - r0) + f3 = r0_2 - r3*r3 + l3x*l3x + l3y*l3y + l3z*l3z + + F = a1*b2*c3 - a1*b3*c2 - a2*b1*c3 + a2*b3*c1 + a3*b1*c2 - a3*b2*c1 + if F > -_EPS_F and F < _EPS_F: + return 1, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + F_2 = F * F + F10 = b1*c2*f3 - b1*c3*f2 - b2*c1*f3 + b2*c3*f1 + b3*c1*f2 - b3*c2*f1 + F11 = -b1*c2*d3 + b1*c3*d2 + b2*c1*d3 - b2*c3*d1 - b3*c1*d2 + b3*c2*d1 + F20 = -a1*c2*f3 + a1*c3*f2 + a2*c1*f3 - a2*c3*f1 - a3*c1*f2 + a3*c2*f1 + F21 = a1*c2*d3 - a1*c3*d2 - a2*c1*d3 + a2*c3*d1 + a3*c1*d2 - a3*c2*d1 + F30 = a1*b2*f3 - a1*b3*f2 - a2*b1*f3 + a2*b3*f1 + a3*b1*f2 - a3*b2*f1 + F31 = -a1*b2*d3 + a1*b3*d2 + a2*b1*d3 - a2*b3*d1 - a3*b1*d2 + a3*b2*d1 + + qa = (F11*F11 + F21*F21 + F31*F31) - F_2 + qb = 2.0 * ((F10*F11 + F20*F21 + F30*F31) - r0*F_2) + qc = (F10*F10 + F20*F20 + F30*F30) - r0_2*F_2 + + R_a = 0.0; R_b = 0.0; nR = 0 + if qa > -1e-300 and qa < 1e-300: + if not (qb > -1e-300 and qb < 1e-300): + R_a = -qc / qb; nR = 1 + else: + disc = qb*qb - 4.0*qa*qc + if disc >= 0.0: + sq = np.sqrt(disc) + R_a = (-qb + sq) / (2.0*qa) + R_b = (-qb - sq) / (2.0*qa) + nR = 2 + elif disc > -_TOL_IMAG: + R_a = -qb / (2.0*qa); nR = 1 + + if nR == 0: + return 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + if nR == 1: + x = (F10 + R_a*F11)/F + l0x + y = (F20 + R_a*F21)/F + l0y + z = (F30 + R_a*F31)/F + l0z + return 0, 1, x, y, z, R_a, 0.0, 0.0, 0.0, 0.0 + + if abs(R_a) > abs(R_b): + R_a, R_b = R_b, R_a + x0 = (F10 + R_a*F11)/F + l0x + y0 = (F20 + R_a*F21)/F + l0y + z0 = (F30 + R_a*F31)/F + l0z + x1 = (F10 + R_b*F11)/F + l0x + y1 = (F20 + R_b*F21)/F + l0y + z1 = (F30 + R_b*F31)/F + l0z + max_ball_rad = max(max(r0, r1), max(r2, r3)) + + if R_a < 0.0 or R_b < 0.0: + if R_a > 0.0 or abs(R_a) < max_ball_rad: + if R_b > 0.0 or abs(R_b) < max_ball_rad: + return 0, 2, x0, y0, z0, R_a, x1, y1, z1, R_b + return 0, 1, x0, y0, z0, R_a, 0.0, 0.0, 0.0, 0.0 + elif R_b > 0.0 or abs(R_b) < max_ball_rad: + return 0, 1, x1, y1, z1, R_b, 0.0, 0.0, 0.0, 0.0 + return 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + return 0, 2, x0, y0, z0, R_a, x1, y1, z1, R_b + + return _awCore + + +_ACCEL_APPLIED = False + + +def _patchEverywhere(name, orig, replacement): + """Rebind ``name`` to ``replacement`` in every already-imported vorpy module that + currently holds ``orig`` (vorpy does ``from ... import name`` widely, so each such + binding must be replaced). Returns the number of bindings swapped.""" + n = 0 + for modname, mod in list(sys.modules.items()): + if mod is None or not modname.startswith('vorpy'): + continue + if getattr(mod, name, None) is orig: + setattr(mod, name, replacement) + n += 1 + return n + + +def accelerateVorpy(): + """Monkeypatch vorpy's two pure-Python hot functions with faster drop-ins, + everywhere they are bound: ``calc_vert`` (the 4-sphere Apollonius solve, replaced + by the compiled closed-form kernel) and ``calc_dist`` (a 3-D Euclidean distance + that vorpy implements with per-call numpy array allocation + a Python ``sum``, + called millions of times). Idempotent. Returns the number of bindings replaced.""" + global _ACCEL_APPLIED + if _ACCEL_APPLIED: + return 0 + + # calc_vert -> compiled closed-form kernel + import vorpy.src.calculations.vert as _vmod + vert_orig = _vmod.calc_vert + core = _importNumbaKernel() + + def calcVertFast(locs, rads): + # locs is a 4-list of ball centres, rads a 4-list of radii (vorpy always + # calls with [locs[i] for i in balls]). Hand the kernel their 16 scalars + # directly instead of materialising a numpy array per call. + l0, l1, l2, l3 = locs + deg, code, ax, ay, az, ar, bx, by, bz, br = core( + l0[0], l0[1], l0[2], l1[0], l1[1], l1[2], + l2[0], l2[1], l2[2], l3[0], l3[1], l3[2], + rads[0], rads[1], rads[2], rads[3]) + if deg: + return vert_orig(locs, rads) + if code == 0: + return None, None, None, None + if code == 1: + return [ax, ay, az], ar, None, None + return [ax, ay, az], ar, [bx, by, bz], br + + # calc_dist -> allocation-free 3-D distance (falls back for non-3-D inputs) + import vorpy.src.calculations.calcs as _cmod + dist_orig = _cmod.calc_dist + + def calcDistFast(l0, l1): + if len(l0) == 3: + d0 = l0[0] - l1[0] + d1 = l0[1] - l1[1] + d2 = l0[2] - l1[2] + return sqrt(d0 * d0 + d1 * d1 + d2 * d2) + return dist_orig(l0, l1) + + # write_verts -> no-op: vorpy dumps every vertex to a .txt at the end of + # find_net_verts via pandas iterrows (slow over tens of thousands of rows); we + # build our own arrays from net.verts and never read that file. + import vorpy.src.output as _omod + wv_orig = _omod.write_verts + + def writeVertsNoop(net): + return None + + n = _patchEverywhere('calc_vert', vert_orig, calcVertFast) + n += _patchEverywhere('calc_dist', dist_orig, calcDistFast) + n += _patchEverywhere('write_verts', wv_orig, writeVertsNoop) + _ACCEL_APPLIED = True + return n + + +# vorpy prints a "\r...finding vertices: verts - %" line to stdout as it +# traces the network; we parse it to drive a ProDy progress bar and otherwise swallow +# vorpy's chatter. +_PROG_RE = re.compile(r'finding vertices:\s*(\d+)\s*verts\s*-\s*([\d.]+)\s*%') + + +class _ProgressSink(io.StringIO): + """stdout replacement that discards vorpy's output but extracts its vertex-finding + percentage and forwards it to ``on_progress(n_verts, percent)``.""" + + def __init__(self, on_progress): + super().__init__() + self._cb = on_progress + + def write(self, s): + if '%' not in s: # skip the regex on vorpy's non-progress chatter + return len(s) + m = _PROG_RE.search(s) + if m is not None: + try: + self._cb(int(m.group(1)), float(m.group(2))) + except Exception: + pass + return len(s) + + +@contextlib.contextmanager +def _captureVorpy(on_progress=None): + """Contain vorpy's on-disk vertex writes in a throwaway temp directory (cwd + restored afterwards) and capture its stdout: parsed for progress when + ``on_progress`` is given, silently discarded otherwise.""" + cwd = os.getcwd() + tmp = tempfile.mkdtemp(prefix='prody_awvor_') + sink = _ProgressSink(on_progress) if on_progress is not None else io.StringIO() + try: + os.chdir(tmp) + with contextlib.redirect_stdout(sink): + yield + finally: + os.chdir(cwd) + + +def _cacheKey(coords, vdw_radii, max_vert): + """Content hash identifying a tessellation input: the ball centres, their radii + and the ``max_vert`` cutoff (the only inputs that change the diagram). A cache is + reused only when all three match, so editing the structure -- or, via ``r1``, + ``max_vert`` -- transparently forces a recompute.""" + h = hashlib.sha1() + h.update(np.ascontiguousarray(coords, dtype=np.float64).tobytes()) + h.update(np.ascontiguousarray(vdw_radii, dtype=np.float64).tobytes()) + h.update(np.float64(max_vert).tobytes()) + return 'v{0}:{1}'.format(_CACHE_VERSION, h.hexdigest()) + + +def _loadCache(path, key): + """Return the cached ``(simplices, neighbors, vertices, clearances)`` at ``path`` + if it exists and was built for ``key``, else ``None`` (caller recomputes). Any + read/format error is swallowed and treated as a miss.""" + if not path or not os.path.isfile(path): + return None + try: + with np.load(path, allow_pickle=False) as d: + if str(d['key']) != key: + return None + return (d['simplices'], d['neighbors'], d['vertices'], d['clearances']) + except Exception: + return None + + +def _saveCache(path, key, result): + """Persist a tessellation ``result`` under ``key`` to ``path`` (compressed npz), + writing to a temporary file and renaming so a crashed write can't leave a + half-written cache. Failures warn but do not abort the calculation.""" + simplices, neighbors, vertices, clearances = result + tmp = path + '.tmp.npz' + try: + d = os.path.dirname(path) + if d and not os.path.isdir(d): + os.makedirs(d) + np.savez_compressed(tmp, key=np.array(key), simplices=simplices, + neighbors=neighbors, vertices=vertices, + clearances=clearances) + os.replace(tmp, path) + except Exception as exc: + from prody import LOGGER + LOGGER.warn('Could not write additively-weighted Voronoi cache to {0}: {1}' + .format(path, exc)) + try: + if os.path.isfile(tmp): + os.remove(tmp) + except OSError: + pass + + +def resolveCachePath(cache, output_path=None, title=None): + """Resolve the ``weighted_cache`` argument of :func:`.calcChannels` to a concrete + ``.npz`` path or ``None``. ``False``/``None`` (and ``""``) disable caching; a + string is used verbatim; ``True`` auto-derives a path next to ``output_path`` + (e.g. ``aw.pqr`` -> ``aw.awvoronoi.npz``) or, lacking an ``output_path``, from the + structure ``title`` in the current working directory.""" + if not cache: + return None + if isinstance(cache, str): + return cache + if output_path: + return os.path.splitext(output_path)[0] + '.awvoronoi.npz' + title = (title or 'structure').replace(' ', '_') + return os.path.join(os.getcwd(), title + '.awvoronoi.npz') + + +def buildAwTessellation(coords, vdw_radii, max_vert=8.0, accelerate=True, + target_atoms=500, cache=None): + """Build the additively-weighted (Apollonius) Voronoi network of the balls + ``(coords, vdw_radii)`` with vorpy and adapt it to the channel finder's arrays. + + Large *and spatially extended* structures are split into local boxes (each padded + with a correctness halo), tessellated one after another and merged by global 4-ball + signature. Splitting only engages where a box core would be comfortably larger than + its halo (see ``_planBoxes``); a compact globular structure -- even a few thousand + atoms -- stays on the single-pass path, since there the halo would re-cover the + whole structure and decomposition would be a net loss. + + :arg coords: (N,3) ball centres (atom coordinates), in Angstrom. + :arg vdw_radii: (N,) van der Waals radii aligned with ``coords``. + :arg max_vert: probe distance / maximum empty-sphere radius. Vertices whose + clearance exceeds this are dropped. ~8 A keeps every channel-relevant vertex + while pruning irrelevant large voids. Default 8.0. + :arg accelerate: apply the compiled ``calc_vert`` kernel first. Default True. + :arg target_atoms: approximate number of balls per box core; also the threshold + below which the structure is tessellated in a single pass. Default 500. + :arg cache: optional path to an ``.npz`` cache file. When given and it already + holds a diagram built for this exact ``(coords, vdw_radii, max_vert)``, it is + loaded and vorpy is skipped; otherwise the freshly built diagram is written + there for next time. A mismatching or unreadable file is ignored and + overwritten. ``None`` (default) disables caching. + + :returns: ``(simplices, neighbors, vertices, clearances)`` where + ``simplices`` is (M,4) int (the 4 tangent atoms per vertex), + ``neighbors`` is (M,4) int (face-adjacent vertices sharing 3 atoms, -1 padded), + ``vertices`` is (M,3) float (Voronoi-vertex positions), and + ``clearances`` is (M,) float (exact additively-weighted clearance per vertex). + """ + coords = np.asarray(coords, dtype=float) + vdw_radii = np.asarray(vdw_radii, dtype=float) + N = len(coords) + from prody import LOGGER + + # Reuse a previously saved diagram for this exact input if one is cached. + key = _cacheKey(coords, vdw_radii, max_vert) + if cache: + cached = _loadCache(cache, key) + if cached is not None: + LOGGER.info('Reusing additively-weighted Voronoi diagram cached at {0} ' + '({1} vertices); skipping vorpy.'.format(cache, len(cached[2]))) + return cached + + # Each box must be padded with a correctness halo (see _buildAwDecomposed). A box + # is only worth splitting out if its core is meaningfully larger than that halo, + # otherwise core+halo re-tessellates almost the whole structure at a net loss. So + # decomposition only engages for large/elongated systems; a compact globular + # protein (even a few thousand atoms) stays on the single-pass path. + halo = _haloWidth(max_vert, vdw_radii) + boxes = _planBoxes(coords, target_atoms, halo) if N > target_atoms else None + if boxes and len(boxes) > 1: + # Only decompose if it meaningfully shrinks the per-box work. Every box is + # padded with the halo, so for a compact (globular) structure each box's + # core+halo re-covers most of the atoms and splitting is a net loss (extra + # overlap work, no gain). In that case fall back to the single pass, which is + # both faster here and shows a fine-grained progress bar. + loads = _boxLoads(coords, boxes, halo) + worst = max(ch for _, ch in loads) + if worst > 0.6 * N: + LOGGER.info('Weighted tessellation: skipping decomposition (largest of {0} ' + 'boxes would still span {1}/{2} atoms once haloed); using a ' + 'single pass.'.format(len(boxes), worst, N)) + boxes = None + if not boxes or len(boxes) == 1: + if accelerate: + accelerateVorpy() + prog = _MonoProgress() + verts = _runVorpy(coords, vdw_radii, max_vert, accelerate, on_progress=prog) + prog.finish() + balls = list(verts['balls']) + result = _assemble(balls, list(verts['loc']), list(verts['rad']), + list(verts['dub']) if 'dub' in verts else None) + else: + result = _buildAwDecomposed(coords, vdw_radii, max_vert, accelerate, + boxes, halo) + + if cache: + _saveCache(cache, key, result) + LOGGER.info('Cached additively-weighted Voronoi diagram ({0} vertices) to {1}.' + .format(len(result[2]), cache)) + return result + + +def _haloWidth(max_vert, vdw_radii): + """Halo padding a box needs so that every vertex it owns (clearance <= max_vert) + has all its defining/blocking balls present: clearance + ball radius + margin.""" + return float(max_vert) + float(np.max(vdw_radii)) + 1.0 + + +def _boxLoads(coords, boxes, halo): + """For each box return ``(core_count, core_plus_halo_count)``: how many atoms it + owns and how many it must actually tessellate (core plus halo).""" + loads = [] + for lo, hi in boxes: + core = int(np.all((coords >= lo) & (coords < hi), axis=1).sum()) + near = int(np.all((coords >= lo - halo) & (coords < hi + halo), axis=1).sum()) + loads.append((core, near)) + return loads + + +class _MonoProgress(object): + """Throttled adapter turning vorpy's vertex-finding percentage into periodic ProDy + log lines for the single-pass (mono) path. Uses plain ``LOGGER.info`` (not the + ``progress``/``update`` bar) so it neither assumes a TTY nor mutates the logger's + verbosity state, and so it can never suppress later messages.""" + + _msg = 'Additively-weighted Voronoi diagram' + _interval = 3.0 # seconds between log lines + + def __init__(self): + from prody import LOGGER + self._logger = LOGGER + self._last_t = 0.0 + + def __call__(self, n_verts, percent): + # Throttle by wall time (vorpy calls this per edge). Report the vertex COUNT, + # which is exact and monotonic; vorpy's percentage is an estimate that saturates + # at 100% while work continues, so it is only a secondary hint. Time-based + # throttling keeps the line moving even past that bogus 100%. + now = time.perf_counter() + if now - self._last_t >= self._interval: + self._last_t = now + self._logger.info('{0}: {1} vertices ({2}%)'.format( + self._msg, n_verts, min(percent, 100))) + + def finish(self): + # Nothing to restore: LOGGER.info leaves no state. The caller's LOGGER.report + # ('... tessellation constructed in %.2fs') closes out the stage. + pass + + +def _runVorpy(coords, vdw_radii, max_vert, accelerate, on_progress=None): + """Run vorpy's additively-weighted vertex finder on the given balls and return + its ``verts`` DataFrame (columns: balls, loc, rad, dub). When ``on_progress`` is + given it is called with ``(n_verts, percent)`` as vorpy traces the network.""" + if accelerate: + accelerateVorpy() + + from vorpy.src.network import Network + + coords = np.asarray(coords, dtype=float) + vdw_radii = np.asarray(vdw_radii, dtype=float) + settings = {'surf_res': 0.2, 'surf_col': 'plasma', 'surf_scheme': 'mean', + 'max_vert': float(max_vert), 'box_size': 1.5, 'net_type': 'aw', + 'build_type': 'all', 'num_splits': None, 'print_metrics': False, + 'ball_type': 'mol', 'sys_dir': os.getcwd(), 'foam_box': None, + 'atom_rad': None, 'scheme_factor': 'log'} + locs = [row.copy() for row in coords] + rads = [float(x) for x in vdw_radii] + + with _captureVorpy(on_progress): + net = Network(locs=locs, rads=rads, settings=settings, sort_balls=True) + net.find_verts() + + return net.verts + + +_BIG = 1.0e6 # sentinel bound: outer leaf faces extend to +/- this so that every + # Voronoi vertex (including those outside the atom bounding box) is + # owned by exactly one box core. +_SPLIT_HALO_FACTOR = 5.0 # only split an axis whose extent exceeds this * halo, so + # each child core is comfortably wider than its halo shell. + + +def _planBoxes(coords, target, halo): + """Partition space into axis-aligned leaf boxes each owning <= ``target`` centres, + by recursive median bisection along the widest data axis. An axis is only split + when its extent exceeds ``_SPLIT_HALO_FACTOR * halo``; recursion stops otherwise, + so compact structures collapse to a single box (and take the single-pass path). + + Returns a list of ``(lo, hi)`` half-open bound pairs (float (3,) arrays) that tile + *all* of space: internal cut planes are finite, while any leaf face on the global + boundary is pushed out to +/-``_BIG`` so exterior vertices are still owned. Every + centre falls in exactly one box.""" + lo0 = coords.min(axis=0) + hi0 = coords.max(axis=0) + min_extent = _SPLIT_HALO_FACTOR * halo + boxes = [] + + def rec(idx, lo, hi, lo_out, hi_out): + # Effective bounds: outer faces open out to the +/-_BIG sentinel. + eff_lo = np.where(lo_out, -_BIG, lo) + eff_hi = np.where(hi_out, _BIG, hi) + extent = hi - lo + # Widest axis that is both long enough to be worth splitting and actually + # spans the points in this cell. + order = np.argsort(extent)[::-1] + axis = None + for a in order: + if extent[a] > min_extent: + axis = int(a) + break + if len(idx) <= target or axis is None: + boxes.append((eff_lo, eff_hi)) + return + split = float(np.median(coords[idx, axis])) + left = idx[coords[idx, axis] < split] + right = idx[coords[idx, axis] >= split] + # guard against a degenerate median that fails to divide the points + if split <= lo[axis] or split >= hi[axis] or len(left) == 0 or len(right) == 0: + boxes.append((eff_lo, eff_hi)) + return + hl = hi.copy(); hl[axis] = split + lr = lo.copy(); lr[axis] = split + lo_out_r = lo_out.copy(); lo_out_r[axis] = False + hi_out_l = hi_out.copy(); hi_out_l[axis] = False + rec(left, lo, hl, lo_out, hi_out_l) + rec(right, lr, hi, lo_out_r, hi_out) + + ones = np.ones(3, dtype=bool) + rec(np.arange(len(coords)), lo0, hi0, ones.copy(), ones.copy()) + return boxes + + +def _buildAwDecomposed(coords, vdw_radii, max_vert, accelerate, boxes, halo): + """Tessellate each box (core + halo) independently, keep only the vertices whose + location falls in the box core, and merge them by global 4-ball signature. + + The halo width is chosen so a ball can invalidate a kept vertex (clearance <= + max_vert) only if its centre lies within clearance + its radius <= max_vert + + max_radius of the vertex; since the vertex sits in the core, that halo around the + core contains every ball a core vertex depends on (+1 A margin).""" + from prody import LOGGER + + tasks, cores, haloed = [], [], [] + for lo, hi in boxes: + core = np.all((coords >= lo) & (coords < hi), axis=1) + near = np.all((coords >= lo - halo) & (coords < hi + halo), axis=1) + if not core.any(): + continue + gidx = np.nonzero(near)[0] + tasks.append((coords[gidx].copy(), vdw_radii[gidx].copy(), gidx, + lo.copy(), hi.copy(), float(max_vert), bool(accelerate))) + cores.append(int(core.sum())) + haloed.append(len(gidx)) + + LOGGER.info('Weighted tessellation: {0} boxes, {1}-{2} atoms each ({3}-{4} with ' + 'halo).'.format(len(tasks), min(cores), max(cores), + min(haloed), max(haloed))) + + results = [] + for t in tasks: + results.append(_awBoxTask(t)) + LOGGER.info('Additively-weighted Voronoi diagram: box {0}/{1} done' + .format(len(results), len(tasks))) + + # Merge box outputs, deduplicating by global 4-ball signature. Ownership-by-core- + # location already assigns each vertex to a single box, so collisions are rare and + # identical; keep the first seen. + seen = {} + balls, locs, rads = [], [], [] + for out in results: + for sig, x, y, z, rad in out: + if sig in seen: + continue + seen[sig] = True + balls.append(sig) + locs.append((x, y, z)) + rads.append(rad) + + return _assemble(balls, locs, rads, None) + + +def _awBoxTask(task): + """Tessellate one box's balls and return the core-owned vertices as + ``(global_signature, x, y, z, clearance)`` records.""" + coords, vdw_radii, gidx, lo, hi, max_vert, accelerate = task + verts = _runVorpy(coords, vdw_radii, max_vert, accelerate) + out = [] + for bl, loc, rad in zip(verts['balls'], verts['loc'], verts['rad']): + x, y, z = float(loc[0]), float(loc[1]), float(loc[2]) + if not (lo[0] <= x < hi[0] and lo[1] <= y < hi[1] and lo[2] <= z < hi[2]): + continue + sig = tuple(sorted(int(gidx[int(b)]) for b in bl)) + out.append((sig, x, y, z, float(rad))) + return out + + +def _assemble(balls, locs, rads, dubs=None): + """Build ``(simplices, neighbors, vertices, clearances)`` from per-vertex lists of + 4-ball tuples, locations and additively-weighted clearances. ``balls`` indices must + already be in the final (global) numbering.""" + from prody import LOGGER + + if dubs is not None: + ndub = int(sum(1 for d in dubs if d != 0)) + if ndub: + LOGGER.warn('vorpy AW network has {0} doublet vertices; their face ' + 'adjacency may be incomplete.'.format(ndub)) + + M = len(balls) + simplices = np.array([sorted(int(b) for b in bl) for bl in balls], + dtype=np.int64).reshape(M, 4) if M else \ + np.zeros((0, 4), dtype=np.int64) + vertices = np.array([np.asarray(l, dtype=float) for l in locs], + dtype=float).reshape(M, 3) if M else np.zeros((0, 3), float) + clearances = np.array([float(r) for r in rads], dtype=float) + + # Face adjacency: two vertices are neighbours iff they share 3 of their 4 atoms. + # Map each triangular face (sorted 3-atom tuple) -> the vertices carrying it; + # a face shared by exactly two vertices links them. + face_to_verts = {} + for i in range(M): + b0, b1, b2, b3 = simplices[i] + for face in ((b0, b1, b2), (b0, b1, b3), (b0, b2, b3), (b1, b2, b3)): + face_to_verts.setdefault(face, []).append(i) + + neighbors = np.full((M, 4), -1, dtype=np.int64) + for i in range(M): + b0, b1, b2, b3 = simplices[i] + for k, face in enumerate(((b0, b1, b2), (b0, b1, b3), + (b0, b2, b3), (b1, b2, b3))): + owners = face_to_verts[face] + if len(owners) == 2: + neighbors[i, k] = owners[0] if owners[1] == i else owners[1] + # len 1 -> boundary face (-1); len >2 -> degenerate, left as -1 + + return simplices, neighbors, vertices, clearances diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 3b9fde515..ba1435ef0 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -1,38 +1,89 @@ # -*- coding: utf-8 -*- -"""This module called CaviFinder and defines functions for calculating channels, tunnels and pores -within protein structure. +"""This module is called CaviTracer and defines functions for calculating +channels, tunnels, and surface cavities within protein structure. """ -__author__ = 'Karolina Mikulska-Ruminska', 'Eryk Trzcinski' -__credits__ = ['Karolina Mikulska-Ruminska', 'Eryk Trzcinski'] +__author__ = 'Karolina Mikulska-Ruminska', 'Jan Brezovsky', 'Eryk Trzcinski' +__credits__ = ['Karolina Mikulska-Ruminska', 'Jan Brezovsky', 'Eryk Trzcinski'] __email__ = ['karolamik@fizyka.umk.pl'] +import logging +from contextlib import contextmanager + import numpy as np -from numpy import * from prody import LOGGER, PY3K from prody.atomic import Atomic from prody.utilities import checkCoords, getCoords, isListLike from prody.proteins import writePDB, parsePDB, parsePQR from prody.ensemble import Ensemble -from prody.measure import calcCenter - -import multiprocessing -from .fixer import * -from .compare import * -from prody.measure import calcTransformation, calcDistance, calcRMSD, superpose +from prody.measure import calcCenter, calcTransformation, calcDistance, calcRMSD, superpose __all__ = ['getVmdModel', 'calcChannels', 'calcChannelsMultipleFrames', - 'getChannelParameters', 'getChannelAtoms', 'showChannels', 'showCavities', - 'showSurfaceCavities', 'selectChannelBySelection', 'getChannelResidueNames', + 'getChannelParameters', 'getChannelAtoms', 'showChannels', + 'showCavities', 'showSurfaceCavities', 'selectChannelBySelection', + 'getChannelResidueNames', 'calcChannelSurfaceOverlaps', 'calcSurfaceCavities', 'calcSurfaceCavitiesMultipleFrames', 'getSurfaceCavityParameters', 'getSurfaceCavityResidueNames', 'selectSurfaceCavityBySelection', - 'calcSurfaceCavityOverlaps','getSurfaceCavityResidueNamesMultipleFrames', - 'getSurfaceCavityParametersMultipleFrames', 'getChannelParametersMultipleFrames', + 'calcSurfaceCavityOverlaps', + 'getSurfaceCavityResidueNamesMultipleFrames', + 'getSurfaceCavityParametersMultipleFrames', + 'getChannelParametersMultipleFrames', '_reportAtomsInputComposition', 'getChannelResidueNamesMultipleFrames'] +# Sampling of the enclosure test used to strip the moat (see +# ChannelCalculator.calcEnclosure). These are constants, not knobs: the enclosure +# of a point depends on how many directions are sampled and how far they are +# followed, so min_enclosure is only meaningful against a fixed sampling. Adding +# rays lowers every enclosure, since more directions find more of the thin ways +# out of a channel, and so invalidates the threshold rather than refining it. +ENCLOSURE_RAYS = 32 +ENCLOSURE_RANGE = 25.0 +ENCLOSURE_STEP = 0.75 +# One radius for every atom, on the scale of a heavy-atom vdW radius. Enclosure is +# a burial heuristic, so resolving 1.52 A from 1.7 A would only shift every value +# by a little and be absorbed by min_enclosure; a single radius means a single +# tree and a plain nearest-neighbour test. +ENCLOSURE_RADIUS = 1.7 + +_OVERLAP_OFFSET_CACHE = {} + +@contextmanager +def _warningsDelivered(): + """Let WARNING records through for the duration of the block. + + Importing ProDy installs a logging filter that drops every WARNING record from + the package logger (prody.dynamics.adaptive2, at module scope), so every + ``LOGGER.warn`` in ProDy is silently discarded. Whether that filter should exist + at all is a question for the package as a whole; until it is settled, this + module at least delivers its own warnings. + + Offending filters are found by asking them, rather than by importing the class + and matching on type: any filter that rejects a synthetic WARNING record is + detached for the block and reinstated afterwards. That keeps this working if the + filter is renamed or moved, and it leaves the filter in force everywhere else, + so nothing outside this module changes behaviour.""" + logger = LOGGER._logger + probe = logging.LogRecord(logger.name, logging.WARNING, __file__, 0, + '', (), None) + muting = [f for f in list(logger.filters) + if not (f.filter(probe) if hasattr(f, 'filter') else f(probe))] + for f in muting: + logger.removeFilter(f) + try: + yield + finally: + for f in muting: + logger.addFilter(f) + + +def _warn(message): + """``LOGGER.warn``, but actually emitted. See :func:`_warningsDelivered`.""" + with _warningsDelivered(): + LOGGER.warn(message) + def checkAndImport(package_name): """Check for package and import it if possible and return **True**. @@ -51,43 +102,156 @@ def checkAndImport(package_name): if PY3K: import importlib.util if importlib.util.find_spec(package_name) is None: - LOGGER.warn("Package " + str(package_name) + " is not installed. Please install it to use this function.") + _warn("Package " + str(package_name) + " is not installed. " + "Please install it to use this function.") return False else: try: __import__(package_name) except ImportError: - LOGGER.warn("Package " + str(package_name) + " is not installed. Please install it to use this function.") + _warn("Package " + str(package_name) + " is not installed. " + "Please install it to use this function.") return False return True -def getVmdModel(vmd_path, atoms, representation='NewCartoon'): - """Generates a 3D model of molecular structures using VMD and returns it as an Open3D TriangleMesh. +def _getOverlapSphereOffsets(radius, resolution): + """Return integer voxel offsets inside a sphere of a given radius.""" + + key = (round(float(radius), 3), float(resolution)) + + if key in _OVERLAP_OFFSET_CACHE: + return _OVERLAP_OFFSET_CACHE[key] + + n = int(np.ceil(radius / resolution)) + grid = np.arange(-n, n + 1, dtype=int) + + dx, dy, dz = np.meshgrid(grid, grid, grid, indexing='ij') + offsets = np.vstack((dx.ravel(), dy.ravel(), dz.ravel())).T + xyz = offsets.astype(float) * resolution + mask = np.sum(xyz * xyz, axis=1) <= radius * radius + offsets = offsets[mask] + _OVERLAP_OFFSET_CACHE[key] = offsets + + return offsets + - This function creates a temporary PDB file from the provided atomic data and uses VMD (Visual Molecular Dynamics) - to render this data into an STL file, which is then loaded into Open3D as a TriangleMesh. The function handles - the creation and cleanup of temporary files and manages the subprocess call to VMD. +def _surfaceFromPqrWorker(args): + """Create voxelized FIL surface for one PQR file.""" + + pqr_file, resolution = args + atoms = parsePQR(pqr_file) + fil = atoms.select('resname FIL') + + if fil is None: + return set() + + coords = fil.getCoords() + radii = fil.getRadii() + surface = set() + + for center, radius in zip(coords, radii): + center_idx = np.rint(center / resolution).astype(int) + offsets = _getOverlapSphereOffsets(float(radius), resolution) + voxels = offsets + center_idx + surface.update(map(tuple, voxels)) + + return surface + + +def _reportAtomsInputComposition(atoms): + """Report the composition of atoms supplied for channel analysis. + + This function checks whether the input atomic structure contains only + protein atoms or also includes water, non-water HETATM records, or other + non-protein components. If non-protein atoms are present, a warning is + issued indicating that all supplied atoms will be included in the channel + calculation. + + The function does not modify or filter the input structure. To analyze only + the protein, the user should provide an appropriate ProDy selection, for + example ``atoms.select('protein')``. """ + + if not isinstance(atoms, Atomic): + raise TypeError( + "atoms must be a ProDy Atomic object, such as an AtomGroup " + "or Selection") + + protein = atoms.select('protein') + water = atoms.select('water') + hetero = atoms.select('hetero and not water') + other = atoms.select('not protein and not hetero') + nonprotein = atoms.select('not protein') + + if nonprotein is None: + LOGGER.info("The atoms supplied to calcChannels contain protein atoms only.") + return + + components = [] + + if water is not None: + components.append( + "water: {0} atoms in {1} residues".format( + water.numAtoms(), + len(np.unique(water.getResindices())))) + + if hetero is not None: + components.append( + "non-water hetero components: {0} atoms " + "(resnames: {1})".format( + hetero.numAtoms(), + ", ".join(sorted(np.unique(hetero.getResnames()))))) + + if other is not None: + components.append( + "other non-protein components: {0} atoms " + "(resnames: {1})".format( + other.numAtoms(), + ", ".join(sorted(np.unique(other.getResnames()))))) + + _warn("The atoms supplied to calcChannels() contain non-protein components: " + "{0}. All supplied atoms except waters will be used for channel analysis. " + "To analyze only the protein structure, provide an appropriate " + "selection, for example atoms.select('protein').".format( + "; ".join(components))) + + +def getVmdModel(vmd_path, atoms, representation='NewCartoon'): + """Generates a 3D model of molecular structures using VMD and returns + it as an Open3D TriangleMesh. + + This function creates a temporary PDB file from the provided atomic data + and uses VMD (Visual Molecular Dynamics) to render this data into an STL + file, which is then loaded into Open3D as a TriangleMesh. The function + handles the creation and cleanup of temporary files and manages the + subprocess call to VMD. To install Open3D use: - conda install open3d (for Anaconda users; version open3d-0.19.0 was used during the developement) - or pip install open3d + conda install open3d (for Anaconda users; version open3d-0.19.0 was used + during the development) or pip install open3d + + If problem with `ipykernel.comm.Comm` class appeared while using getVmdModel + please update dash: python -m pip install -U dash - :param vmd_path: Path to the VMD executable. This is required to run VMD and execute the TCL script. + :arg vmd_path: Path to the VMD executable. This is required to run VMD and + execute the TCL script. :type vmd_path: str - :param atoms: Atomic data to be written to a PDB file. This should be an object or data structure - that is compatible with the `writePDB` function. + :arg atoms: Atomic data to be written to a PDB file. This should be an + object or data structure that is compatible with the `writePDB` function. :type atoms: object - :raises ImportError: If required libraries ('subprocess', 'pathlib', 'tempfile', 'open3d') are not installed, - an ImportError is raised, specifying which libraries are missing. + :raises ImportError: If required libraries ('subprocess', 'pathlib', + 'tempfile', 'open3d') are not installed, an ImportError is raised, + specifying which libraries are missing. - :raises ValueError: If the STL file is not created or is empty, or if the STL file cannot be read as a TriangleMesh, + :raises ValueError: If the STL file is not created or is empty, or if the + STL file cannot be read as a TriangleMesh, a ValueError is raised. - :returns: An Open3D TriangleMesh object representing the 3D model generated from the PDB data. + :returns: An Open3D TriangleMesh object representing the 3D model generated + from the PDB data. :rtype: open3d.geometry.TriangleMesh Example usage: @@ -100,7 +264,8 @@ def getVmdModel(vmd_path, atoms, representation='NewCartoon'): if not checkAndImport(name): missing.append(name) if errorMsg is None: - errorMsg = 'To run getVmdModel, please install {0}'.format(missing[0]) + errorMsg = 'To run getVmdModel, ' \ + 'please install {0}'.format(missing[0]) else: errorMsg += ', ' + name @@ -125,7 +290,8 @@ def getVmdModel(vmd_path, atoms, representation='NewCartoon'): rep_key = representation.lower() if rep_key not in representation_map: raise ValueError( - "representation must be one of: 'NewCartoon', 'VDW', 'Surf', 'QuickSurf', or 'CPK'") + "representation must be one of: 'NewCartoon', 'VDW', 'Surf', " \ + "'QuickSurf', or 'CPK'") representation_style = representation_map[rep_key] if PY3K: @@ -143,7 +309,8 @@ def getVmdModel(vmd_path, atoms, representation='NewCartoon'): if PY3K: output_path = temp_script_path.parent / "output.stl" else: - output_path = os.path.join(os.path.dirname(temp_script.name), "output.stl") + output_path = os.path.join(os.path.dirname(temp_script.name), + "output.stl") vmd_script = """ set file_path [lindex $argv 0] @@ -165,7 +332,8 @@ def getVmdModel(vmd_path, atoms, representation='NewCartoon'): temp_script.write(vmd_script.encode('utf-8')) - command = [vmd_path, '-e', str(temp_script_path), '-args', str(temp_pdb_path), str(output_path)] + command = [vmd_path, '-e', str(temp_script_path), '-args', + str(temp_pdb_path), str(output_path)] try: if PY3K: @@ -175,7 +343,7 @@ def getVmdModel(vmd_path, atoms, representation='NewCartoon'): if returncode != 0: LOGGER.info("VMD exited with status " + str(returncode) + ".") except Exception as e: - LOGGER.warn("An unexpected error occurred: " + str(e)) + _warn("An unexpected error occurred: " + str(e)) finally: if os.path.exists(temp_script_path): os.unlink(temp_script_path) @@ -198,33 +366,41 @@ def getVmdModel(vmd_path, atoms, representation='NewCartoon'): def showChannels(channels, model=None, surface=None): - """Visualizes the channels, and optionally, the molecular model and surface, using Open3D. + """Visualizes the channels, and optionally, the molecular model and + surface, using Open3D. - This function renders a 3D visualization of molecular channels based on their spline representations. - It can also display a molecular model (e.g., the protein structure) and a surface (e.g., cavity surface) - in the same visualization. The function utilizes the Open3D library to create and render the 3D meshes. + This function renders a 3D visualization of molecular channels based on + their spline representations. It can also display a molecular model (e.g., + the protein structure) and a surface (e.g., cavity surface) in the same + visualization. The function utilizes the Open3D library to create and + render the 3D meshes. To install Open3D use: - conda install open3d (for Anaconda users; version open3d-0.19.0 was used during the developement) - or pip install open3d + conda install open3d (for Anaconda users; version open3d-0.19.0 was used + during the development) or pip install open3d - :arg channels: A list of channel objects or a single channel object. Each channel should have a - `get_splines()` method that returns two CubicSpline objects: one for the centerline and one for the radii. + :arg channels: A list of channel objects or a single channel object. Each + channel should have a `getSplines()` method that returns two + CubicSpline objects: one for the centerline and one for the radii. :type channels: list or single channel object - :arg model: An optional Open3D TriangleMesh object representing the molecular model, such as a protein. - If provided, this model will be rendered in the visualization. + :arg model: An optional Open3D TriangleMesh object representing the + molecular model, such as a protein. If provided, this model will be + rendered in the visualization. Model can be generated using getVmdModel() function. :type model: open3d.geometry.TriangleMesh, optional - :arg surface: An optional list containing the surface data. The list should have two elements: + :arg surface: An optional list containing the surface data. The list should + have two elements: - `points`: The coordinates of the vertices on the surface. - - `simp`: The simplices that define the surface (e.g., triangles or tetrahedra). - If provided, the surface will be rendered as a wireframe overlay in the visualization. + - `simp`: The simplices that define the surface (e.g., triangles or + tetrahedra). + If provided, the surface will be rendered as a wireframe overlay in the + visualization. :type surface: list (with two numpy arrays), optional - :raises ImportError: If the Open3D library is not installed, an ImportError is raised, - prompting the user to install Open3D. + :raises ImportError: If the Open3D library is not installed, an ImportError + is raised, prompting the user to install Open3D. :returns: None. This function only renders the visualization. @@ -243,7 +419,11 @@ def create_mesh_from_spline(centerline_spline, radius_spline, n=5): centers = centerline_spline(t) radii = radius_spline(t) - spheres = [o3d.geometry.TriangleMesh.create_sphere(radius=r, resolution=20).translate(c) for r, c in zip(radii, centers)] + spheres = [ + o3d.geometry.TriangleMesh.create_sphere(radius=r, + resolution=20).translate(c) + for r, c in zip(radii, centers) + ] mesh = spheres[0] for sphere in spheres[1:]: mesh += sphere @@ -253,7 +433,7 @@ def create_mesh_from_spline(centerline_spline, radius_spline, n=5): if not isinstance(channels, list): channels = [channels] - channel_meshes = [create_mesh_from_spline(*channel.get_splines()) for channel in channels] + channel_meshes = [create_mesh_from_spline(*channel.getSplines()) for channel in channels] meshes_to_visualize = [o3d.geometry.TriangleMesh.create_coordinate_frame(size=0.1, origin=[0, 0, 0])] if model is not None: @@ -284,7 +464,8 @@ def create_mesh_from_spline(centerline_spline, radius_spline, n=5): triangles.sort(axis=1) triangles_tuple = [tuple(tri) for tri in triangles] - unique_triangles, counts = np.unique(triangles_tuple, return_counts=True, axis=0) + unique_triangles, counts = np.unique(triangles_tuple, + return_counts=True, axis=0) surface_triangles = unique_triangles[counts == 1] @@ -309,27 +490,31 @@ def create_mesh_from_spline(centerline_spline, radius_spline, n=5): def showCavities(surface, show_surface=False): """Visualizes the cavities within a molecular surface using Open3D. - This function displays a 3D visualization of cavities detected in a molecular structure. - It uses the Open3D library to render the cavities as a triangle mesh. Optionally, it can also - display the molecular surface as a wireframe overlay. + This function displays a 3D visualization of cavities detected in a + molecular structure. + It uses the Open3D library to render the cavities as a triangle mesh. + Optionally, it can also display the molecular surface as a wireframe + overlay. To install Open3D use: - conda install open3d (for Anaconda users; version open3d-0.19.0 was used during the developement) - or pip install open3d + conda install open3d (for Anaconda users; version open3d-0.19.0 was used + during the development) or pip install open3d :arg surface: A list containing three elements: - - `points`: The coordinates of the vertices (atoms) in the molecular structure. + - `points`: The coordinates of the vertices (atoms) in the molecular + structure. - `surf_simp`: The simplices that define the molecular surface. - `simp_cavities`: The simplices corresponding to the detected cavities. :type surface: list (with three numpy arrays) - :arg show_surface: A boolean flag indicating whether to display the molecular surface - as a wireframe overlay in the visualization. If True, the surface will be displayed - in addition to the cavities. Default is False. + :arg show_surface: A boolean flag indicating whether to display the + molecular surface + as a wireframe overlay in the visualization. If True, the surface will + be displayed in addition to the cavities. Default is False. :type show_surface: bool - :raises ImportError: If the Open3D library is not installed, an ImportError is raised, - prompting the user to install Open3D. + :raises ImportError: If the Open3D library is not installed, an ImportError + is raised, prompting the user to install Open3D. :returns: None @@ -353,7 +538,8 @@ def showCavities(surface, show_surface=False): sorted([tetra[0], tetra[2], tetra[3]]), sorted([tetra[1], tetra[2], tetra[3]])]) - surface_triangles = np.unique(np.array(triangles), axis=0, return_counts=True)[0] + surface_triangles = np.unique(np.array(triangles), axis=0, + return_counts=True)[0] mesh = o3d.geometry.TriangleMesh() mesh.vertices = o3d.utility.Vector3dVector(points) @@ -378,7 +564,8 @@ def showCavities(surface, show_surface=False): triangles.sort(axis=1) triangles_tuple = [tuple(tri) for tri in triangles] - unique_triangles, counts = np.unique(triangles_tuple, return_counts=True, axis=0) + unique_triangles, counts = np.unique(triangles_tuple, + return_counts=True, axis=0) surface_triangles = unique_triangles[counts == 1] @@ -421,7 +608,7 @@ def showSurfaceCavities(surface, cavities=None, model=None, show_surface=False, coordinates of the supplied pseudoatoms instead of from `surface[2]` or `cavities`. - :param surface: Surface data returned by :func:`calcSurfaceCavities`. + :arg surface: Surface data returned by :func:`calcSurfaceCavities`. Required for `mode='tetra'`, for `mode='smooth'` when `cavity_atoms` is not provided, and for displaying the molecular surface when `show_surface=True`. The expected list contains: @@ -434,42 +621,43 @@ def showSurfaceCavities(surface, cavities=None, model=None, show_surface=False, :type surface: list or None - :param cavities: List of :class:`Cavity` objects returned by + :arg cavities: List of :class:`Cavity` objects returned by :func:`calcSurfaceCavities`. Required when `mode='smooth'` and `cavity_atoms` is not provided, because the function uses `cavity.tetrahedra` to select the corresponding Voronoi vertices. :type cavities: list or None - :param model: Optional Open3D `TriangleMesh` representing the protein or + :arg model: Optional Open3D `TriangleMesh` representing the protein or another molecular model. The model can be generated with :func:`getVmdModel`. :type model: open3d.geometry.TriangleMesh, or None - :param show_surface: If `True`, display the molecular surface wireframe + :arg show_surface: If `True`, display the molecular surface wireframe derived from `surface[1]` in addition to the cavity representation. This requires `surface` to be provided. Default is `False`. :type show_surface: bool - :param mode: Visualization mode used when `cavity_atoms` is not provided. + :arg mode: Visualization mode used when `cavity_atoms` is not provided. Accepted values are `'tetra'` and `'smooth'`. Default is `'tetra'`. :type mode: str - :param alpha: Alpha value used for alpha-shape surface reconstruction in + :arg alpha: Alpha value used for alpha-shape surface reconstruction in `mode='smooth'` and when visualizing `cavity_atoms`. Smaller values produce tighter surfaces, while larger values may connect more distant points and generate broader surfaces. Default is 4.0. :type alpha: float - :param smoothing: Number of Taubin smoothing iterations applied to the + :arg smoothing: Number of Taubin smoothing iterations applied to the reconstructed cavity mesh. If `0` or `None`, no smoothing is applied. Default is 0. :type smoothing: int or None - :param cavity_atoms: Optional pseudoatom representation of surface + :arg cavity_atoms: Optional pseudoatom representation of surface cavities. This can be either a path to a PDB/PQR file or a parsed ProDy `AtomGroup`, or an Open3D `TriangleMesh` generated, for example, with :func:`getVmdModel`. - :type cavity_atoms: str, :class:`.AtomGroup`, open3d.geometry.TriangleMesh, or None + :type cavity_atoms: str, :class:`.AtomGroup`, open3d.geometry.TriangleMesh, + or None Examples: p = parsePDB('1tqn') @@ -522,8 +710,8 @@ def showSurfaceCavities(surface, cavities=None, model=None, show_surface=False, cavity_atoms = parsePDB(cavity_atoms) if not hasattr(cavity_atoms, 'getCoords'): - raise TypeError("cavity_atoms must be a PDB/PQR filename, a ProDy AtomGroup, " - "or an Open3D TriangleMesh.") + raise TypeError("cavity_atoms must be a PDB/PQR filename, a " + "ProDy AtomGroup,or an Open3D TriangleMesh.") resnums = np.unique(cavity_atoms.getResnums()) @@ -564,7 +752,8 @@ def showSurfaceCavities(surface, cavities=None, model=None, show_surface=False, sorted([tetra[0], tetra[2], tetra[3]]), sorted([tetra[1], tetra[2], tetra[3]])]) - surface_triangles = np.unique(np.array(triangles), axis=0, return_counts=True)[0] + surface_triangles = np.unique(np.array(triangles), axis=0, + return_counts=True)[0] cavity_mesh = o3d.geometry.TriangleMesh() cavity_mesh.vertices = o3d.utility.Vector3dVector(points) cavity_mesh.triangles = o3d.utility.Vector3iVector(surface_triangles) @@ -610,7 +799,8 @@ def showSurfaceCavities(surface, cavities=None, model=None, show_surface=False, triangles.sort(axis=1) triangles_tuple = [tuple(tri) for tri in triangles] - unique_triangles, counts = np.unique(triangles_tuple, return_counts=True, axis=0) + unique_triangles, counts = np.unique(triangles_tuple, + return_counts=True, axis=0) surface_triangles = unique_triangles[counts == 1] lines = [] @@ -626,98 +816,330 @@ def showSurfaceCavities(surface, cavities=None, model=None, show_surface=False, o3d.visualization.draw_geometries(meshes_to_visualize) +def calcChannels(atoms, output_path=None, separate=False, start_point=None, + restrict_channels_to_start_point=True, start_point_search=3.0, + r1=3, r2=0.9, min_depth=5, + min_volume=None, max_volume=None, max_depth=None, bottleneck=0.0, + sparsity=1, min_tetrahedra=None, max_tetrahedra=None, cavities_only=False, + diagram="homogenized", max_deviation=0.1, truncate_at_surface=True, + similarity=0.8, route_tolerance=1.0, min_enclosure=0.70, max_peel_depth=None, + weighted_cache=True, weighted_mouth_depth=2.5, edge_cost=None): + """Computes and identifies channels within a molecular structure using + Voronoi and Delaunay tessellations. + + This function analyzes the provided atomic structure to detect channels, + which are voids or pathways within the molecular structure. It employs + Voronoi and Delaunay tessellations to identify these regions, then filters + and refines the detected channels based on various parameters such as the + minimum depth and bottleneck size. The results can be saved to a PQR file + (PDB is optional) if an output path is provided. The `separate` parameter + controls whether each detected channel is saved to a separate file or if all + channels are saved in a single file. -def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3, r2=1.25, min_depth=10, - min_volume=None, max_volume=None, max_depth=None, bottleneck=1, sparsity=15, - min_tetrahedra=None, max_tetrahedra=None, cavities_only=False): - """Computes and identifies channels within a molecular structure using Voronoi and Delaunay tessellations. + The implementation is inspired by the methods described in the following + publications: + "MOLE 2.0: advanced approach for analysis of biomacromolecular channels" by + D. Sehnal, et al., published in J Chemoinform, 5 (39) 2013. + + "CAVER: Algorithms for Analyzing Dynamics of Tunnels in Macromolecules". by + A. Pavelka, et al., published in IEEE ACM T COMPUT BI, (13) 2016. - This function analyzes the provided atomic structure to detect channels, which are voids or pathways - within the molecular structure. It employs Voronoi and Delaunay tessellations to identify these regions, - then filters and refines the detected channels based on various parameters such as the minimum depth - and bottleneck size. The results can be saved to a PQR file (PDB is optional) if an output path is provided. - The `separate` parameter controls whether each detected channel is saved to a separate file or if all - channels are saved in a single file. + "Software Tools for Identification, Visualization and Analysis of Protein + Tunnels and Channels". by J. Brezovsky, et al., Biotechnol Adv (31) 2013. - The implementation is inspired by the methods described in the publication: - "MOLE 2.0: advanced approach for analysis of biomacromolecular channels" by D. Sehnal, et al., published in - J Chemoinform, 5 (39) 2013. - :param atoms: An object representing the molecular structure, typically containing atomic coordinates - and element types. + :arg atoms: An object representing the molecular structure, typically + containing atomic coordinates and element types. :type atoms: `Atoms` object - :param output_path: Optional path to save the resulting channels and associated data in PQR (or PDB) format. - If None, results are not saved. Default is None. + :arg output_path: Optional path to save the resulting channels and + associated data in PQR (or PDB) format. If None, results are not saved. + Default is None. :type output_path: str or None - :param separate: If True, each detected channel is saved to a separate PDB file. If False, all channels - are saved in a single PDB file. Default is False. + :arg separate: If True, each detected channel is saved to a separate PDB + file. If False, all channels are saved in a single PDB file. Default is + False. :type separate: bool - :param start_point: Optional starting point for channel search. This can be either a 3D coordinate point or - an atomic selection/AtomGroup. If the 3D coordinate point will be provided, the algorithm will use the - tetrahedron whose Voronoi vertex is closest to this point as the starting tetrahedron (overriding - the default automatic seed selection based on the deepest tetrahedron). Coordinates must be given in Å. - If an atomic selection is provided, its geometric center is used as the starting point. - :type start_point: list, tuple, or ndarray (length 3), :class:`.Atomic`, or None - - :param r1: The first radius threshold used during the deletion of simplices, which is used to define - the outer surface of the channels. Default is 3. + :arg start_point: Optional starting point for channel search. This can be + either a 3D coordinate point or an atomic selection/AtomGroup. If the + 3D coordinate point will be provided, the algorithm seeds the channel + search near this point (overriding the default automatic seed selection + based on the deepest tetrahedron); see ``start_point_search`` for how + the seed tetrahedron itself is picked. Coordinates must be given in Å. + If an atomic selection is provided, its geometric center is used as the + starting point. + :type start_point: list, tuple, or ndarray (length 3), :class:`.Atomic`, or None + + :arg restrict_channels_to_start_point: Only used when ``start_point`` is + provided. If True (default), the channel search is restricted to the + single cavity whose closest tetrahedron is globally nearest to + ``start_point``, so channels are computed only for the region around + that point instead of one channel bundle per detected cavity. If False, + ``start_point`` merely overrides the seed (starting) tetrahedron of + every cavity and channels are still computed for all cavities. + :type restrict_channels_to_start_point: bool + + :arg start_point_search: Only used when ``start_point`` is provided. Radius, + in Angstrom, of the neighbourhood of ``start_point`` searched for the seed + tetrahedron. The tetrahedron nearest ``start_point`` is often a tight one, + and since every channel of the cavity starts there, its inscribed radius + caps all of their bottlenecks and appears as one shared bottleneck at the + joint beginning of the bundle. Seeded instead is the widest tetrahedron within + ``start_point_search`` of ``start_point`` that belongs to the same cavity, is + no shallower than the nearest one (so the seed cannot drift out towards the + mouth) and is reachable from it through that neighbourhood (so it stays in the + void the point sits in rather than crossing a wall). Default is 3.0; use 0 to + seed the nearest tetrahedron as-is. + :type start_point_search: float + + :arg r1: The first radius threshold used during the deletion of simplices, + which is used to define the outer surface of the channels. Default is 3 :type r1: float - :param r2: The second radius threshold used to define the inner surface of the channels. Default is 1.25. + :arg r2: The second radius threshold used to define the inner surface of + the channels. Default is 0.9. + + Below about 1.2 Angstrom the probe is smaller than a water molecule, and + then the structure must carry explicit hydrogens. Without them, every + carbon keeps its full vdW radius while the space its hydrogens occupied is + left empty, and a sub-water probe is small enough to thread those + interstices: the interior percolates into a sponge and the channel count + can rise several-fold. At 1.2 Angstrom and above, protonated and + unprotonated structures give the same channels, and an X-ray file may + be used as it comes. A warning is issued for the unsafe combination. :type r2: float - :param min_depth: The minimum depth a cavity must have to be considered as a channel. Default is 10. - :type min_depth: int + :arg min_depth: The minimum depth, in Angstrom, a cavity must reach to be + considered as a channel. Depth is the geodesic distance from the cavity's + surface opening to its farthest point along the Voronoi network, so it is a + physical length independent of the tessellation density. Default is 5. + :type min_depth: float - :param max_depth: Maximum cavity depth. Cavities deeper than this value are trimmed to the specified depth. - Default is None. - :type max_depth: int + :arg max_depth: Maximum cavity depth, in Angstrom. Portions of a cavity deeper + than this value are trimmed away. Default is None (no trimming). + :type max_depth: float - :param bottleneck: The minimum allowed bottleneck size (narrowest point) for the channels. Default is 1. + :arg bottleneck: Acts as secondary filter following channel identification. + The minimum allowed bottleneck size (narrowest point) for the channels. + It it critical when diagram=simple, as it partially corrects for wrong + diagram topology. Default is 0.0, no filtering applied. :type bottleneck: float - :param min_volume: Minimum volume required for a channel/cavity to be retained. Default is None. + :arg min_volume: Minimum volume required for a channel/cavity to be + retained. Default is None. :type min_volume: float - :param max_volume: Maximum volume allowed for a channel/cavity to be retained. Default is None. + :arg max_volume: Maximum volume allowed for a channel/cavity to be + retained. Default is None. :type max_volume: float - :param sparsity: The sparsity parameter controls the sampling density when analyzing the molecular surface. - A higher value results in fewer sampling points. Default is 15. - :type sparsity: int + :arg sparsity: Size of a channel surface opening (mouth), in Angstrom: how far + apart two exits must lie to count as separate openings. It is one quantity, + reached by whichever branch of the search is running, and the two branches + are mutually exclusive. With ``truncate_at_surface`` True (the default) it + is a floor on the radius of a reported opening, so two channels leaving + closer than ``sparsity`` are treated as sharing that opening and are merged + if they also share a corridor (see ``similarity``); being applied *after* + the search, it can only merge channels there, never hide one, and is a + reporting preference rather than part of the geometry. With + ``truncate_at_surface`` False it is instead the spacing at which exit + tetrahedra are sampled as channel termini, and it does then decide which + channels are found at all. Either way a higher value reports fewer channels. + It has no effect on the cavities, which are found from the exit tetrahedra + before any thinning. Default is 1. + :type sparsity: float + + :arg diagram: + "homogenized" (default) - every atom is substituted by a set of + homogeneous balls whose common radius equals the smallest van der Waals + radius present in the structure, before building the Voronoi and + Delaunay tessellations. This yields an accurate estimate of the + additively weighted Voronoi diagram from an ordinary one, as done in + MolAxis and CAVER 3 + "simple" - the original atoms are used with their individual van der + Waals radii directly. This is very inaccurate and should be avoided in + almost all cases. + "weighted" - the true additively-weighted (Apollonius) Voronoi diagram is + built directly from the atoms with their individual van der Waals radii, + using the third-party ``vorpy`` package (with a compiled kernel when + ``numba`` is available). This is the exact diagram the "homogenized" mode + approximates, at a higher cost (~ 100x slower, for 4000 atoms). To use this + approach install ``vorpy`` library using ``pip install vorpy3``. + :type diagram: str + + :arg max_deviation: Maximum tolerated deviation, in Angstrom, between the + union surface of the substitute balls and the original van der Waals + surface when ``diagram = homogenized`` . It controls the trade-off + between surface accuracy and the number of balls generated: an atom + whose radius exceeds the smallest radius (``rho``) by more than + ``max_deviation`` is filled with several balls, otherwise it is kept as + a single ``rho`` ball. Default is 0.1. Guideline values: + + * ``0.1`` fine accurate surface with minimal errors, but on 13-15x + more balls than original + * ``0.15`` in heavy-atom-only structures it starts filling carbon, which + is otherwise left as a single ball with a uniform ~0.18 A inset). + * ``0.2`` speed optimized ; e.g. carbon fills to ~15 balls when + hydrogens are present (``rho``=1.2), resulting roughly to ~8 times + more balls. Without hydrogens, carbons are single balls. + + Only used when ``diagram = homogenized``. + :type max_deviation: float + + :arg truncate_at_surface: If True (default), surface (exit) tetrahedra are + made *absorbing*: a channel may end at one, but no channel may pass + through one. A mouth is a surface tetrahedron a probe of the traversal + radius ``r2`` can leave through. This forbids the cheapest path from + surfacing at one mouth, running along the outside and re-entering at + another - a surface hop, not a tunnel - which the width-rewarding cost + would otherwise prefer, since surface grooves are the widest space + available. Enforcing it during the search (rather than cutting the + winning path afterwards) is what keeps genuine narrow interior corridors + in the output: cut afterwards, such a corridor loses the cheapest-path + race to the groove leading to the same mouth and is never enumerated at + all. If False, paths run freely to their end tetrahedra, surface hops + included. + :type truncate_at_surface: bool + + :arg similarity: Only used when ``truncate_at_surface`` is True. Fraction + (0-1) of the **longer** of two channels, measured in Angstrom along its + centerline, that must run within ``route_tolerance`` of the other one for + the two to count as the same corridor. Two channels are merged (cheapest + kept) only when they take the same corridor **and** leave through the + same opening (see ``sparsity``); a corridor that forks near the surface + and exits twice is one tunnel, but two different corridors to one opening, + or one corridor reaching two openings, are two tunnels. The comparison is + geometric rather than a shared prefix of tetrahedra, so it is unaffected + by *where* two routes diverge (variants that split and rejoin still count + as one), by containment (a long route is not deleted as the "duplicate" of + a short one it happens to start with), and by ``max_deviation`` (a + tetrahedron count is not mesh-invariant; Angstrom are). ``1.0`` merges + only routes that coincide along their whole length; ``0.0`` merges every + channel that shares an opening. Default is 0.8. + :type similarity: float + + :arg route_tolerance: Only used when ``truncate_at_surface`` is True. How far + apart, in Angstrom, two centerlines may drift and still count as the same + corridor when computing ``similarity``. Larger values merge more + aggressively (nearby parallel routes read as one tunnel); smaller values + report finer route variants separately. Default is 1.0. + :type route_tolerance: float + + :arg min_enclosure: Fraction of directions that must be blocked by protein for + a tetrahedron to count as interior, in ``[0, 1]``. Default is 0.70. + + Once the r1 surface is built it is eroded inward with the r2 probe, to + strip the shell of true exterior that an r1 probe bridges over rather than + enters (the "moat"); that shell would otherwise join the cavity and offer + wide, low-cost routes along the outside of the protein. Erosion continues + while the tetrahedra at the front are *open*, meaning that fewer than + ``min_enclosure`` of the directions leaving them run into protein within + :data:`ENCLOSURE_RANGE` Angstrom, and halts at the first buried layer. + + Bounding the erosion by size instead does not work. A count of tetrahedron + layers is not mesh-invariant, as a layer is one tetrahedron thick and + tetrahedra shrink as ``max_deviation`` is lowered. A depth in Angstrom is + not ``r1``-invariant, as the moat is as deep as ``r1 - r2`` in a concavity + but vanishes on a flat face, so a depth that clears it where it is thick + also marches down the channel mouths and erodes the channels themselves. + Testing burial locally instead leaves ``r1`` to decide only where the + erosion starts, not where it stops, so results are independent of it, and + ``r1`` is left doing the one job it should: capping the mouths. + + It is bounded from both sides, and the window is narrow. + + Too low and the moat is not fully stripped. Too high and the erosion + never meets a layer buried enough to stop it, so it percolates down the + channels and eats the cavity. The default sits at the floor, which is + the safe end: over-peeling deletes real channels, whereas + under-peeling shows up as surface-riding routes that can be recognised. + :type min_enclosure: float + + :arg max_peel_depth: Optional hard cap, **in Angstrom**, on how far the peel + above may advance from the r1 surface. ``None`` (default) is uncapped, and + the enclosure test alone decides where erosion stops. Set it only as a + backstop on a structure where the peel misbehaves; it is deliberately not + tied to ``r1``, since a cap that scales with ``r1`` reintroduces exactly + the ``r1`` dependence that ``min_enclosure`` exists to remove. + :type max_peel_depth: float or None + + :arg weighted_cache: Cache the raw additively-weighted Voronoi diagram to disk so + that re-running ``diagram="weighted"`` on the same structure skips the + expensive vorpy tessellation (which dominates the ~10 min run time). The + diagram only depends on the atoms and ``r1``, so re-runs that change only + ``r2``, ``bottleneck``, ``sparsity``, ``start_point`` etc. reuse it. ``True`` + (default) caches next to ``output_path`` (or, absent one, under the structure + title in the current directory); pass a path to place it explicitly, or + ``False`` to disable. The cache is keyed by content, so editing the structure + or ``r1`` transparently forces a recompute. Only used for ``diagram="weighted"``. + :type weighted_cache: bool or str + + :arg weighted_mouth_depth: Only used for ``diagram="weighted"``. The additively- + weighted (Apollonius) tessellation is not a clean simplicial complex, so it + leaves false interior boundary faces that the pipeline would misread as + surface openings, truncating channels to stubs. To repair this a *homogenized* + diagram of the same atoms is built as an interior/exterior oracle, and only + exit tetrahedra whose Voronoi vertex lies within ``weighted_mouth_depth`` + Angstrom (geodesic distance below the molecular surface) are treated as mouths. + Default 2.5 (the value at which the recovered channels match the + ``"homogenized"`` result); ``None`` disables the relabeling. + :type weighted_mouth_depth: float or None + + :arg edge_cost: How each Voronoi edge is priced in the Dijkstra tunnel search. + ``"integral"`` prices each edge by the integral of its clearance profile + along the edge, which is mesh-invariant. ``"bottleneck"`` uses the legacy + ``length / (gate**2 + eps)``, charging the whole edge at its single + narrowest point, whose value (and the routing it produces) drifts as + ``max_deviation`` coarsens. ``None`` (default) selects ``"integral"`` for + ``diagram="homogenized"``/``"simple"`` (straight edges, where the integral + is exact) and ``"bottleneck"`` for ``diagram="weighted"``. ``"integral"`` + is rejected for ``diagram="weighted"`` (Apollonius edges are arcs the + straight-chord integral cannot price). The reported bottleneck radius is + unaffected by this choice. + :type edge_cost: str or None :returns: A tuple containing two elements: - - `channels`: A list of detected channels, where each channel is an object containing information - about its path and geometry. - - `surface`: A list containing additional information for further visualization, including - the atomic coordinates, simplices defining the surface, and merged cavities. + - `channels`: A list of detected channels, where each channel is an + object containing information about its path and geometry. + - `surface`: A list containing additional information for further + visualization, including the atomic coordinates, simplices defining + the surface, and merged cavities. :rtype: tuple (list, list) This function performs the following steps: - 1. **Selection and Filtering:** Selects non-hetero atoms from the protein, calculates van der Waals radii, - and performs 3D Delaunay triangulation and Voronoi tessellation on the coordinates. - 2. **State Management:** Creates and updates different stages of channel detection of the protein structure - to filter out simplices based on the given radii. - 3. **Surface Layer Calculation:** Determines the surface and second-layer simplices from the filtered results. - 4. **Cavity and Channel Detection:** Finds and filters cavities based on their depth and calculates channels - using Dijkstra's algorithm. - 5. **Visualization and Saving:** Generates meshes for the detected channels, filters them by bottleneck size, - and either saves the results to a PDB file or visualizes them based on the specified parameters. + 1. **Selection and Filtering:** Selects non-hetero atoms from the protein + and calculates van der Waals radii. When ``homogenize`` is True, each + atom is replaced by homogeneous balls of the smallest radius present + so that an ordinary tessellation approximates the additively weighted + Voronoi diagram. It then performs 3D Delaunay triangulation and Voronoi + tessellation on the resulting coordinates. + 2. **State Management:** Creates and updates different stages of channel + detection of the protein structure to filter out simplices based on the + given radii. + 3. **Surface Layer Calculation:** Determines the surface and second-layer + simplices from the filtered results. + 4. **Cavity and Channel Detection:** Finds and filters cavities based on + their depth and calculates channels using Dijkstra's algorithm. + 5. **Visualization and Saving:** Generates meshes for detected channels, + filters them by bottleneck size, and either saves the results to a PDB + file or visualizes them based on the specified parameters. Example usage: channels, surface = calcChannels(atoms, output_path="channels", separate=True) - channels, surface = calcChannels(atoms, output_path="all_channels.pdb", start_point=[-22.312, -20.065, -11.144]) + channels, surface = calcChannels(atoms, output_path="all_channels.pdb", + start_point=[-22.312, -20.065, -11.144]) start_sel = protein.select('resid 212 309 483') - channels, surface = calcChannels(atoms, output_path="all_channels.pdb", start_point=start_sel) + channels, surface = calcChannels(atoms, output_path="all_channels.pdb", + start_point=start_sel) To save the results as PDB file: - channels, surface = calcChannels(atoms, output_path="channels.pdb", separate=False, r1=3, r2=1.25, min_depth=10, - bottleneck=1, sparsity=15) """ + channels, surface = calcChannels(atoms, output_path="channels.pdb", + separate=False, r1=3, r2=0.9, min_depth=5, + bottleneck=1, sparsity=3) """ required = ['heapq', 'collections', 'scipy', 'pathlib', 'warnings'] missing = [] @@ -735,8 +1157,6 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 errorMsg = ', '.join(errorMsg.split(', ')[:-1]) + ' and ' + errorMsg.split(', ')[-1] raise ImportError(errorMsg) - from scipy.spatial import Voronoi, Delaunay - if PY3K: from pathlib import Path else: @@ -763,67 +1183,224 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 LOGGER.info("Using user-provided start_point for channel seed: [{:.3f}, {:.3f}, {:.3f}] Å" .format(start_point[0], start_point[1], start_point[2])) - calculator = ChannelCalculator(atoms, r1, r2, min_depth, bottleneck, sparsity) + LOGGER.timeit('_prody_calcChannels') + + # Edge-cost mode for the Dijkstra routing (buildSparseGraph). Default: the + # mesh-invariant profile integral for the straight-edged homogenized/simple + # diagrams, the legacy l/(d^2+b) for weighted (Apollonius edges are arcs, + # where the straight-chord integral is only approximate; weighted is + # experimental). An explicit value overrides the per-diagram default. + if edge_cost is None: + edge_cost = 'bottleneck' if diagram == 'weighted' else 'integral' + elif edge_cost not in ('integral', 'bottleneck'): + raise ValueError("edge_cost must be 'integral', 'bottleneck' or None, " + "got {0!r}".format(edge_cost)) + elif edge_cost == 'integral' and diagram == 'weighted': + raise ValueError("edge_cost='integral' is only valid for the straight-edge " + "diagrams ('homogenized'/'simple'); the weighted " + "(Apollonius) diagram has arc edges the straight-chord " + "integral cannot price. Use edge_cost='bottleneck' (the " + "default for diagram='weighted') or None.") - atoms = atoms.select('not hetero and noh') # Excluding hydrogens + _reportAtomsInputComposition(atoms) + atoms = atoms.select('not water') # water is excluded from the selection + calculator = ChannelCalculator(atoms, r2=r2, sparsity=sparsity, + route_tolerance=route_tolerance, + edge_cost=edge_cost) + + elements = np.char.upper(np.asarray(atoms.getElements(), dtype=str)) + has_hydrogens = bool(np.any(elements == 'H')) + + # An experimental structure generally carries no hydrogens -- X-ray and cryo-EM + # alike, since neither resolves them except at the very highest resolutions -- and + # its carbons keep their full vdW radius, so the ~0.6 A the missing H occupied is + # left as void, around every heavy atom at once, including buried contacts that + # never come apart. That is usually harmless, and is often defended as standing in + # for thermal motion: a probe of water size cannot enter those interstices anyway, + # and protonated and unprotonated runs agree from about 1.2 A upwards. Below that + # the probe is small enough to thread them and the interior percolates into a + # sponge rather than merely widening. Those routes are fictitious, not the real + # ones made wider. So a sub-water probe needs real hydrogens; above it, take the + # file as it comes. + if not has_hydrogens and r2 < 1.2: + _warn("structure has no hydrogens and r2={0:.2f} is below 1.2 A: the space " + "left by the missing H is then wide enough for the probe to pass, and " + "channels will be found through interstices that do not exist in the " + "real protein (their number can rise several-fold). Either add " + "hydrogens, or raise r2 to 1.2 A or more, where protonated and " + "unprotonated structures give the same channels.".format(r2)) + + if diagram == "simple": + # 'simple' builds an *unweighted* Delaunay of the atom centres, i.e. it + # ignores the differences between atomic radii. That approximation is worst + # when the radius spread is largest -- which is exactly when hydrogens (small + # vdW) are present -- so warn there and steer the user to a radius-aware mode. + # With H absent the heavy-atom radii are much closer, so 'simple' is more + # defensible and matches the heavy-atom-only input most tools accept (at the + # cost of over-large empty space where the missing H would sit). + if has_hydrogens: + _warn("diagram='simple' with hydrogens present: the unweighted " + "Voronoi diagram ignores radius differences, which are largest when H " + "are present, so its topology and clearances are significantly less " + "accurate. Consider diagram='homogenized' (or 'weighted'), which " + "account for per-atom radii.") + coords = atoms.getCoords() - - vdw_radii = calculator.get_vdw_radii(atoms.getElements()) - - dela = Delaunay(coords) - voro = Voronoi(coords) - - s_prt = State(dela.simplices, dela.neighbors, voro.vertices) + vdw_radii = calculator.getVdwRadii(atoms.getElements()) + # Burial is a property of the protein, not of the tessellation, so the enclosure + # test that strips the moat runs against the real atoms rather than the balls the + # diagram happens to be built on. Homogenization would otherwise make it a + # function of max_deviation, which min_enclosure must not be. + atom_coords = coords + # For diagram="weighted" only: a homogenized-surface depth oracle used to relabel + # the additively-weighted diagram's leaky surface mouths (see getSurfaceCavities). + mouth_oracle = None + + if diagram == "homogenized": + LOGGER.timeit('_prody_channels_homogenize') + coords, vdw_radii = calculator.homogenizeAtoms(coords, vdw_radii, max_deviation) + LOGGER.report("Substituted {0} atoms with {1} homogeneous balls of radius {2:.2f} A in %.2fs.".format( + atoms.numAtoms(), len(coords), float(vdw_radii[0])), '_prody_channels_homogenize') + + LOGGER.timeit('_prody_channels_tessellation') + if diagram == "weighted": + # True additively-weighted (Apollonius) Voronoi network via the third-party + # vorpy package: van der Waals radii are baked into the diagram exactly, + # instead of being approximated by homogenising atoms into uniform balls. + # buildAwTessellation returns the same (simplices, neighbors, vertices) triple + # a scipy Delaunay would, so the downstream erosion/cavity pipeline is + # untouched. Because every AW vertex is equidistant (additively) to its 4 + # tangent atoms, the sum-based clearance test in deleteSimplices3d reduces + # exactly to the per-atom clearance, so the returned clearances are not needed. + if not checkAndImport('vorpy'): + raise ImportError('diagram="weighted" requires the vorpy package for the ' + 'additively-weighted (Apollonius) Voronoi diagram. Install vorpy, or ' + 'use diagram="homogenized"/"simple".') + # The compiled calc_vert kernel needs numba; without it the weighted path + # still works but is ~5x slower, so fall back with a warning rather than fail. + accelerate = checkAndImport('numba') + if not accelerate: + _warn('numba is not installed; the additively-weighted tessellation ' + 'will run without the compiled kernel and may be very slow.') + from ._vorpy_aw import buildAwTessellation, resolveCachePath + try: + title = atoms.getTitle() + except Exception: + title = None + cache_path = resolveCachePath(weighted_cache, output_path, title) + simplices, neighbors, verts, _ = buildAwTessellation( + coords, vdw_radii, max_vert=max(2.0 * r1, 8), accelerate=accelerate, + cache=cache_path) + LOGGER.report('Additively-weighted (Apollonius) tessellation of {0} atoms ' + 'constructed in %.2fs.'.format(len(coords)), + '_prody_channels_tessellation') + # The AW->simplicial mapping leaves false interior boundary faces that the + # pipeline would misread as surface mouths (collapsing channels to stubs). + # Build a homogenized diagram of the same atoms as an interior/exterior depth + # oracle; getSurfaceCavities then keeps only exit tetrahedra within + # weighted_mouth_depth Angstrom (geodesic) of the true molecular surface. + if weighted_mouth_depth is not None: + LOGGER.timeit('_prody_channels_mouth_oracle') + mouth_oracle = calculator.buildSurfaceDepthOracle( + coords, vdw_radii, r1, max_deviation, weighted_mouth_depth) + LOGGER.report('Homogenized surface oracle (weighted mouth relabeling) ' + 'built in %.2fs.', '_prody_channels_mouth_oracle') + else: + from scipy.spatial import Delaunay + # We deliberately do NOT joggle/jitter the input (no QJ), unlike CAVER, + # which perturbs by ~0.001 A to dodge the cospherical "T5" degeneracy it + # reports in nearly every structure. scipy's default Qhull options + # (Qbb Qc Qz) merge cospherical facets instead of joggling, so the + # circumcenters stay finite even at degeneracies (verified: 0 NaN/inf + # across millions of tetrahedra, even under heavy homogenized refinement). + # Not joggling keeps the pipeline exactly reproducible (homogenizeAtoms is + # a fixed Fibonacci lattice, and nothing here uses an RNG). The only + # residual is coincident circumcenters at true degeneracies, handled where + # it matters by the twin-tetrahedron guard in _edgeBottleneck. + dela = Delaunay(coords) + # circumcenters straight from the Delaunay paraboloid lifting, so we + # skip the redundant second Qhull pass (scipy Voronoi). Numerically identical + # to voro.vertices for points in general position. + simplices = dela.simplices + neighbors = dela.neighbors + verts = calculator.calcCircumcenters(dela) + LOGGER.report('Delaunay tessellation of {0} points constructed in %.2fs.'.format( + len(coords)), '_prody_channels_tessellation') + + LOGGER.timeit('_prody_channels_surface') + s_prt = State(simplices, neighbors, verts) if PY3K: - s_tmp = State(*s_prt.get_state()) + s_tmp = State(*s_prt.getState()) s_prv = State(None, None, None) else: - s_tmp = apply(State, s_prt.get_state()) + s_tmp = apply(State, s_prt.getState()) s_prv = State(None, None, None) while True: - s_prv.set_state(*s_tmp.get_state()) + s_prv.setState(*s_tmp.getState()) if PY3K: - #s_tmp.set_state(*calculator.delete_simplices3d(coords, *(s_tmp.get_state() + [vdw_radii, r1, True]))) - s_tmp.set_state(*calculator.delete_simplices3d(coords, *(s_tmp.get_state() + tuple([vdw_radii, r1, True])))) + s_tmp.setState(*calculator.deleteSimplices3d(coords, *(s_tmp.getState() + tuple([vdw_radii, r1, True])))) else: - tmp_state = calculator.delete_simplices3d(coords, *(s_tmp.get_state() + [vdw_radii, r1, True])) - s_tmp.set_state(*tmp_state) + tmp_state = calculator.deleteSimplices3d(coords, *(s_tmp.getState() + [vdw_radii, r1, True])) + s_tmp.setState(*tmp_state) if s_tmp == s_prv: break - s_srf = State(*s_tmp.get_state()) - #s_inr = State(*calculator.delete_simplices3d(coords, *(s_srf.get_state() + [vdw_radii, r2, False]))) - s_inr = State(*calculator.delete_simplices3d(coords, *(s_srf.get_state() + tuple([vdw_radii, r2, False])))) - - l_first_layer_simp, l_second_layer_simp = calculator.surface_layer(s_srf.simp, s_inr.simp, s_srf.neigh) - s_clr = State(*calculator.delete_section(l_first_layer_simp, *s_inr.get_state())) - - c_cavities = calculator.find_groups(s_clr.neigh) - c_surface_cavities = calculator.get_surface_cavities(c_cavities, s_clr.simp, l_second_layer_simp, s_clr, coords, vdw_radii, sparsity) - - calculator.find_deepest_tetrahedra(c_surface_cavities, s_clr.neigh) + s_srf = State(*s_tmp.getState()) + + # Moat removal: erode the r1 surface inward with the r2 probe, stripping the shell + # of true exterior that a large r1 probe bridges over instead of entering (it would + # otherwise join the cavity and offer wide, low-cost routes along the outside). + # Erosion stops where the tetrahedra stop being open to the solvent, which is a + # local criterion, so neither the mesh nor r1 sets how deep the peel goes. + s_srf = State(*calculator.peelSurfaceByEnclosure( + coords, *(s_srf.getState() + tuple([vdw_radii, r2, atom_coords, + min_enclosure, max_peel_depth])))) + + s_inr = State(*calculator.deleteSimplices3d(coords, *(s_srf.getState() + tuple([vdw_radii, r2, False])))) + + l_first_layer_simp, l_second_layer_simp = calculator.surfaceLayer(s_srf.simp, s_inr.simp, s_srf.neigh) + s_clr = State(*calculator.deleteSection(l_first_layer_simp, *s_inr.getState())) + LOGGER.report('Surface and inner simplices filtered in %.2fs.', '_prody_channels_surface') + + LOGGER.timeit('_prody_channels_cavities') + c_cavities = calculator.findGroups(s_clr.neigh) + c_surface_cavities = calculator.getSurfaceCavities(c_cavities, s_clr.simp, + l_second_layer_simp, + s_clr, coords, + vdw_radii, sparsity, + mouth_oracle) + + calculator.findDeepestTetrahedra(c_surface_cavities, s_clr.neigh, s_clr.verti, + coords, s_clr.simp) if start_point is not None: - calculator.set_starting_tetrahedra_from_point(c_surface_cavities, s_clr.verti, start_point) - - c_filtered_cavities = calculator.filter_cavities(c_surface_cavities, min_depth) + c_surface_cavities = calculator.setStartingTetrahedraFromPoint( + c_surface_cavities, s_clr.verti, start_point, coords, vdw_radii, + s_clr.simp, s_clr.neigh, restrict_channels_to_start_point, + start_point_search) + + c_filtered_cavities = calculator.filterCavities(c_surface_cavities, min_depth) + LOGGER.report('{0} surface cavities detected and filtered in %.2fs.'.format( + len(c_filtered_cavities)), '_prody_channels_cavities') if cavities_only: if max_depth is not None: - calculator.trim_cavities_by_depth(c_filtered_cavities, max_depth) + calculator.trimCavitiesByDepth(c_filtered_cavities, max_depth) if min_tetrahedra is not None or max_tetrahedra is not None: - c_filtered_cavities = calculator.filter_cavities_by_tetrahedra(c_filtered_cavities, min_tetrahedra, max_tetrahedra) + c_filtered_cavities = calculator.filterCavitiesByTetrahedra( + c_filtered_cavities, min_tetrahedra, max_tetrahedra) calculator.calculate_cavity_volumes(c_filtered_cavities, s_clr.simp, coords) if min_volume is not None or max_volume is not None: - c_filtered_cavities = calculator.filter_cavities_by_volume(c_filtered_cavities, min_volume, max_volume) + c_filtered_cavities = calculator.filterCavitiesByVolume( + c_filtered_cavities, min_volume, max_volume) - merged_cavities = calculator.merge_cavities(c_filtered_cavities, s_clr.simp) + merged_cavities = calculator.mergeCavities(c_filtered_cavities, s_clr.simp) # Early-return for the calcSurfaceCavities function: if cavities_only: @@ -842,21 +1419,38 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 else: LOGGER.info("Saving multiple surface cavities to directory " + str(output_path.parent) + ".") - calculator.save_cavities_to_pdb(c_filtered_cavities, s_clr.verti, output_path, separate) - + calculator.saveCavitiesToPdb(c_filtered_cavities, s_clr.verti, + output_path, separate) + + LOGGER.report('Surface cavity calculation completed in %.2fs.', '_prody_calcChannels') return c_filtered_cavities, [coords, s_srf.simp, merged_cavities, s_clr.simp, s_clr.verti] + LOGGER.timeit('_prody_channels_pathfinding') + # build the weighted adjacency matrix once for the whole cleared + # state, then run a single multi-target Dijkstra per cavity (scipy csgraph), + # instead of one heap Dijkstra per (seed, exit) pair. + simplices, neighbors, vertices = s_clr.getState() + graph = calculator.buildSparseGraph(simplices, neighbors, vertices, coords, + vdw_radii) for cavity in c_filtered_cavities: - #calculator.dijkstra(cavity, *(s_clr.get_state() + [coords, vdw_radii])) - calculator.dijkstra(cavity, *(s_clr.get_state() + tuple([coords, vdw_radii]))) - - calculator.filter_channels_by_bottleneck(c_filtered_cavities, bottleneck) + calculator.dijkstra(cavity, graph, simplices, neighbors, vertices, + coords, vdw_radii, + truncate_at_surface, similarity) + LOGGER.report('Channel pathfinding (graph Dijkstra) over {0} cavities completed in %.2fs.'.format( + len(c_filtered_cavities)), '_prody_channels_pathfinding') + + calculator.filterChannelsByBottleneck(c_filtered_cavities, bottleneck) if min_volume is not None or max_volume is not None: - calculator.filter_channels_by_volume(c_filtered_cavities, min_volume, max_volume) + calculator.filterChannelsByVolume(c_filtered_cavities, min_volume, + max_volume) channels = [channel for cavity in c_filtered_cavities for channel in cavity.channels] - + # Order channels by ascending Dijkstra cost so that channel 0 is the best + # tunnel (a short path through wide tetrahedra). This ordering drives both + # the returned list and the channel numbering in the saved PQR/PDB files. + channels.sort(key=lambda ch: ch.cost if ch.cost is not None else float('inf')) + no_of_channels = len(channels) LOGGER.info("Detected " + str(no_of_channels) + " channels.") @@ -873,55 +1467,68 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 LOGGER.info("Saving results to " + str(output_path) + ".") else: LOGGER.info("Saving multiple results to directory " + str(output_path.parent) + ".") - calculator.save_channels_to_pdb(c_filtered_cavities, output_path, separate) + calculator.saveChannelsToPdb(channels, output_path, separate) else: LOGGER.info("No output path given.") - + + LOGGER.report('Channel calculation completed in %.2fs.', '_prody_calcChannels') return channels, [coords, s_srf.simp, merged_cavities, s_clr.simp] -def calcChannelsMultipleFrames(atoms, trajectory=None, output_path=None, separate=False, start_point=None, **kwargs): - """Compute channels for each frame in a given trajectory or multi-model PDB file. - - This function calculates the channels for each frame in a trajectory or for each model - in a multi-model PDB file. The `kwargs` can include parameters necessary for channel calculation. - If the `separate` parameter is set to True, each detected channel will be saved in a separate PDB file. - - :param atoms: Atomic data or object containing atomic coordinates and methods for accessing them. +def calcChannelsMultipleFrames(atoms, trajectory=None, output_path=None, + separate=False, start_point=None, **kwargs): + """Compute channels for each frame in a given trajectory or multi-model + PDB file. + + This function calculates the channels for each frame in a trajectory or for + each model in a multi-model PDB file. The `kwargs` can include parameters + necessary for channel calculation. If the `separate` parameter is set to + True, each detected channel will be saved in a separate PDB file. + + :arg atoms: Atomic data or object containing atomic coordinates and methods + for accessing them. :type atoms: object - :param trajectory: Trajectory object containing multiple frames or a multi-model PDB file. + :arg trajectory: Trajectory object containing multiple frames or a + multi-model PDB file. :type trajectory: Atomic or Ensemble object - :param output_path: Optional path to save the resulting channels and associated data in PDB format. - If a directory is specified, each frame/model will have its results saved in separate files. - If None, results are not saved. Default is None. + :arg output_path: Optional path to save the resulting channels and + associated data in PDB format. If a directory is specified, each + frame/model will have its results saved in separate files. If None, + results are not saved. Default is None. :type output_path: str or None - :param separate: If True, each detected channel is saved to a separate PDB file for each frame/model. - If False, all channels for each frame/model are saved in a single file. Default is False. + :arg separate: If True, each detected channel is saved to a separate PDB + file for each frame/model. + If False, all channels for each frame/model are saved in a single file. + Default is False. :type separate: bool - :param start_point: Optional starting point for channel search. If provided, the algorithm will use - the tetrahedron whose Voronoi vertex is closest to this point as the starting tetrahedron (overriding - the default automatic seed selection based on the deepest tetrahedron). Coordinates must be given in Å. + :arg start_point: Optional starting point for channel search. If provided, + the algorithm will use the tetrahedron whose Voronoi vertex is closest + to this point as the starting tetrahedron (overriding the default automatic + seed selection based on the deepest tetrahedron). Coordinates must be given in Å. :type start_point: list, tuple, or ndarray (length 3), or None - :param kwargs: Additional parameters required for channel calculation. This can include parameters such as - radius values (r1, r2), minimum depth (min_depth), bottleneck values, etc. + :arg kwargs: Additional parameters required for channel calculation. This can + include parameters such as radius values (r1, r2), minimum depth (min_depth), + bottleneck values, etc. See the available parameters in calcChannels(). :type kwargs: dict - :returns: List of channels and surfaces computed for each frame or model. Each entry in the list corresponds - to a specific frame or model. + :returns: List of channels and surfaces computed for each frame or model. + Each entry in the list corresponds to a specific frame or model. :rtype: list of lists Example usage: - channels_all, surfaces_all = calcChannelsMultipleFrames(atoms, trajectory=traj, output_path="channels.pdb", - separate=False, r1=3, r2=1.25, min_depth=10, bottleneck=1, sparsity=15) + channels_all, surfaces_all = calcChannelsMultipleFrames(atoms, trajectory=traj, + output_path="channels.pdb", separate=False, r1=3, + r2=0.9, min_depth=5, bottleneck=1, sparsity=3) - channels_all, surfaces_all = calcChannelsMultipleFrames(atoms, trajectory=traj, output_path="channels.pdb", - separate=False, start_point=[-10.353, -0.133, 5.608]) """ + channels_all, surfaces_all = calcChannelsMultipleFrames(atoms, trajectory=traj, + output_path="channels.pdb", separate=False, + start_point=[-10.353, -0.133, 5.608]) """ if PY3K: @@ -996,7 +1603,8 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, output_path=None, separat return channels_all, surfaces_all -def calcSurfaceCavitiesMultipleFrames(atoms, trajectory=None, output_path=None, separate=False, **kwargs): +def calcSurfaceCavitiesMultipleFrames(atoms, trajectory=None, output_path=None, + separate=False, **kwargs): """Compute surface cavities for each frame in a trajectory or multi-model PDB. This function calculates surface cavities for each frame of a trajectory or @@ -1006,33 +1614,33 @@ def calcSurfaceCavitiesMultipleFrames(atoms, trajectory=None, output_path=None, passed directly to :func:`calcSurfaceCavities` and can include parameters controlling cavity detection, filtering, and output generation. - :param atoms: Atomic object containing the molecular structure. For trajectory + :arg atoms: Atomic object containing the molecular structure. For trajectory analysis, this object provides the reference topology and is updated with coordinates from each frame. For multi-model PDB files, the individual coordinate sets are analyzed one by one. :type atoms: :class:`.Atomic` - :param trajectory: Optional trajectory or ensemble object containing multiple + :arg trajectory: Optional trajectory or ensemble object containing multiple coordinate frames. If provided, surface cavities are calculated for each selected trajectory frame. If not provided, the function attempts to use multiple coordinate sets stored in `atoms`. :type trajectory: :class:`.Atomic`, :class:`.Ensemble`, or trajectory-like object - :param output_path: Optional filename used to save detected surface cavities. + :arg output_path: Optional filename used to save detected surface cavities. If provided, one output file is generated for each frame/model by appending the frame/model index to the file name. If `None`, results are returned but not written in the folder. Default is `None`. :type output_path: str or None - :param separate: If `True`, each detected surface cavity is saved as a separate + :arg separate: If `True`, each detected surface cavity is saved as a separate PQR/PDB file for each frame/model. If `False`, all cavities detected in a given frame/model are saved in a single file. Default is `False`. :type separate: bool - :param kwargs: Additional parameters passed to :func:`calcSurfaceCavities`. + :arg kwargs: Additional parameters passed to :func:`calcSurfaceCavities`. These can include `r1`, `r2`, `min_depth`, `max_depth`, `min_tetrahedra`, `max_tetrahedra`, `min_volume`, `max_volume`, - `sparsity`, `start_frame`, and `stop_frame`. + `start_frame`, and `stop_frame`. :type kwargs: dict :returns: Two lists: @@ -1044,10 +1652,12 @@ def calcSurfaceCavitiesMultipleFrames(atoms, trajectory=None, output_path=None, Example usage: protein = parsePDB('1tqn').select('protein') - cavities_all, surfaces_all = calcSurfaceCavitiesMultipleFrames(protein, trajectory=traj, output_path="surface_cavities", - r1=4.5, r2=2.0, min_depth=2, max_depth=3, min_volume=50) + cavities_all, surfaces_all = calcSurfaceCavitiesMultipleFrames(protein, + trajectory=traj, output_path="surface_cavities", + r1=4.5, r2=2.0, min_depth=1.5, max_depth=2.5, min_volume=50) - cavities_all, surfaces_all = calcSurfaceCavitiesMultipleFrames(protein, start_frame=0, stop_frame=10, r1=4.5, r2=2.0) """ + cavities_all, surfaces_all = calcSurfaceCavitiesMultipleFrames(protein, start_frame=0, + stop_frame=10, r1=4.5, r2=2.0) """ if PY3K: if not checkAndImport('pathlib'): @@ -1095,9 +1705,12 @@ def calcSurfaceCavitiesMultipleFrames(atoms, trajectory=None, output_path=None, if output_path: cavities, surface = calcSurfaceCavities(atoms_copy, - output_path=str(output_path) + "{0}.pqr".format(j0), separate=separate, **kwargs) + output_path=str(output_path) + "{0}.pqr".format(j0), + separate=separate, **kwargs) else: - cavities, surface = calcSurfaceCavities(atoms_copy, separate=separate, **kwargs) + cavities, surface = calcSurfaceCavities(atoms_copy, + separate=separate, + **kwargs) cavities_all.append(cavities) surfaces_all.append(surface) @@ -1138,7 +1751,8 @@ def calcSurfaceCavitiesMultipleFrames(atoms, trajectory=None, output_path=None, def parseParameters(channels, **kwargs): - """Extracts and returns the lengths, bottlenecks, and volumes of each channel in a given list of channels. """ + """Extracts and returns the lengths, bottlenecks, and volumes of each + channel in a given list of channels. """ lengths = [] bottlenecks = [] @@ -1158,24 +1772,28 @@ def parseParameters(channels, **kwargs): def getChannelParameters(channels, **kwargs): - """Extracts and returns the lengths, bottlenecks, and volumes of each channel in a given list of channels. + """Extracts and returns the lengths, bottlenecks, and volumes of each + channel in a given list of channels. - This functaaion iterates through a list of channel objects, extracting the length, bottleneck, - and volume of each channel. These values are collected into separate lists, which are returned - as a tuple for further use. + This function iterates through a list of channel objects, extracting the + length, bottleneck, and volume of each channel. These values are collected + into separate lists, which are returned as a tuple for further use. - :arg channels: A list of channel objects, where each channel has attributes `length`, `bottleneck`, - and `volume`. These attributes represent the length of the channel, the minimum radius - (bottleneck) along its path, and the total volume of the channel, respectively. + :arg channels: A list of channel objects, where each channel has attributes + `length`, `bottleneck`,and `volume`. These attributes represent the + length of the channel, the minimum radius (bottleneck) along its path, + and the total volume of the channel, respectively. :type channels: list - :arg param_file_name: The files with parameters will be saved in a text file with the provided name. - Use one word which will be added to '_Parameters_All_channels.txt' sufix. - If further analysis will be performed with selectChannelBySelection() function, the preferable + :arg param_file_name: The files with parameters will be saved in a text + file with the provided name. Use one word which will be added to + '_Parameters_All_channels.txt' suffix. If further analysis will be + performed with selectChannelBySelection() function, the preferable param_file_name is PDB+chain for example: '1bbhA'. :type param_file_name: str - :returns: Three lists containing the lengths, bottlenecks, and volumes of the channels. + :returns: Three lists containing the lengths, bottlenecks, and volumes of + the channels. :rtype: tuple (list, list, list) Example usage: @@ -1187,7 +1805,9 @@ def getChannelParameters(channels, **kwargs): try: results_L_B_V = parseParameters(channels, **kwargs) lengths, bottlenecks, volumes = results_L_B_V - LOGGER.info("Channel {0}: \t{1} \t{2} \t{3}".format('ID', 'Volume [ų]', 'Length [Å]', 'Bottleneck [Å]')) + LOGGER.info("Channel {0}: \t{1} \t{2} \t{3}".format('ID', 'Volume [ų]', + 'Length [Å]', + 'Bottleneck [Å]')) for i in range(len(lengths)): LOGGER.info("channel {0}: \t{1} \t\t{2} \t\t{3}".format(i, np.round(volumes[i],2), np.round(lengths[i], 2), np.round(bottlenecks[i], 2))) return results_L_B_V @@ -1198,7 +1818,9 @@ def getChannelParameters(channels, **kwargs): results = parseParameters(channels[nr_i], param_file_name=safe_param_file_name + str(nr_i)) multi_model_param.append(results) - LOGGER.info("Channel {0}: \t{1} \t{2} \t{3}".format('ID', 'Volume [ų]', 'Length [Å]', 'Bottleneck [Å]')) + LOGGER.info("Channel {0}: \t{1} \t{2} \t{3}".format('ID', 'Volume [ų]', + 'Length [Å]', + 'Bottleneck [Å]')) for frame_nr, frame in enumerate(multi_model_param): lengths, bottlenecks, volumes = frame LOGGER.info("Frame {0}".format(frame_nr)) @@ -1262,7 +1884,7 @@ def parseSurfaceCavityParameters(cavities, **kwargs): tetrahedra_counts.append(tetrahedra_count) if param_file_name is not None: - lines.append("{0}_cavity{1}: {2:.3f} {3} {4}\n".format( + lines.append("{0}_cavity{1}: {2:.3f} {3:.2f} {4}\n".format( param_file_name, nr_cav, volume, depth, tetrahedra_count)) if param_file_name is not None: @@ -1300,7 +1922,9 @@ def getSurfaceCavityParameters(cavities, **kwargs): results_V_D_T = parseSurfaceCavityParameters(cavities, **kwargs) volumes, depths, tetrahedra_counts = results_V_D_T - LOGGER.info("Cavity {0}: \t{1} \t{2} \t{3}".format('ID', 'Volume [ų]', 'Depth [Å]', 'Tetrahedra count')) + LOGGER.info("Cavity {0}: \t{1} \t{2} \t{3}".format('ID', 'Volume [ų]', + 'Depth [Å]', + 'Tetrahedra count')) for i in range(len(volumes)): LOGGER.info("cavity {0}: \t{1} \t\t{2} \t\t{3}".format(i, np.round(volumes[i], 2), np.round(depths[i], 2), @@ -1315,7 +1939,9 @@ def getSurfaceCavityParameters(cavities, **kwargs): param_file_name=safe_param_file_name + str(nr_i)) multi_model_param.append(results) - LOGGER.info("Cavity {0}: \t{1} \t{2} \t{3}".format('ID', 'Volume [ų]', 'Depth [Å]', 'Tetrahedra count')) + LOGGER.info("Cavity {0}: \t{1} \t{2} \t{3}".format('ID', 'Volume [ų]', + 'Depth [Å]', + 'Tetrahedra count')) for frame_nr, frame in enumerate(multi_model_param): volumes, depths, tetrahedra_counts = frame @@ -1357,30 +1983,35 @@ def getSurfaceCavityParametersMultipleFrames(cavities_all, **kwargs): def getChannelAtoms(channels, protein=None, num_samples=5): - """Generates an AtomGroup object representing the atoms along the paths of the given channels - and optionally combines them with an existing protein structure. - - This function takes a list of channel objects and generates atomic representations of the - channels based on their centerline splines and radius splines. The function samples points - along each channel's centerline and assigns atom positions at these points with corresponding - radii, creating a list of PDB-formatted lines. These lines are then converted into an AtomGroup - object using the ProDy library. If a protein structure is provided, it is combined with the - generated channel atoms by merging their respective PDB streams. - - :param channels: A list of channel objects. Each channel has a method `get_splines()` that + """Generates an AtomGroup object representing the atoms along the paths of + the given channels and optionally combines them with an existing protein + structure. + + This function takes a list of channel objects and generates atomic + representations of the channels based on their centerline splines and + radius splines. The function samples points along each channel's centerline + and assigns atom positions at these points with corresponding radii, + creating a list of PDB-formatted lines. These lines are then converted + into an AtomGroup object using the ProDy library. If a protein structure is + provided, it is combined with the generated channel atoms by merging their + respective PDB streams. + + :arg channels: A list of channel objects. Each channel has a method + `getSplines()` that returns the centerline spline and radius spline of the channel. :type channels: list - :param protein: An optional AtomGroup object representing a protein structure. If provided, - it will be combined with the generated channel atoms. + :arg protein: An optional AtomGroup object representing a protein structure. + If provided, it will be combined with the generated channel atoms. :type protein: prody.atomic.AtomGroup or None - :param num_samples: The number of atom samples to generate along each segment of the channel. - More samples result in a finer representation of the channel. Default is 5. + :arg num_samples: The number of atom samples to generate along each segment + of the channel. More samples result in a finer representation of the + channel. Default is 5. :type num_samples: int - :returns: An AtomGroup object representing the combined atoms of the channels and the protein, - if a protein is provided. + :returns: An AtomGroup object representing the combined atoms of the + channels and the protein, if a protein is provided. :rtype: prody.atomic.AtomGroup Example usage: @@ -1406,13 +2037,15 @@ def convert_lines_to_atomic(atom_lines): channels = [channels] for channel in channels: - centerline_spline, radius_spline = channel.get_splines() + centerline_spline, radius_spline = channel.getSplines() samples = len(channel.tetrahedra) * num_samples t = np.linspace(centerline_spline.x[0], centerline_spline.x[-1], samples) centers = centerline_spline(t) radii = radius_spline(t) - for i, (x, y, z, radius) in enumerate(zip(centers[:, 0], centers[:, 1], centers[:, 2], radii), start=atom_index): + for i, (x, y, z, radius) in enumerate(zip(centers[:, 0], centers[:, 1], + centers[:, 2], radii), + start=atom_index): pdb_lines.append("ATOM %5d H FIL T 1 %8.3f%8.3f%8.3f%6.2f%6.2f\n" % (i, x, y, z, 1.00, radius)) if protein is not None: @@ -1444,17 +2077,19 @@ def getChannelResidueNames(atoms, channels, **kwargs): :arg atoms: an Atomic object from which residues are selected :type atoms: :class:`.Atomic`, :class:`.LigandInteractionsTrajectory` - :param channels: A list of channel objects. Each channel has a method `get_splines()` that - returns the centerline spline and radius spline of the channel. + :arg channels: A list of channel objects. Each channel has a method + `getSplines()` that returns the centerline spline and radius spline of + the channel. :type channels: list :arg distA: Residues will be provided based on this value. default is 4 [Ang] :type distA: int, float - :arg residues_file_name: The file with residues will be saved in a text file with the provided name. - Use one word which will be added to '_Residues_All_channels.txt' sufix. - If further analysis will be performed with selectChannelBySelection() function, the preferable + :arg residues_file_name: The file with residues will be saved in a text + file with the provided name. Use one word which will be added to + '_Residues_All_channels.txt' sufix. If further analysis will be + performed with selectChannelBySelection() function, the preferable residues_file_name is PDB+chain for example: '1bbhA'. :type residues_file_name: str @@ -1659,7 +2294,8 @@ def getSurfaceCavityResidueNames(atoms, cavities, surface, **kwargs): :type distA: int, float :arg residues_file_name: The file with residues will be saved in a text file - with the provided name. The suffix '_Residues_All_surface_cavities.txt' will be added. + with the provided name. The suffix '_Residues_All_surface_cavities.txt' + will be added. :type residues_file_name: str :arg one_letter_aa: Whether to apply one-letter code to residue names. @@ -1739,15 +2375,20 @@ def getSurfaceCavityResidueNames(atoms, cavities, surface, **kwargs): return selected_residues_cav -def getSurfaceCavityResidueNamesMultipleFrames(atoms, cavities_all, surfaces_all, trajectory=None, **kwargs): - """Provides residue names for surface cavities calculated for multiple frames/models. +def getSurfaceCavityResidueNamesMultipleFrames(atoms, cavities_all, + surfaces_all, + trajectory=None, **kwargs): + """Provides residue names for surface cavities calculated for multiple + frames/models. This function is a multi-frame wrapper for :func:`getSurfaceCavityResidueNames`. - For each model or trajectory frame, the atomic coordinates are matched with the - corresponding surface cavity prediction. Thus, cavities calculated for frame/model ``i`` - are analyzed against the protein coordinates from frame/model ``i``. + For each model or trajectory frame, the atomic coordinates are matched with + the corresponding surface cavity prediction. Thus, cavities calculated + for frame/model ``i`` are analyzed against the protein coordinates from + frame/model ``i``. - This function should be used with results returned by :func:`calcSurfaceCavitiesMultipleFrames`. + This function should be used with results returned by + :func:`calcSurfaceCavitiesMultipleFrames`. :arg atoms: an Atomic object from which residues are selected. :type atoms: :class:`.Atomic` @@ -1794,7 +2435,8 @@ def getSurfaceCavityResidueNamesMultipleFrames(atoms, cavities_all, surfaces_all if trajectory is None: # multi-model PDB - for frame_pos, (cavities, surface) in enumerate(zip(cavities_all, surfaces_all)): + for frame_pos, (cavities, surface) in enumerate(zip(cavities_all, + surfaces_all)): model_index = start_frame + frame_pos atoms.setACSIndex(model_index) @@ -1832,8 +2474,8 @@ def getSurfaceCavityResidueNamesMultipleFrames(atoms, cavities_all, surfaces_all def selectChannelBySelection(atoms, residue_sele, **kwargs): - """Select PQR files with channels that are having FIL residues within certain distance (distA) from - selected residue (temporarly one residue). + """Select PQR files with channels that are having FIL residues within + certain distance (distA) from selected residue (temporarily one residue). If not all files should be included use pqr_files to provide the new list. For example: pqr_files = [file for file in os.listdir('.') if file.startswith('7lafA_') and file.endswith('.pqr')] @@ -1847,22 +2489,23 @@ def selectChannelBySelection(atoms, residue_sele, **kwargs): :type residue_sele: str :arg pqr_files: list of PQR files to analyze - default is False, which means that all .pqr files from the current directory will be analyzed. + default is False, which means that all .pqr files from the + current directory will be analyzed. :type pqr_files: bool or list :arg folder_name: The name of the folder to which PDBs will be extracted :type folder_name: str - :arg distA: non-zero value, maximal distance from selected region to channel (FIL atoms) - default is 5 + :arg distA: non-zero value, maximal distance from selected region to + channel (FIL atoms). Default is 5. :type distA: int, float - :arg residues_file: File with residues forming the channel created by getChannelResidues() - default is False + :arg residues_file: File with residues forming the channel created by + getChannelResidues(), default is False :type residues_file: bool - :arg param_file: File with residues forming the channel created by getChannelParameters() - default is False + :arg param_file: File with residues forming the channel created by + getChannelParameters(). Default is False. :type param_file: bool """ try: @@ -1876,7 +2519,6 @@ def selectChannelBySelection(atoms, residue_sele, **kwargs): 'with `getCoords` method') import os, shutil - import numpy as np pqr_files = kwargs.pop('pqr_files', False) distA = kwargs.pop('distA', 5) @@ -1887,8 +2529,10 @@ def selectChannelBySelection(atoms, residue_sele, **kwargs): object_name = kwargs.pop('object_name', 'channel') residues_suffix = kwargs.pop('residues_suffix', '_Residues_All_channels.txt') parameters_suffix = kwargs.pop('parameters_suffix', '_Parameters_All_channels.txt') - selected_residues_output = kwargs.pop('selected_residues_output', 'Selected_channel_residues.txt') - selected_parameters_output = kwargs.pop('selected_parameters_output', 'Selected_channel_parameters.txt') + selected_residues_output = kwargs.pop('selected_residues_output', + 'Selected_channel_residues.txt') + selected_parameters_output = kwargs.pop('selected_parameters_output', + 'Selected_channel_parameters.txt') copied_files_list = [] @@ -1912,7 +2556,7 @@ def selectChannelBySelection(atoms, residue_sele, **kwargs): else: pass - # Extract paramaters and/or residues with channel selection + # Extract parameters and/or residues with channel selection if residues_file == True: selected_residues = [] for file in copied_files_list: @@ -1998,43 +2642,63 @@ def selectSurfaceCavityBySelection(atoms, residue_sele, **kwargs): kwargs.setdefault('object_name', 'cavity') kwargs.setdefault('residues_suffix', '_Residues_All_surface_cavities.txt') kwargs.setdefault('parameters_suffix', '_Parameters_All_surface_cavities.txt') - kwargs.setdefault('selected_residues_output', 'Selected_surface_cavity_residues.txt') - kwargs.setdefault('selected_parameters_output', 'Selected_surface_cavity_parameters.txt') + kwargs.setdefault('selected_residues_output', + 'Selected_surface_cavity_residues.txt') + kwargs.setdefault('selected_parameters_output', + 'Selected_surface_cavity_parameters.txt') return selectChannelBySelection(atoms, residue_sele, **kwargs) def calcChannelSurfaceOverlaps(**kwargs): - """Calculate overlapping parts of the predicted channels, tunnels, and pores denote as 'FIL' atoms. - Results are normalized within [0,1]. + """Calculate overlapping parts of the predicted channels, tunnels, and + pores denote as 'FIL' atoms. Results are normalized within [0,1]. :arg resolution: Surface sampling resolution. default is 0.5 :type resolution: float + + :arg max_proc: Maximum number of parallel processes used to voxelize individual + PQR files. If 1, files are processed serially. If None, all available CPU + cores are used. Default is 2. + :type max_proc: int or None :arg output_file_name: The name of the PDB file with overlapping surfaces. :type output_file_name: str - :arg pqr_files: File with residues forming the channel created by getChannelResidues() - default is False (then all the files from the current directory will be analyzed) - when providing a list, only the PDBs from list will be analyzed - when providing str, it will be treated as a folder path + :arg pqr_files: File with residues forming the channel created by + getChannelResidues(). Default is False, then all the files from the + current directory will be analyzed. When providing a list, only the + PDBs from the list will be analyzed. When providing str, it will be + treated as a folder path. :type pqr_files: bool, list or str Example usage: - calcChannelSurfaceOverlaps() - all the files in the current directory will be analyzed + calcChannelSurfaceOverlaps() - all the files in the current directory will + be analyzed + + from pathlib import Path + pqr_files = [str(f) for f in Path(".").glob("channels_*.pqr")] + calcChannelSurfaceOverlaps(pqr_files=pqr_files, + output_file_name='results.pdb', max_proc=4) + - files with the "channels_" prefix will be selected from the current folder + and analyzed using four parallel processes. - calcChannelSurfaceOverlaps(pqr_files='./DATA', output_file_name='results.pdb') - only files from - the DATA folder will be analyzed and results will be saved as results.pdb + calcChannelSurfaceOverlaps(pqr_files='./DATA', output_file_name='results.pdb') + - only files from the DATA folder will be analyzed and results will be saved + as results.pdb list_of_files = ['file1.pqr', 'file2.pqr', 'file3.pqr', ..] - calcChannelSurfaceOverlaps(pqr_files=list_of_files, output_file_name='results.pdb') - files from - the list will be analyzed and results will be saved as results.pdb + calcChannelSurfaceOverlaps(pqr_files=list_of_files, output_file_name='results.pdb') + - files from the list will be analyzed and results will be saved as results.pdb """ import os + import multiprocessing + from collections import Counter resolution = kwargs.pop('resolution', 0.5) + max_proc = kwargs.pop('max_proc', 2) pqr_files = kwargs.pop('pqr_files', False) if pqr_files == False or pqr_files is None: @@ -2042,7 +2706,8 @@ def calcChannelSurfaceOverlaps(**kwargs): pqr_files = [file for file in os.listdir('.') if file.endswith('.pqr')] elif isinstance(pqr_files, str): # folder path - pqr_files = [file for file in os.listdir(pqr_files) if file.endswith('.pqr')] + if os.path.isdir(pqr_files): + pqr_files = [os.path.join(pqr_files, file) for file in os.listdir(pqr_files) if file.endswith('.pqr')] elif isinstance(pqr_files, list): # list of PQRs pqr_files = [file for file in pqr_files if file.endswith('.pqr')] @@ -2050,71 +2715,64 @@ def calcChannelSurfaceOverlaps(**kwargs): raise ValueError('Please provide list with PQR files, folder path, or nothing to analyze PQRs in the current folder') output_file_name = kwargs.pop('output_file_name','overlap_regions.pdb') + + if len(pqr_files) == 0: + LOGGER.info("No PQR files found.") + return None + if os.path.exists(output_file_name): - os.rename(output_file_name, output_file_name+'-old') + os.rename(output_file_name, output_file_name + '-old') - def loadPDBdata(filepath): - """Parse a PQR file and return a list of atom dictionaries for lines containing 'FIL'.""" - atoms_set = [] - FILatoms = parsePQR(filepath).select('resname FIL') - - if FILatoms == None: - pass - else: - for nr_i, i in enumerate(FILatoms): - FILatoms_coords = FILatoms.getCoords()[nr_i] - FILBetas_value = FILatoms.getRadii()[nr_i] - atoms_set.append({ - 'x': float(FILatoms_coords[0]), - 'y': float(FILatoms_coords[1]), - 'z': float(FILatoms_coords[2]), - 'radius': float(FILBetas_value) - }) - return atoms_set - - def create_surface(atoms, resolution=resolution): - """Create a 3D grid representing the surface occupied by the atoms.""" - surface = {} - Zr = 0 - for atom in atoms: - x, y, z, radius = atom['x'], atom['y'], atom['z'], atom['radius'] - for i in np.arange(x - radius, x + radius, resolution): - for j in np.arange(y - radius, y + radius, resolution): - for k in np.arange(z - radius, z + radius, resolution): - if (i - x) ** 2 + (j - y) ** 2 + (k - z) ** 2 <= radius ** 2: - key = (round(i, Zr), round(j, Zr), round(k, Zr)) - surface[key] = surface.get(key, 0) + 1 - return surface - - def merge_surfaces(surfaces): - """Merge multiple surfaces and calculate overlap counts.""" - merged_surface = {} - for surface in surfaces: - for key in surface: - merged_surface[key] = merged_surface.get(key, 0) + 1 - return merged_surface - - def write_merge_surf_pdb(merged_surface, filename, nr_pdbs): - """Write the merged surface into a PDB file.""" - with open(filename, 'w') as file: - atom_id = 1 - for (x, y, z), count in merged_surface.items(): - norm_count = count/nr_pdbs - file.write("ATOM {:5d} H FIL T 1 {:8.3f}{:8.3f}{:8.3f}{:6.2f} 1.00\n".format(atom_id, x, y, z, norm_count)) - atom_id += 1 - - surfaces = [] - for nr_pdbs,pqr_file in enumerate(pqr_files): - LOGGER.info("Processing file: {0}".format(pqr_file)) - atoms = loadPDBdata(pqr_file) - if atoms: - surface = create_surface(atoms, resolution=resolution) - surfaces.append(surface) - - nr_pdbs = nr_pdbs+1 - merged_surface = merge_surfaces(surfaces) - write_merge_surf_pdb(merged_surface, output_file_name, nr_pdbs) + if max_proc is None: + max_proc = multiprocessing.cpu_count() + + max_proc = int(max_proc) + if max_proc < 1: + max_proc = 1 + max_proc = min(max_proc, len(pqr_files)) + + LOGGER.info("Number of PQR files: {0}".format(len(pqr_files))) + LOGGER.info("Resolution: {0}".format(resolution)) + LOGGER.info("max_proc: {0}".format(max_proc)) + + merged_surface = Counter() + tasks = [(pqr_file, resolution) for pqr_file in pqr_files] + + if max_proc > 1: + LOGGER.info("Calculating overlaps using {0} processes.".format(max_proc)) + chunksize = max(1, len(tasks) // (max_proc * 4)) + with multiprocessing.Pool(processes=max_proc) as pool: + for surface in pool.imap_unordered(_surfaceFromPqrWorker, tasks, + chunksize=chunksize): + merged_surface.update(surface) + + else: + for pqr_file in pqr_files: + LOGGER.info("Processing file: {0}".format(pqr_file)) + surface = _surfaceFromPqrWorker((pqr_file, resolution)) + merged_surface.update(surface) + + with open(output_file_name, 'w') as out: + atom_id = 1 + + for (ix, iy, iz), count in merged_surface.items(): + x = ix * resolution + y = iy * resolution + z = iz * resolution + + norm_count = float(count) / float(len(pqr_files)) + + out.write("ATOM {:5d} H FIL T 1 {:8.3f}{:8.3f}{:8.3f}{:6.2f} 1.00\n" + .format(atom_id, x, y, z, norm_count)) + + atom_id += 1 + + LOGGER.info("Overlap written to: {0}".format(output_file_name)) + LOGGER.info("Number of occupied overlap voxels: {0}".format(len(merged_surface))) + + return output_file_name + def calcSurfaceCavityOverlaps(**kwargs): """Calculate overlapping regions of surface cavities represented as FIL atoms. @@ -2128,6 +2786,11 @@ def calcSurfaceCavityOverlaps(**kwargs): :arg resolution: surface sampling resolution. Default is 0.5. :type resolution: float + + :arg max_proc: Maximum number of parallel processes used to voxelize individual + PQR files. If 1, files are processed serially. If None, all available CPU + cores are used. Default is 2. + :type max_proc: int or None :arg output_file_name: name of the output PDB file with overlapping cavity regions. Default is ``'surface_cavity_overlap_regions.pdb'``. @@ -2143,73 +2806,94 @@ def calcSurfaceCavityOverlaps(**kwargs): return calcChannelSurfaceOverlaps(**kwargs) -def calcSurfaceCavities(atoms, output_path=None, r1=4.5, r2=2.0, min_depth=2, max_depth=3, - min_tetrahedra=None, max_tetrahedra=None, min_volume=50, max_volume=None, - sparsity=15, separate=False): - """Calculate surface cavities (pockets) on protein surface using CaviTracer approach. +def calcSurfaceCavities(atoms, output_path=None, r1=4.5, r2=2.0, min_depth=1.5, + max_depth=2.5, min_tetrahedra=None, max_tetrahedra=None, + min_volume=50, max_volume=None, sparsity=None, + separate=False): + """Calculate surface cavities (pockets) on protein surface using CaviTracer + approach. - :param atoms: An object representing the molecular structure, typically containing atomic coordinates - and element types. + :arg atoms: An object representing the molecular structure, typically + containing atomic coordinates and element types. :type atoms: `Atoms` object - :param output_path: Optional path to save the resulting cavities and associated data in PQR (or PDB) format. - If None, results are not saved. Default is None. + :arg output_path: Optional path to save the resulting cavities and + associated data in PQR (or PDB) format. If None, results are not saved. + Default is None. :type output_path: str or None - :param separate: If True, each detected cavity is saved to a separate PQR file. If False, all cavities - are saved in a single PQR file. Default is False. + :arg separate: If True, each detected cavity is saved to a separate PQR + file. If False, all cavities are saved in a single PQR file. Default is + False. :type separate: bool - :param r1: The first radius threshold used during the deletion of simplices, which is used to define - the outer surface of the cavities. Default is 4.5. + :arg r1: The first radius threshold used during the deletion of simplices, + which is used to define the outer surface of the cavities. Default is 4.5. :type r1: float - :param r2: The second radius threshold used to define the inner surface of the cavities. Default is 2. + :arg r2: The second radius threshold used to define the inner surface of + the cavities. Default is 2. :type r2: float - :param min_depth: The minimum depth a cavity must have to be considered as a cavity. Default is 2. - :type min_depth: int - - :param max_depth: Maximum cavity depth. Cavities deeper than this value are trimmed to the specified depth. - Default is 3. - :type max_depth: int - - :param sparsity: The sparsity parameter controls the sampling density when analyzing the molecular surface. - A higher value results in fewer sampling points. Default is 15. + :arg min_depth: The minimum depth, in Angstrom, a cavity must reach to be + considered. Depth is the geodesic distance from the surface opening along + the Voronoi network, a physical length independent of tessellation density. + Default is 1.5. + :type min_depth: float + + :arg max_depth: Maximum cavity depth, in Angstrom. Portions of a cavity deeper + than this value are trimmed away, keeping the shallow surface shell that + defines a pocket. Default is 2.5. + :type max_depth: float + + :arg sparsity: Deprecated and ignored; accepted only so that existing calls + keep working. It never affected surface cavities: it thinned the sampling + of the mouth (exit) tetrahedra used as termini by the *channel* search, + and no cavity property reads that thinned set. Cavity extent, depth, + volume and filtering are all derived from the unthinned exit tetrahedra, + so passing 1 or 15 returns the same cavities. :type sparsity: int - :param min_tetrahedra: Minimum number of tetrahedra required for a cavity to be retained. - Smaller cavities are discarded. Default is None. + :arg min_tetrahedra: Minimum number of tetrahedra required for a cavity to + be retained. Smaller cavities are discarded. Default is None. :type min_tetrahedra: int - :param max_tetrahedra: Maximum number of tetrahedra allowed for a cavity to be retained. - Larger cavities are discarded. Default is None. + :arg max_tetrahedra: Maximum number of tetrahedra allowed for a cavity to + be retained. Larger cavities are discarded. Default is None. :type max_tetrahedra: int - :param min_volume: Minimum volume required for a channel/cavity to be retained. Default is 50. + :arg min_volume: Minimum volume required for a channel/cavity to be + retained. Default is 50. :type min_volume: float - :param max_volume: Maximum volume allowed for a channel/cavity to be retained. Default is None. + :arg max_volume: Maximum volume allowed for a channel/cavity to be + retained. Default is None. :type max_volume: float :returns: A tuple containing two elements: - - `cavities`: A list of detected cavities, where each channel is an object containing information - about its path and geometry. - - `surface`: A list containing additional information for further visualization, including - the atomic coordinates, simplices defining the surface, and merged cavities. + - `cavities`: A list of detected cavities, where each channel is an + object containing information about its path and geometry. + - `surface`: A list containing additional information for further + visualization, including the atomic coordinates, simplices defining + the surface, and merged cavities. :rtype: tuple (list, list) This function performs the following steps: - 1. **Selection and Filtering:** Selects non-hetero atoms from the protein, calculates van der Waals radii, - and performs 3D Delaunay triangulation and Voronoi tessellation on the coordinates. - 2. **Surface and Interior Filtering:** Iteratively removes simplices based on the user-defined radii - (`r1` and `r2`) to distinguish the molecular surface from the internal void space. - 3. **Surface Cavity Identification:** Detects connected void regions and identifies those that remain - connected to the protein surface, corresponding to surface-accessible cavities and pockets. - 4. **Depth Calculation and Filtering:** Estimates cavity depth using a graph-based traversal from the cavity - openings, identifies the deepest tetrahedra, and filters cavities according to the specified depth criteria. - 5. **Output Generation:** Optionally trims cavities exceeding the specified maximum depth, saves detected - cavities to PDB/PQR files, and returns cavity objects together with the surface representation for further + 1. **Selection and Filtering:** Selects non-hetero atoms from the protein, + calculates van der Waals radii, and performs 3D Delaunay triangulation + and Voronoi tessellation on the coordinates. + 2. **Surface and Interior Filtering:** Iteratively removes simplices based + on the user-defined radii (`r1` and `r2`) to distinguish the molecular + surface from the internal void space. + 3. **Surface Cavity Identification:** Detects connected void regions and + identifies those that remain connected to the protein surface, + corresponding to surface-accessible cavities and pockets. + 4. **Depth Calculation and Filtering:** Estimates cavity depth using a + graph-based traversal from the cavity openings, identifies the deepest + tetrahedra, and filters cavities according to the specified depth criteria. + 5. **Output Generation:** Optionally trims cavities exceeding the specified + maximum depth, saves detected cavities to PDB/PQR files, and returns + cavity objects together with the surface representation for further analysis and visualization. Example usage: @@ -2217,30 +2901,58 @@ def calcSurfaceCavities(atoms, output_path=None, r1=4.5, r2=2.0, min_depth=2, ma protein = p.select('protein') cavities, surface = calcSurfaceCavities(protein, output_path='test_surf_cav.pqr') """ - + if sparsity is not None: + _warn("sparsity is deprecated in calcSurfaceCavities and is " + "ignored. It thinned the mouth tetrahedra sampled as termini " + "by the channel search; cavities are built from the unthinned " + "ones, so it never changed them.") + + # No peel (min_enclosure=0). The enclosure peel strips the shell of true + # exterior that a large r1 probe bridges over instead of entering, because it + # offers a channel wide, low-cost routes along the outside of the protein. A + # surface cavity *is* that shell: a pocket is shallow and open by definition, + # so the peel deletes these cavities cavities, surface = calcChannels( - atoms, + atoms, output_path=output_path, separate=separate, - r1=r1, r2=r2, + r1=r1, r2=r2, min_depth=min_depth, max_depth=max_depth, min_volume=min_volume, max_volume=max_volume, min_tetrahedra=min_tetrahedra, max_tetrahedra=max_tetrahedra, - sparsity=sparsity, cavities_only=True) + min_enclosure=0.0, cavities_only=True) return cavities, surface class Channel: - def __init__(self, tetrahedra, centerline_spline, radius_spline, length, bottleneck, volume): + def __init__(self, tetrahedra, centerline_spline, radius_spline, length, + bottleneck, volume, cost=None): self.tetrahedra = tetrahedra self.centerline_spline = centerline_spline self.radius_spline = radius_spline self.length = length self.bottleneck = bottleneck self.volume = volume - - def get_splines(self): + # cost: accumulated Dijkstra path weight (sum of l / (d**2 + b) edge + # costs) from the seed to the exit. Lower is better - a short path + # through wide tetrahedra. Set by dijkstra(); None when not computed. + self.cost = cost + # curvature: path length / straight-line end-to-end distance + # (dimensionless, >= 1; 1.0 == perfectly straight). + self.curvature = self._computeCurvature() + + def _computeCurvature(self): + """Path length divided by straight-line end-to-end distance.""" + x = self.centerline_spline.x + start = np.asarray(self.centerline_spline(x[0])) + end = np.asarray(self.centerline_spline(x[-1])) + straight = float(np.linalg.norm(end - start)) + if straight <= 1e-9: + return float('nan') + return float(self.length / straight) + + def getSplines(self): return self.centerline_spline, self.radius_spline @@ -2257,12 +2969,12 @@ def __eq__(self, other): np.array_equal(self.neigh, other.neigh) and np.array_equal(self.verti, other.verti)) - def set_state(self, simplices, neighbors, vertices): + def setState(self, simplices, neighbors, vertices): self.simp = simplices self.neigh = neighbors self.verti = vertices - def get_state(self): + def getState(self): return self.simp, self.neigh, self.verti class Cavity: @@ -2275,98 +2987,413 @@ def __init__(self, tetrahedra, is_connected_to_surface): self.tetrahedra_depths = {} self.volume = 0.0 - def make_surface(self): + def makeSurface(self): self.is_connected_to_surface = True - def set_exit_tetrahedra(self, exit_tetrahedra, end_tetrahedra): + def setExitTetrahedra(self, exit_tetrahedra, end_tetrahedra): self.exit_tetrahedra = exit_tetrahedra self.end_tetrahedra = end_tetrahedra - def set_starting_tetrahedron(self, tetrahedron): + def setStartingTetrahedron(self, tetrahedron): self.starting_tetrahedron = tetrahedron - def set_depth(self, depth): + def setDepth(self, depth): self.depth = depth - def add_channel(self, channel): + def addChannel(self, channel): self.channels.append(channel) +def _rowsIsin(a, b): + """Boolean mask marking which rows of 2D integer array ``a`` occur as a row + in 2D array ``b`` (exact, order-sensitive match). + + Uses a void-dtype view so each row is treated as a single hashable scalar, + turning an O(len(a) x len(b)) row-by-row scan into an O(len(a) + len(b)) + hashed membership test. + """ + a = np.ascontiguousarray(a) + b = np.ascontiguousarray(b) + if a.size == 0 or b.size == 0: + return np.zeros(a.shape[0], dtype=bool) + if a.dtype != b.dtype: + b = b.astype(a.dtype) + va = a.view(np.dtype((np.void, a.dtype.itemsize * a.shape[1]))).ravel() + vb = b.view(np.dtype((np.void, b.dtype.itemsize * b.shape[1]))).ravel() + return np.isin(va, vb) + + class ChannelCalculator: - def __init__(self, atoms, r1=3, r2=1.25, min_depth=10, bottleneck=1, sparsity=15): + def __init__(self, atoms, r2=0.9, sparsity=1, route_tolerance=1.0, + edge_cost='integral'): + # Only the parameters the class actually consults are held here. r1, + # min_depth and bottleneck are stages of the pipeline, applied to the + # tessellation and to the finished channels by calcChannels; keeping copies + # of them on the calculator suggested it filtered by them, which it does + # not. self.atoms = atoms - self.r1 = r1 self.r2 = r2 - self.min_depth = min_depth - self.bottleneck = bottleneck self.sparsity = sparsity - - def sphere_fit(self, vertices, tetrahedron, vertice, vdw_radii, r): - center = vertice - d_sum = sum(np.linalg.norm(center - vertices[atom]) for atom in tetrahedron) - r_sum = sum(r + vdw_radii[atom] for atom in tetrahedron) - - return d_sum >= r_sum + self.route_tolerance = route_tolerance + # 'integral' (clearance-profile integral) or 'bottleneck' (l/(d^2+b)); + # the Dijkstra edge weight in buildSparseGraph. Resolved per diagram by + # calcChannels (weighted defaults to 'bottleneck'). + self.edge_cost = edge_cost + # Filled once by buildSparseGraph and read by the channel geometry: + # the per-simplex Voronoi-vertex clearance (the spline knots) and the + # per-edge gate clearance on each shared Delaunay face (the reported + # bottleneck and, later, the Dijkstra cost). Cached so the two consumers + # share one definition of width instead of recomputing it apart. + self._vertex_clearance = None + self._edge_bottleneck = None + + def sphereFit(self, points, simplices, vertices, vdw_radii, r, rows=None): + """Sum-based clearance test: for each tetrahedron, decide whether a probe + of radius ``r`` fits at its Voronoi vertex. + + Compares the sum of the distances from the Voronoi vertex to the four + atom centres against the sum of ``r + vdw_radius`` over the same four + atoms. Summing over the four atoms instead of testing each one is the + tangent-sphere test exactly when the vertex is equidistant from all four, + and a mild relaxation of it otherwise; the erosion defaults are + calibrated against that behaviour. + + :arg points: coordinates of all atoms, shape ``(n_atoms, 3)``. + :type points: :class:`~numpy.ndarray` + + :arg simplices: atom indices of each tetrahedron, shape ``(n, 4)``. + :type simplices: :class:`~numpy.ndarray` + + :arg vertices: Voronoi vertex of each tetrahedron, shape ``(n, 3)``. + :type vertices: :class:`~numpy.ndarray` + + :arg vdw_radii: van der Waals radius of each atom. + :type vdw_radii: :class:`~numpy.ndarray` + + :arg r: probe radius in Angstrom. + :type r: float + + :arg rows: optional boolean mask over the tetrahedra. Rows outside it are + reported as ``False`` without paying for the distance computation. + :meth:`deleteSimplices3d` uses it to restrict the surface pass to the + boundary shell. + :type rows: :class:`~numpy.ndarray`, optional + + :returns: boolean array of length ``len(simplices)``, ``True`` where the + probe fits. + :rtype: :class:`~numpy.ndarray` + """ + fits = np.zeros(len(simplices), dtype=bool) + if rows is None: + rows = slice(None) + elif not rows.any(): + return fits + + atom_coords = points[simplices[rows]] # (m, 4, 3) + d_sum = np.linalg.norm( + atom_coords - vertices[rows][:, None, :], axis=2).sum(axis=1) + r_sum = (r + vdw_radii[simplices[rows]]).sum(axis=1) + fits[rows] = d_sum >= r_sum + + return fits + + def deleteSimplices3d(self, points, simplices, neighbors, vertices, + vdw_radii, r, surface): + """Delete the tetrahedra that fail the :meth:`sphereFit` probe test and + return the compacted tessellation. + + Which side of the test is deleted depends on the pass. The ``surface`` + pass erodes from the outside in: a boundary tetrahedron the probe fits + into is open to the solvent, hence exterior, and goes. Only tetrahedra on + the boundary (those with a ``-1`` neighbour) are reachable from outside, + so the test is restricted to that shell (~n^(2/3) rows) instead of being + evaluated over every tetrahedron on each erosion iteration, and the + caller iterates to convergence. The inner pass instead drops every + tetrahedron too tight for the probe, leaving the space it can actually + occupy. + + :arg points: coordinates of all atoms, shape ``(n_atoms, 3)``. + :type points: :class:`~numpy.ndarray` + + :arg simplices: atom indices of each tetrahedron, shape ``(n, 4)``. + :type simplices: :class:`~numpy.ndarray` + + :arg neighbors: index of the tetrahedron opposite each vertex, ``-1`` on + the boundary, shape ``(n, 4)``. + :type neighbors: :class:`~numpy.ndarray` + + :arg vertices: Voronoi vertex of each tetrahedron, shape ``(n, 3)``. + :type vertices: :class:`~numpy.ndarray` + + :arg vdw_radii: van der Waals radius of each atom. + :type vdw_radii: :class:`~numpy.ndarray` + + :arg r: probe radius in Angstrom. + :type r: float + + :arg surface: ``True`` for one erosion step of the surface pass, ``False`` + for the inner pass. + :type surface: bool + + :returns: the surviving ``(simplices, neighbors, vertices)``, with + neighbour indices remapped to the new numbering and deleted + neighbours set to ``-1``. + :rtype: tuple + """ + simplices = np.asarray(simplices) + neighbors = np.asarray(neighbors) + vertices = np.asarray(vertices) + + n = len(simplices) + if n == 0: + return simplices, neighbors, vertices + + if surface: + # Erode: a boundary tetrahedron the probe fits into is open to the + # solvent, hence exterior. Only the boundary shell is reachable from + # outside, so only it is tested. + boundary = (neighbors == -1).any(axis=1) + should_delete = self.sphereFit(points, simplices, vertices, + vdw_radii, r, rows=boundary) + else: + # Carve: drop every tetrahedron too tight for the probe, anywhere, + # leaving the space it can actually occupy. + should_delete = ~self.sphereFit(points, simplices, vertices, + vdw_radii, r) + + keep = ~should_delete + simp = simplices[keep] + neigh = neighbors[keep].copy() + verti = vertices[keep] + + # Remap neighbour indices from the old numbering to the compacted one in a + # single pass (deleted neighbours -> -1), replacing the previous + # O(len(deleted) x len(neigh)) decrement loop. + new_index = np.full(n, -1, dtype=neigh.dtype) + new_index[keep] = np.arange(keep.sum(), dtype=neigh.dtype) + neigh = np.where(neigh == -1, -1, new_index[neigh]) - def delete_simplices3d(self, points, simplices, neighbors, vertices, vdw_radii, r, surface): - simp, neigh, verti, deleted = [], [], [], [] - - for i, tetrahedron in enumerate(simplices): - should_delete = (-1 in neighbors[i] and self.sphere_fit(points, tetrahedron, vertices[i], vdw_radii, r)) if surface else not self.sphere_fit(points, tetrahedron, vertices[i], vdw_radii, r) - - if should_delete: - deleted.append(i) - else: - simp.append(simplices[i]) - neigh.append(neighbors[i]) - verti.append(vertices[i]) - - simp = np.array(simp) - neigh = np.array(neigh) - verti = np.array(verti) - deleted = np.array(deleted) - - mask = np.isin(neigh, deleted) - neigh[mask] = -1 - - for i in reversed(deleted): - mask = (neigh > i) & (neigh != -1) - neigh[mask] -= 1 - return simp, neigh, verti - def delete_section(self, simplices_subset, simplices, neighbors, vertices, reverse=False): - simp, neigh, verti, deleted = [], [], [], [] - - for i, tetrahedron in enumerate(simplices): - match = any((simplices_subset == tetrahedron).all(axis=1)) - if reverse: - if match: - simp.append(tetrahedron) - neigh.append(neighbors[i]) - verti.append(vertices[i]) - else: - deleted.append(i) - else: - if match: - deleted.append(i) - else: - simp.append(tetrahedron) - neigh.append(neighbors[i]) - verti.append(vertices[i]) + def calcEnclosure(self, query, centers, tree=None): + """Fraction of the directions seen from each point of ``query`` that are + blocked by an atom within :data:`ENCLOSURE_RANGE` Angstrom. + + A local, probe-independent measure of burial: a point in the open solvent + sees sky in most directions, a point inside a channel is surrounded + whatever the channel's width. + + Rays are marched outwards and a ray is dropped as soon as it is blocked, + which is what keeps this affordable: in a buried region most directions + hit protein within the first few Angstrom, and only the few that escape + are followed the whole way out. Marching costs + ``ENCLOSURE_RAYS x steps`` tree queries per point and so is all but + insensitive to how many atoms there are, whereas testing every atom in + range against every ray costs a multiple of the atom count, and with rays + this sparse nearly all of that work is wasted on atoms that lie near no + ray at all. + + Pass the real atoms here, not the balls of a homogenized diagram: burial + is a property of the protein, not of the tessellation. Atoms are all given + the same :data:`ENCLOSURE_RADIUS`, so one tree and one plain + nearest-neighbour test suffice. This is a burial heuristic and not a + surface calculation, and the alternative -- a per-atom radius, which no + nearest-neighbour query can express -- buys nothing that + ``min_enclosure`` cannot absorb. + + :data:`ENCLOSURE_RAYS` is fixed rather than exposed, because it is part of + the definition of the quantity and not an accuracy knob. Adding rays is + not a free refinement: more directions discover more of the thin escape + routes out of a channel, so the enclosure of every point drifts downwards + and a threshold calibrated at one ray count does not carry over to + another. + + :arg query: points to evaluate, ``(n, 3)``. + :arg centers: atom centres. + :arg tree: optional prebuilt :class:`~scipy.spatial.cKDTree` over + ``centers``, to avoid rebuilding it on every call. + :returns: ``n`` fractions in ``[0, 1]``.""" + from scipy.spatial import cKDTree + + query = np.asarray(query, dtype=float) + if len(query) == 0: + return np.empty(0) + if tree is None: + tree = cKDTree(np.asarray(centers, dtype=float)) + + i = np.arange(ENCLOSURE_RAYS) + 0.5 + phi = np.arccos(1 - 2 * i / ENCLOSURE_RAYS) + theta = np.pi * (1 + 5 ** 0.5) * i # Fibonacci sphere + directions = np.stack([np.sin(phi) * np.cos(theta), + np.sin(phi) * np.sin(theta), + np.cos(phi)], axis=1) + + blocked = np.zeros((len(query), ENCLOSURE_RAYS), dtype=bool) + live = np.ones_like(blocked) + for step in np.arange(ENCLOSURE_STEP, + ENCLOSURE_RANGE + ENCLOSURE_STEP, ENCLOSURE_STEP): + point, ray = np.nonzero(live) + if not len(point): + break + samples = query[point] + directions[ray] * step + hit = tree.query(samples, distance_upper_bound=ENCLOSURE_RADIUS, + workers=-1)[0] <= ENCLOSURE_RADIUS + blocked[point[hit], ray[hit]] = True + live[point[hit], ray[hit]] = False + + return blocked.mean(axis=1) + + def peelSurfaceByEnclosure(self, points, simplices, neighbors, vertices, + vdw_radii, r, atom_coords, min_enclosure, + max_depth=None): + """Erode the surface inward with a probe of radius ``r``, stopping where + the tetrahedra stop being open to the solvent. + + This removes the "moat": the shell of true exterior that lies inside the + ``r1`` surface, because an ``r1`` probe cannot enter the concavities it + bridges over. Left in place the moat joins the cavity and offers wide, + cheap routes along the outside of the protein. + + Neither obvious way of bounding the erosion works. A count of tetrahedron + layers is not mesh-invariant, since a layer is one tetrahedron thick and + tetrahedra shrink as the tessellation is refined. A depth in Angstrom is + not ``r1``-invariant, since the moat has no constant thickness: it is as + deep as ``r1 - r`` inside a concavity and vanishes on a flat face, so a + depth large enough to clear it where it is thick also marches down the + channel mouths and erodes the channels themselves. At ``r1 = 10``, a + reasonable setting for a porin or a ribosome, that leaves almost nothing. + + The rule used here is local instead. A boundary tetrahedron is stripped + only while it is *open*, that is while its enclosure is below + ``min_enclosure`` (see :meth:`calcEnclosure`). The moat is open by + construction and goes; erosion halts by itself at the first buried layer. + ``r1`` then decides only where the erosion starts, not where it stops, so + the result no longer depends on it, and ``r1`` is left doing the one job + it should: capping the mouths. + + Since enclosure is a static field, the peel is really "delete the + outside-connected component of ``{enclosure < min_enclosure}``". A + threshold above the enclosure of a channel interior (empirically about + 0.93, as a channel is itself an escape direction) therefore percolates + along the channels and erodes the cavity away entirely. ``max_depth`` is + available as a hard backstop, but the default threshold leaves a wide + margin and the failure mode is loud - no channels at all - rather than a + plausible-looking result with the real channels missing. + + Note that the probe test and the enclosure test deliberately run against + different spheres. ``points`` and ``vdw_radii`` are the balls the diagram + is built on, and the probe test has to use them or its geometry stops + agreeing with :meth:`deleteSimplices3d`. ``atom_coords`` are the real + atoms, and the enclosure test has to use those, or burial would depend on + the tessellation. Under ``diagram="homogenized"`` the two are not the same + set, as some 4700 atoms become some 33000 equal balls, which would make + the enclosure test both slower and a function of ``max_deviation``. Under + ``"simple"`` and ``"weighted"`` they coincide. + + :arg atom_coords: the real atoms, for the enclosure test. + :arg min_enclosure: fraction of directions that must be blocked for a + tetrahedron to count as interior and stop the erosion. ``<= 0`` + returns the state unchanged. + :arg max_depth: optional cap, in Angstrom, on how far the front may + advance from the initial surface. ``None`` (default) is uncapped. + :returns: ``(simplices, neighbors, vertices)``, compacted.""" + from scipy.spatial import cKDTree + + simplices = np.asarray(simplices) + neighbors = np.asarray(neighbors) + vertices = np.asarray(vertices) + + if min_enclosure <= 0 or len(simplices) == 0: + return simplices, neighbors, vertices + + boundary = (neighbors == -1).any(axis=1) + if not boundary.any(): + return simplices, neighbors, vertices + + # Fixed for the whole peel, so the cap bounds the total advance of the + # front rather than its advance per pass. + surface = cKDTree(vertices[boundary]) if max_depth is not None else None + atoms = cKDTree(atom_coords) + # Enclosure is a property of a point, not of the shrinking mesh, so a + # tetrahedron re-examined on a later pass is never re-traced. Tetrahedra + # are renumbered by the compaction below, but the four balls they are + # built on are not, so those index the cache. + traced = {} + + while True: + n = len(simplices) + if n == 0: + break + boundary = np.nonzero((neighbors == -1).any(axis=1))[0] + if not len(boundary): + break + + # The cheap tests first: does the probe fit (the same sum-based test + # as deleteSimplices3d), and are we still inside the optional cap? + ball_coords = points[simplices[boundary]] + d_sum = np.linalg.norm( + ball_coords - vertices[boundary][:, None, :], axis=2).sum(axis=1) + r_sum = (r + vdw_radii[simplices[boundary]]).sum(axis=1) + candidate = d_sum >= r_sum + if max_depth is not None: + candidate &= surface.query(vertices[boundary])[0] <= max_depth + candidates = boundary[candidate] + if not len(candidates): + break + + # Ray tracing runs only on what survived those, and only once each. + keys = [tuple(key) for key in simplices[candidates]] + fresh = [i for i, key in enumerate(keys) if key not in traced] + if fresh: + values = self.calcEnclosure(vertices[candidates[fresh]], + atom_coords, tree=atoms) + for i, value in zip(fresh, values): + traced[keys[i]] = value + enclosure = np.array([traced[key] for key in keys]) + + should_delete = np.zeros(n, dtype=bool) + should_delete[candidates[enclosure < min_enclosure]] = True + if not should_delete.any(): + break + + keep = ~should_delete + simplices = simplices[keep] + neigh = neighbors[keep].copy() + vertices = vertices[keep] + + new_index = np.full(n, -1, dtype=neigh.dtype) + new_index[keep] = np.arange(keep.sum(), dtype=neigh.dtype) + neighbors = np.where(neigh == -1, -1, new_index[neigh]) + + return simplices, neighbors, vertices + + def deleteSection(self, simplices_subset, simplices, neighbors, vertices, + reverse=False): + simplices = np.asarray(simplices) + neighbors = np.asarray(neighbors) + vertices = np.asarray(vertices) + + n = len(simplices) + if n == 0: + return simplices, neighbors, vertices + + # Which rows of `simplices` also appear in `simplices_subset` (exact, + # order-sensitive row match via hashed membership instead of the former + # O(n x len(subset)) per-row scan). + matches = _rowsIsin(simplices, np.asarray(simplices_subset)) + keep = matches if reverse else ~matches + + simp = simplices[keep] + neigh = neighbors[keep].copy() + verti = vertices[keep] + + new_index = np.full(n, -1, dtype=neigh.dtype) + new_index[keep] = np.arange(keep.sum(), dtype=neigh.dtype) + neigh = np.where(neigh == -1, -1, new_index[neigh]) - simp, neigh, verti = map(np.array, [simp, neigh, verti]) - deleted = np.array(deleted) - - mask = np.isin(neigh, deleted) - neigh[mask] = -1 - - for i in reversed(deleted): - neigh = np.where((neigh > i) & (neigh != -1), neigh - 1, neigh) - return simp, neigh, verti - def get_vdw_radii(self, atoms): + def getVdwRadii(self, atoms): vdw_radii_dict = { 'H': 1.20, 'HE': 1.40, 'LI': 1.82, 'BE': 1.53, 'B': 1.92, 'C': 1.70, 'N': 1.55, 'O': 1.52, 'F': 1.47, 'NE': 1.54, 'NA': 2.27, 'MG': 1.73, @@ -2382,45 +3409,191 @@ def get_vdw_radii(self, atoms): return np.array([vdw_radii_dict[atom] for atom in atoms]) - def surface_layer(self, shape_simplices, filtered_simplices, shape_neighbors): - surface_simplices, surface_neighbors = [], [] - interior_simplices = [] - - for i in range(len(shape_simplices)): - if -1 in shape_neighbors[i]: - surface_simplices.append(shape_simplices[i]) - surface_neighbors.append(shape_neighbors[i]) - else: - interior_simplices.append(shape_simplices[i]) - - surface_simplices = np.array(surface_simplices) - surface_neighbors = np.array(surface_neighbors) - interior_simplices = np.array(interior_simplices) - - filtered_surface_simplices = surface_simplices[ - np.any(np.all(surface_simplices[:, None] == filtered_simplices, axis=2), axis=1) - ] - filtered_surface_neighbors = surface_neighbors[ - np.any(np.all(surface_simplices[:, None] == filtered_simplices, axis=2), axis=1) - ] - + def _fibonacciSphere(self, n): + """Return ``n`` roughly evenly distributed unit vectors on a sphere using + the Fibonacci (golden spiral) lattice.""" + n = int(np.maximum(1, n)) + indices = np.arange(n) + 0.5 + phi = np.arccos(1.0 - 2.0 * indices / n) + theta = np.pi * (1.0 + 5.0 ** 0.5) * indices + x = np.sin(phi) * np.cos(theta) + y = np.sin(phi) * np.sin(theta) + z = np.cos(phi) + return np.stack([x, y, z], axis=1) + + def _shellPointCount(self, rad, rho, max_deviation): + """Number of equal balls of radius ``rho`` to place on a shell of radius + ``rad`` so that the outer envelope of their union stays within + ``max_deviation`` of ``rad + rho``. + + A ball centered at radius ``rad`` only touches the target sphere of radius + ``rad + rho`` at a single point, so the shell must be sampled densely + enough that the "valleys" between neighbouring balls do not dip more than + ``max_deviation``. Each ball covers a spherical cap of half-angle + ``alpha`` on the ``rad + rho - max_deviation`` sphere; the count is the + number of such caps needed to tile the sphere (with an overlap factor). + """ + r = rad + rho - max_deviation + cos_alpha = (rad * rad + r * r - rho * rho) / (2.0 * rad * r) + cos_alpha = float(np.clip(cos_alpha, -1.0, 1.0)) + if cos_alpha >= 1.0: + return 1 + # The exact number of caps to tile the sphere is 2 / (1 - cos_alpha); we use + # a factor of 4 (a ~2x overlap margin). This is NOT slack to be trimmed: the + # Fibonacci lattice is not an optimal packing and its coverage efficiency + # degrades as the shell (and count) grows, so the factor needed to actually + # hold the max_deviation bound increases with atom size. A single constant + # must therefore be sized for the largest atoms (e.g. metals); 4 keeps the + # measured dip below max_deviation across the whole range, whereas 3 already + # fails for anything larger than the thinnest shell. Lowering it silently + # breaks large-atom accuracy - tune max_deviation instead to change cost. + return int(np.ceil(4.0 / (1.0 - cos_alpha))) + + def homogenizeAtoms(self, coords, vdw_radii, max_deviation=0.2): + """Substitute every atom by a set of homogeneous balls whose common radius + equals the smallest van der Waals radius present in the structure. + + Each atom of radius ``R`` is replaced by a collection of overlapping balls + of radius ``rho = min(vdw_radii)`` arranged on concentric shells (plus a + central ball) so that their union approximates the original atomic sphere + to within ``max_deviation``. Because all resulting balls share the same + radius, an ordinary Voronoi / Delaunay tessellation of their centers yields + an accurate estimate of the additively weighted (power) Voronoi diagram of + the original atoms. This is the approach used by MolAxis and CAVER 3 and + avoids simply discarding the smaller (e.g. hydrogen) atoms. + + Atoms whose radius is within ``max_deviation`` of ``rho`` are kept as a + single ball, so a structure of similarly sized atoms is left essentially + unchanged while a structure containing hydrogens (small ``rho``) fills its + larger atoms with several balls. + + :arg coords: atomic coordinates, shape ``(N, 3)`` + :arg vdw_radii: per-atom van der Waals radii, shape ``(N,)`` + :arg max_deviation: maximum tolerated deviation (in Angstrom) between the + union surface of the substitute balls and the original atomic surface. + Smaller values are more accurate but generate more balls. Default 0.2. + :returns: a tuple ``(new_coords, new_radii)`` where every entry of + ``new_radii`` equals ``rho``. + """ + coords = np.asarray(coords, dtype=float) + vdw_radii = np.asarray(vdw_radii, dtype=float) + + rho = float(np.min(vdw_radii)) + tol = 1e-6 + new_points = [] + + for center, R in zip(coords, vdw_radii): + # Atoms within max_deviation of the smallest radius stay a single ball. + if R - rho <= max_deviation + tol: + new_points.append(center) + continue + + # Central ball plus concentric shells stepped by rho, with the + # outermost shell at (R - rho) so the union surface reaches R. + new_points.append(center) + shell_radii = list(np.arange(rho, R - rho, rho)) + if not shell_radii or (R - rho) - shell_radii[-1] > tol: + shell_radii.append(R - rho) + + for rad in shell_radii: + if rad <= tol: + continue + n = self._shellPointCount(rad, rho, max_deviation) + new_points.extend(center + rad * self._fibonacciSphere(n)) + + new_points = np.array(new_points) + new_radii = np.full(len(new_points), rho) + + return new_points, new_radii + + def buildSurfaceDepthOracle(self, coords, vdw_radii, r1, max_deviation, max_depth): + """Interior/exterior depth oracle for relabeling the additively-weighted + diagram's surface mouths (``diagram="weighted"``). + + The AW tessellation is not a clean simplicial complex: many interior 3-ball + faces are left unpaired and masquerade as surface boundaries, so channels + truncate to stubs. This builds a *homogenized* Voronoi diagram of the same + atoms (a clean simplicial complex), erodes it with an ``r1`` probe to separate + solvent (exterior) from protein (interior), and labels every tetrahedron by its + geodesic distance (A) below the molecular surface (0 = exterior/solvent). + :meth:`getSurfaceCavities` then keeps only AW exit tetrahedra whose Voronoi + vertex maps (via ``find_simplex``) to depth ``<= max_depth`` Angstrom. + + :returns: ``(delaunay, depth, max_depth)`` -- the homogenized + :class:`~scipy.spatial.Delaunay`, its per-tetrahedron geodesic depth in + Angstrom (points outside the hull are treated as depth 0), and the + passed-through threshold. + """ + from scipy.spatial import Delaunay + + hp, hrho = self.homogenizeAtoms(coords, vdw_radii, max_deviation) + delaunay = Delaunay(hp) + centers = self.calcCircumcenters(delaunay) + clearance = (np.linalg.norm(hp[delaunay.simplices] - centers[:, None, :], axis=2) + - hrho[delaunay.simplices]).min(axis=1) + neighbors = delaunay.neighbors + n = len(delaunay.simplices) + + # r1 surface erosion: peel boundary tetrahedra wide enough for the probe, from + # the hull inward, until nothing more can be removed. Survivors == interior. + alive = np.ones(n, dtype=bool) + while True: + dead = np.zeros(n, dtype=bool) + for k in range(4): + col = neighbors[:, k] + dead |= (col == -1) | ((col >= 0) & ~alive[col]) + peel = alive & dead & (clearance >= r1) + if not peel.any(): + break + alive[peel] = False + + # Geodesic depth (A) below the surface: shortest path from the exterior + # (solvent) tetrahedra inward along Voronoi edges. A physical distance, not a + # tetrahedron-layer count, so weighted_mouth_depth is a real Angstrom threshold + # that does not drift with the homogenization density. + degenerate = self._degenerateTetrahedra(delaunay.simplices, centers, hp) + scratch = np.full(n, -1, dtype=np.intp) + depth = self._geodesicDepth(np.arange(n), np.nonzero(~alive)[0], neighbors, + centers, degenerate, scratch) + # Enclosed pockets never reached from the exterior are deep (never a mouth). + reached = depth[np.isfinite(depth)] + depth[~np.isfinite(depth)] = (float(reached.max()) + 5.0) if reached.size \ + else float(max_depth + 1) + + return delaunay, depth, max_depth + + def surfaceLayer(self, shape_simplices, filtered_simplices, shape_neighbors): + shape_simplices = np.asarray(shape_simplices) + shape_neighbors = np.asarray(shape_neighbors) + filtered_simplices = np.asarray(filtered_simplices) + + # Split simplices into those touching the boundary (a -1 neighbour) and + # the interior ones, preserving order. + boundary = (shape_neighbors == -1).any(axis=1) + surface_simplices = shape_simplices[boundary] + surface_neighbors = shape_neighbors[boundary] + interior_simplices = shape_simplices[~boundary] + + # Row-membership tests replace the former (N, M, 4) broadcast temporaries. + surf_keep = _rowsIsin(surface_simplices, filtered_simplices) + filtered_surface_simplices = surface_simplices[surf_keep] + filtered_surface_neighbors = surface_neighbors[surf_keep] + filtered_surface_neighbors = np.unique(filtered_surface_neighbors) filtered_surface_neighbors = filtered_surface_neighbors[filtered_surface_neighbors != 0] - + filtered_interior_simplices = interior_simplices[ - np.any(np.all(interior_simplices[:, None] == filtered_simplices, axis=2), axis=1) - ] + _rowsIsin(interior_simplices, filtered_simplices)] surface_layer_neighbor_simplices = shape_simplices[filtered_surface_neighbors] - + second_layer = filtered_interior_simplices[ - np.any(np.all(filtered_interior_simplices[:, None] == surface_layer_neighbor_simplices, axis=2), axis=1) - ] + _rowsIsin(filtered_interior_simplices, surface_layer_neighbor_simplices)] return filtered_surface_simplices, second_layer - def find_groups(self, neigh, is_cavity=True): + def findGroups(self, neigh, is_cavity=True): x = neigh.shape[0] visited = np.zeros(x, dtype=bool) groups = [] @@ -2446,7 +3619,8 @@ def dfs(tetra_index): return groups - def get_surface_cavities(self, cavities, interior_simplices, second_layer, state, points, vdw_radii, sparsity): + def getSurfaceCavities(self, cavities, interior_simplices, second_layer, + state, points, vdw_radii, sparsity, mouth_oracle=None): surface_cavities = [] for cavity in cavities: @@ -2454,174 +3628,806 @@ def get_surface_cavities(self, cavities, interior_simplices, second_layer, state second_layer_mask = np.isin(interior_simplices[tetrahedra], second_layer).all(axis=1) if np.any(second_layer_mask): - cavity.make_surface() exit_tetrahedra = tetrahedra[second_layer_mask] - end_tetrahedra = self.get_end_tetrahedra(exit_tetrahedra, state.verti, points, vdw_radii, state.simp, sparsity) - cavity.set_exit_tetrahedra(exit_tetrahedra, end_tetrahedra) + if mouth_oracle is not None: + # diagram="weighted": drop the false (buried) mouths the leaky AW + # diagram produces. Keep an exit tetrahedron only if its Voronoi + # vertex lies within max_depth Angstrom (geodesic) of the true + # molecular surface, per a homogenized interior/exterior oracle. + delaunay, depth, max_depth = mouth_oracle + located = delaunay.find_simplex(state.verti[exit_tetrahedra]) + surface_depth = np.where(located >= 0, depth[located.clip(0)], 0.0) + exit_tetrahedra = exit_tetrahedra[surface_depth <= max_depth] + if len(exit_tetrahedra) == 0: + continue + cavity.makeSurface() + end_tetrahedra = self.getEndTetrahedra(exit_tetrahedra, state.verti, points, vdw_radii, state.simp, sparsity) + cavity.setExitTetrahedra(exit_tetrahedra, end_tetrahedra) surface_cavities.append(cavity) return surface_cavities - def merge_cavities(self, cavities, simplices): + def mergeCavities(self, cavities, simplices): + if not cavities: + # No cavities survived filtering (e.g. restrict_channels_to_start_point + # selected a single cavity shallower than min_depth). Return an empty + # (0, 4) slice so the pipeline yields zero channels instead of crashing + # in np.concatenate on an empty list. + return simplices[np.empty(0, dtype=np.intp)] merged_tetrahedra = np.concatenate([cavity.tetrahedra for cavity in cavities]) return simplices[merged_tetrahedra] - def find_deepest_tetrahedra(self, cavities, neighbors): - from collections import deque - + def _degenerateTetrahedra(self, simplices, vertices, points): + # Scale-invariant flatness flag for the geodesic depth graph. A near-flat + # tetrahedron has a runaway circumcenter, so its incident Voronoi edges can be + # astronomically long (measured up to ~1e14 A) and would corrupt a shortest-path + # depth. Flatness is the ratio of circumradius R to the tetrahedron's own atom + # span L: well-shaped cells sit at R/L < ~4 whatever their absolute size - the + # large fat cells that span a wide pore under a big probe included - while flat + # slivers diverge to R/L -> infinity. Using the dimensionless R/L, never an + # absolute length, keeps this correct for porins, ribosome tunnels and + # large-probe (e.g. r1=20) runs alike. Edges touching a flagged tetrahedron are + # dropped from the graph in _geodesicDepth. + apex = points[simplices] # (n, 4, 3) + R = np.linalg.norm(apex[:, 0] - vertices, axis=1) # circumradius + L = np.zeros(len(simplices)) + for a in range(4): + for b in range(a + 1, 4): + L = np.maximum(L, np.linalg.norm(apex[:, a] - apex[:, b], axis=1)) + # Well-shaped cells measure R/L up to ~3.7 (p99.9); flat slivers reach ~1e14. + # Flag above 5 - the ~14-order gap makes the exact cut irrelevant across [5, 100]. + return R / np.maximum(L, 1e-9) > 5.0 + + def _geodesicDepth(self, tetrahedra, sources, neighbors, vertices, degenerate, scratch): + # Shortest-path distance (A) from the source set to every tetrahedron in + # `tetrahedra`, along Voronoi edges weighted by circumcenter-to-circumcenter + # distance (a multi-source Dijkstra over the induced subgraph). Edges touching a + # degenerate (runaway-circumcenter) tetrahedron are dropped. `scratch` is a + # reusable global->local index buffer (-1 outside `tetrahedra`), reset before + # return so the caller can pass it again. Returns the distance aligned to + # `tetrahedra`, np.inf where a tetrahedron is unreachable from every source. + from scipy.sparse import csr_matrix + from scipy.sparse.csgraph import dijkstra + + T = np.asarray(tetrahedra, dtype=np.intp) + m = len(T) + if m == 0: + return np.empty(0) + scratch[T] = np.arange(m) + nb = neighbors[T] # (m, deg) global + safe_nb = np.where(nb < 0, 0, nb) + local = np.where(nb < 0, -1, scratch[safe_nb]) + keep = local >= 0 + keep &= ~degenerate[T][:, None] # drop from a flat tetra + keep &= ~np.where(nb < 0, True, degenerate[safe_nb]) # drop into a flat tetra + keep = keep.ravel() + row = np.repeat(np.arange(m), nb.shape[1])[keep] + col = local.ravel()[keep] + w = np.linalg.norm(vertices[T[row]] - vertices[T[col]], axis=1) + graph = csr_matrix((w, (row, col)), shape=(m, m)) + src = scratch[np.asarray(sources, dtype=np.intp)] + src = src[src >= 0] + scratch[T] = -1 # reset for reuse + if len(src) == 0: + return np.full(m, np.inf) + return dijkstra(graph, indices=src, min_only=True) + + def findDeepestTetrahedra(self, cavities, neighbors, vertices, points, simplices): + # Cavity depth = the geodesic distance (A) from the cavity's openings (its exit + # tetrahedra) to its farthest point, as a shortest path along Voronoi edges. This + # is a physical length independent of tetrahedron size, so min_depth is mesh + # invariant; the old +1-per-tetrahedron layer count grew as the mesh refined. + degenerate = self._degenerateTetrahedra(simplices, vertices, points) + scratch = np.full(neighbors.shape[0], -1, dtype=np.intp) for cavity in cavities: - exit_tetrahedra = cavity.exit_tetrahedra - visited = np.zeros(neighbors.shape[0], dtype=bool) - visited[exit_tetrahedra] = True - queue = deque([(tetra, 0) for tetra in exit_tetrahedra]) - max_depth = -1 - deepest_tetrahedron = None - tetrahedra_depths = {} - - while queue: - current, depth = queue.popleft() - tetrahedra_depths[current] = depth - - if depth > max_depth: - max_depth = depth - deepest_tetrahedron = current - - for neighbor in neighbors[current]: - if neighbor != -1 and not visited[neighbor] and neighbor in cavity.tetrahedra: - visited[neighbor] = True - queue.append((neighbor, depth + 1)) - - cavity.set_starting_tetrahedron(np.array([deepest_tetrahedron])) - cavity.set_depth(max_depth) - cavity.tetrahedra_depths = tetrahedra_depths - - def dijkstra(self, cavity, simplices, neighbors, vertices, points, vdw_radii): - import heapq - - def calculate_weight(current_tetra, neighbor_tetra): - current_vertex = vertices[current_tetra] - neighbor_vertex = vertices[neighbor_tetra] - l = np.linalg.norm(current_vertex - neighbor_vertex) - - d = np.inf - for atom, radius in zip(points[simplices[neighbor_tetra]], vdw_radii[simplices[neighbor_tetra]]): - dist = np.linalg.norm(neighbor_vertex - atom) - radius - if dist < d: - d = dist + tetra = np.asarray(cavity.tetrahedra, dtype=np.intp) + dist = self._geodesicDepth(tetra, cavity.exit_tetrahedra, neighbors, + vertices, degenerate, scratch) + finite = np.isfinite(dist) + if not finite.any(): + # No exit reached any tetrahedron (degenerate cavity); keep it minimal. + cavity.setStartingTetrahedron(np.array([int(tetra[0])])) + cavity.setDepth(0.0) + cavity.tetrahedra_depths = {int(tetra[0]): 0.0} + continue + deepest = int(np.argmax(np.where(finite, dist, -np.inf))) + cavity.setStartingTetrahedron(np.array([int(tetra[deepest])])) + cavity.setDepth(float(dist[deepest])) + cavity.tetrahedra_depths = {int(tetra[k]): float(dist[k]) + for k in np.nonzero(finite)[0]} - b = 1e-3 - return l / (d**2 + b) - - def dijkstra_algorithm(start, goal, tetrahedra_set): - pq = [(0, start)] - distances = {start: 0} - previous = {start: None} + def calcCircumcenters(self, dela): + # per-simplex circumcenters recovered analytically from the Delaunay + # paraboloid lifting, avoiding a second Qhull pass. Identical to scipy + # Voronoi vertices in general position. + eq = dela.equations + scale = dela.paraboloid_scale + centers = -eq[:, :-2] / (2 * scale * eq[:, -2][:, None]) + return centers + + def _edgeBottleneck(self, ci, cj, shared_atoms, points, vdw_radii): + # Minimum clearance along the Voronoi edge - the segment between the two + # circumcenters ci, cj, dual to the Delaunay face the two tetrahedra + # share - measured against that face's atoms. This is the edge bottleneck + # radius: the circumcenters are local clearance maxima, so the tightest + # point of the segment (the gate) generally lies between them and is + # narrower than either endpoint. + # + # min over t in [0, 1] and over the shared atoms of |p(t) - a| - vdw(a), + # with p(t) = ci + t (cj - ci). Because the min over the (t, atom) + # product equals the min of the per-atom minima, each atom reduces to an + # independent clamped point-to-segment distance - a closed form, no + # sampling. The foot clamps to an endpoint when it falls outside the + # segment, so the gate value is always <= both vertex clearances. Exact + # for the straight edges of the homogenized/simple diagrams; for the + # weighted (Apollonius) diagram the true edge is a slight arc and the + # chord is a local approximation. + a = points[shared_atoms] + r = vdw_radii[shared_atoms] + u = cj - ci + uu = float(u @ u) + if uu <= 1e-12: + # Twin tetrahedra: the circumcenters coincide, the edge is a point, + # so return the shared vertex clearance directly. This guard is + # load-bearing, not defensive boilerplate - do not remove it. Two + # near-cospherical tetrahedra (the "T5" degeneracy CAVER's paper + # jitters against) produce coincident circumcenters; measured on real + # structures these are absent at the default max_deviation but do + # appear - a handful per structure - as the diagram is refined + # (max_deviation -> 0.02). scipy/Qhull's cospherical facet merging + # keeps the geometry finite (no NaN/inf circumcenters) so we do not + # need CAVER's jitter/random-rotation precautions, but the coincident + # circumcenters would still divide by ~0 here. Just above the + # threshold the clip below handles it: for a tiny but non-zero edge + # the foot clamps to an endpoint and the gate degrades continuously + # to the endpoint clearance. + return float(np.min(np.linalg.norm(a - ci, axis=1) - r)) + t = np.clip((a - ci) @ u / uu, 0.0, 1.0) + p = ci + t[:, None] * u + return float(np.min(np.linalg.norm(p - a, axis=1) - r)) + + def _edgeBottleneckBatch(self, ci, cj, shared_atoms, points, vdw_radii): + # Vectorized _edgeBottleneck over F edges at once. ``ci``, ``cj`` are + # (F, 3) circumcenters (``ci`` the lower-index endpoint, so the reverse + # edge yields a bitwise-identical gate) and ``shared_atoms`` is (F, 3) + # atom indices of the shared Delaunay face. Same closed-form clamped + # point-to-segment as the scalar version, with the twin-tetrahedron + # (coincident circumcenter) guard applied row-wise. See _edgeBottleneck. + # Returns ``(gate, tstar)``: the min clearance and the parameter t in + # [0, 1] where it is attained (the binding atom's clamped foot), the + # latter used by _edgeCostIntegralBatch to force a quadrature node on + # the pinch. + a = points[shared_atoms] # (F, 3, 3) + r = vdw_radii[shared_atoms] # (F, 3) + u = cj - ci # (F, 3) + uu = np.einsum('fj,fj->f', u, u) # (F,) + twin = uu <= 1e-12 + uu_safe = np.where(twin, 1.0, uu) + diff = a - ci[:, None, :] # (F, 3, 3) + t = np.clip(np.einsum('faj,fj->fa', diff, u) / uu_safe[:, None], + 0.0, 1.0) # (F, 3) + p = ci[:, None, :] + t[:, :, None] * u[:, None, :] + clr = np.linalg.norm(p - a, axis=2) - r # (F, 3) + amin = clr.argmin(axis=1) + rows = np.arange(len(u)) + gate = clr[rows, amin] + tstar = t[rows, amin] + if twin.any(): + # coincident circumcenters: the edge is a point, so the gate is the + # shared vertex clearance measured at ci (== cj) and tstar is moot. + twin_gate = (np.linalg.norm(diff, axis=2) - r).min(axis=1) + gate = np.where(twin, twin_gate, gate) + tstar = np.where(twin, 0.0, tstar) + return gate, tstar + + def _edgeCostIntegralBatch(self, ci, cj, tstar, shared_atoms, points, + vdw_radii, z=2.0, delta=0.3, r_floor=1e-2): + # Price the edge by the integral of its clearance profile r(t)^-z, not by + # the whole length charged at its single narrowest point (the l/gate^2 MOLE + # cost). This profile-integral idea follows CAVER (TCBB'15, Eq. 1), but this + # is NOT a CAVER reimplementation - the quadrature deliberately differs (see + # below). r(t) = min over the 3 shared-face balls of |ci + t (cj-ci) - a| - + # vdw is the exact clearance profile of the straight (homogenized/simple) + # edge. The integral is additive under subdivision, so unlike l/gate^2 the + # cost is mesh-invariant; the vertex-only formula's error grows with edge + # length (measured ~19% -> ~1500% p95 across length bins) and drifts the + # routing as max_deviation coarsens. + # + # The quadrature differs from CAVER's: where CAVER samples a plain uniform + # grid (and takes its grid-minimum as the bottleneck), we take the uniform + # grid linspace(0, 1, K) AND force a node at the exact analytic gate t*, so + # the narrowest point is never missed - and the reported edge bottleneck is + # that exact gate (see _edgeBottleneckBatch), not a grid sample. K = + # ceil(L/delta) is a fixed ARCLENGTH step in Angstrom, so the sample count + # scales with physical length (what makes it mesh-invariant). A short edge + # (L <= delta) collapses to {0, t*, 1}, the three clearances already in + # hand. r(t) is floored at r_floor so the integrand cannot diverge on a + # sub-r2 edge that dips through an atom (the integral's analog of the + # l/(d^2+b) regularizer); traversable edges have r(t) >= r2 >> r_floor and + # are untouched. Exact for straight edges; a chord approximation for the + # weighted (Apollonius) diagram, whose edges are arcs. + chunk = 20000 + u = cj - ci + L = np.linalg.norm(u, axis=1) + cost = np.zeros(len(ci)) + alive = np.nonzero(L > 1e-6)[0] # twin/zero-length -> 0 + K = np.maximum(2, np.ceil(L / delta).astype(int) + 1) + for k in np.unique(K[alive]): + bucket = alive[K[alive] == k] + grid = np.linspace(0.0, 1.0, k) + for s in range(0, len(bucket), chunk): + ii = bucket[s:s + chunk] + nodes = np.sort(np.concatenate( + [np.broadcast_to(grid, (len(ii), k)), tstar[ii, None]], + axis=1), axis=1) # (n, k+1) + p = ci[ii][:, None, :] + nodes[:, :, None] * u[ii][:, None, :] + a = points[shared_atoms[ii]] # (n, 3, 3) + rr = np.linalg.norm(p[:, :, None, :] - a[:, None, :, :], axis=3) \ + - vdw_radii[shared_atoms[ii]][:, None, :] + rr = np.maximum(rr.min(axis=2), r_floor) # (n, k+1) profile + cost[ii] = np.trapz(rr ** (-z), nodes, axis=1) * L[ii] + return cost + + def buildSparseGraph(self, simplices, neighbors, vertices, points, vdw_radii): + # One weighted CSR adjacency matrix for the whole cleared state, built + # from array ops over the (N, deg) neighbour table - no Python loop. + # Edge (tetra -> neigh) weight is l / (d**2 + b), where l is the + # vertex-to-vertex distance and d is the gate clearance on the shared + # Delaunay face (min clearance along the connecting Voronoi edge). The + # gate is the width the cost should see: the clearance between the two + # circumcenters, not the entered node's own vertex clearance, which is a + # local maximum and lets the search prefer a route that is actually + # narrower at a face it never measures. Because the gate is symmetric, + # the cost no longer depends on traversal direction the way the + # entered-node vertex clearance did. + # + # Two edge-cost modes (self.edge_cost): 'bottleneck' is the l/(d**2 + b) + # above (d = gate); 'integral' replaces it on face edges with a clearance- + # profile integral (_edgeCostIntegralBatch), which is mesh-invariant. + # Either way the gate is still cached as the reported edge bottleneck. + from scipy.sparse import csr_matrix + + simplices = np.asarray(simplices) + neighbors = np.asarray(neighbors) + N, deg = neighbors.shape + + tetra_points = points[simplices] + distances = np.linalg.norm(tetra_points - vertices[:, None, :], axis=2) + bottleneck = np.min(distances - vdw_radii[simplices], axis=1) + # Cache the per-simplex vertex clearance: the channel geometry otherwise + # recomputes this identical min over each path's tetrahedra. + self._vertex_clearance = bottleneck + + # Directed edge list straight off the neighbour table. + rows = np.repeat(np.arange(N), deg) + cols = neighbors.ravel() + keep = cols != -1 + rows = rows[keep] + cols = cols[keep] + + # Drive the gate geometry from the lower-index endpoint so an edge's two + # directed copies see identical inputs and get a bitwise identical, + # direction-symmetric gate - i.e. compute each undirected edge once. + lo = np.minimum(rows, cols) + hi = np.maximum(rows, cols) + + # Shared 3 atoms of each Delaunay face, convention-agnostic (set + # intersection, not scipy's opposite-vertex rule, so it holds for the + # weighted diagram too): ``present`` marks which of the lower simplex's + # four atoms also appear in the higher one; a face-adjacent edge has 3. + slo = simplices[lo] + shi = simplices[hi] + present = (slo[:, :, None] == shi[:, None, :]).any(axis=2) # (M, 4) + face = present.sum(axis=1) == 3 + + # Rare non-face links fall back to the entered node's vertex clearance; + # face links (essentially all of them) overwrite it with the gate. + d = bottleneck[cols].astype(float, copy=True) + fi = tstar = shared = None + if face.any(): + fi = np.nonzero(face)[0] + shared = slo[fi][present[fi]].reshape(-1, 3) + d[fi], tstar = self._edgeBottleneckBatch(vertices[lo[fi]], + vertices[hi[fi]], shared, + points, vdw_radii) + + l = np.linalg.norm(vertices[rows] - vertices[cols], axis=1) + b = 1e-3 + weight = l / (d * d + b) # 'bottleneck' cost; non-face fallback + if self.edge_cost == 'integral' and fi is not None: + # 'integral' cost: replace the l/(d^2+b) of every face edge (the + # bottleneck-only MOLE cost) with the clearance-profile integral. The + # l/(d^2+b) fallback stays on the rare non-face links. + # + # No R/L flatness guard is needed here (unlike buildSurfaceDepthOracle, + # which runs on the full diagram): this graph is the *cleared interior* + # state, and a runaway circumcenter means a huge clearance, so the r1 + # erosion has already stripped every flat boundary tetra - measured 0 + # degenerate tetra and max edge ~4 A in the cleared graph. That matters + # because l/(d^2+b) *over*-prices a runaway edge (huge l) so the search + # avoids it for free, whereas the integral would *under*-price it (r(t) + # is huge over most of a runaway edge, so r(t)^-z ~ 0 there). + weight[fi] = self._edgeCostIntegralBatch( + vertices[lo[fi]], vertices[hi[fi]], tstar, shared, points, vdw_radii) + + # Per-edge gate cache (unordered key) read by _pathGates: each face edge + # stored once, keyed (lo, hi). The reported bottleneck reads the same map. + undirected = face & (rows < cols) + self._edge_bottleneck = { + (int(i), int(j)): float(v) + for i, j, v in zip(rows[undirected], cols[undirected], + d[undirected]) + } - while pq: - current_distance, current_tetra = heapq.heappop(pq) - - if current_tetra == goal: - path = [] - while current_tetra is not None: - path.append(current_tetra) - current_tetra = previous[current_tetra] - return path[::-1] - - if current_distance > distances[current_tetra]: + return csr_matrix((weight, (rows, cols)), shape=(N, N)) + + def dijkstra(self, cavity, graph, simplices, neighbors, vertices, points, + vdw_radii, truncate_at_surface=True, similarity=0.8): + # a single multi-target Dijkstra from the seed over the cavity subgraph, + # then every exit path reconstructed from the predecessor tree - + # instead of one heap search per (seed, exit) pair. + # Channel geometry still goes through the current + # process_channel/Channel (Simpson-based volume). + from scipy.sparse.csgraph import dijkstra + from collections import defaultdict + + cavity_tetra = np.asarray(cavity.tetrahedra) + if len(cavity_tetra) == 0: + return + global_to_local = {tetra: i for i, tetra in enumerate(cavity_tetra)} + cavity_graph = graph[np.ix_(cavity_tetra, cavity_tetra)] + + # A tunnel ends at the surface, but the Dijkstra cost has no such term: + # it rewards width, and the widest places are the surface grooves. Left + # free, the cheapest path to a far exit leaves the pocket at one mouth, + # runs along the outside and re-enters at another - which is not a + # tunnel. So a mouth must not *conduct*, only *absorb*: we drop the + # outgoing edges of every mouth before the search, making "a channel + # ends at its first surface contact" a hard constraint of the search + # rather than a cut applied afterwards to the winning path. + # The ordering is the whole point. Truncating after selection cuts a + # path that was itself chosen *because* it ran along the surface, while + # a genuine narrow interior corridor to the same mouth loses the argmin + # to that groove and is never enumerated at all - it vanishes from the + # output even though it is open. Mesh refinement makes the groove + # cheaper, so interior tunnels drop out one by one as max_deviation + # shrinks; forbidding transit removes that dependence entirely. + # A mouth is a surface (exit) tetrahedron a probe of the traversal + # radius r2 can leave through. Note the gate is r2, not bottleneck: + # bottleneck is a reporting filter, and letting it decide which mouths + # absorb would let it silently re-open narrow mouths as transit nodes, + # i.e. change the routes rather than filter them. For the homogenized + # and weighted diagrams every surviving tetrahedron already clears r2 by + # construction (equal radii + equidistant circumcenter collapse the + # sum-based test in deleteSimplices3d to the per-atom clearance), so the + # test is a no-op there; it earns its keep for diagram="simple", where + # unequal radii break that identity. + # Local indices of the tetrahedra a channel is allowed to end at. + terminals_local = [global_to_local[int(t)] + for t in np.asarray(cavity.end_tetrahedra) + if int(t) in global_to_local] + if truncate_at_surface: + exit_tetra = np.asarray(getattr(cavity, 'exit_tetrahedra', + np.empty(0, dtype=np.intp))) + if len(exit_tetra): + verts = vertices[exit_tetra] + atom_pos = points[simplices[exit_tetra]] + atom_rad = vdw_radii[simplices[exit_tetra]] + clearance = (np.linalg.norm(atom_pos - verts[:, None, :], + axis=2) - atom_rad).min(axis=1) + seeds = set(int(s) for s in cavity.starting_tetrahedron) + + # Only the mouths themselves absorb. A tetrahedron that merely lies + # inside a mouth's inscribed ball must NOT be absorbed: the ball's + # radius is the clearance (up to ~2 A) and it reaches inward as + # well as outward, so absorbing on it eats the corridors that + # approach the surface and truncates real tunnels before they + # arrive - measured to delete both known side tunnels at + # max_deviation=0.1 while keeping them at 0.02, i.e. exactly the + # silent, mesh-dependent tunnel loss this design exists to prevent. + # A path can consequently still slip *past* a mouth through a twin + # tetrahedron - a neighbour sharing almost the same circumcenter, + # not itself in the second layer and so still conducting - and + # surface again somewhere else. That leak is real but narrow (the + # twins sit 0.1-0.7 A from a mouth, in the surface shell at depth + # 1-4), and such a path always passes through the exit sphere of a + # channel that is already reported. It is therefore handled in + # _addDedupedChannels, which cuts a path at the first reported exit + # sphere it enters - the point where it truly leaves the protein - + # rather than walling the graph off against every mouth. + absorbing = [global_to_local[int(t)] + for t, c in zip(exit_tetra, clearance) + if c >= self.r2 and int(t) in global_to_local + and int(t) not in seeds] + if absorbing: + # Zero the mouths' rows: edges *into* a mouth survive (a + # channel may end there), edges *out of* it are gone. + cavity_graph = cavity_graph.tolil() + for i in absorbing: + cavity_graph.rows[i] = [] + cavity_graph.data[i] = [] + cavity_graph = cavity_graph.tocsr() + # Every mouth is a terminus; the dedup decides which of them are + # one opening. See the comment at the target loop below. + terminals_local = absorbing + + candidates = [] + + for start_global in cavity.starting_tetrahedron: + if start_global not in global_to_local: + continue + start_local = global_to_local[start_global] + # directed=True: edge (u -> v) keeps weight l / (d_v**2 + b), i.e. + # clearance of the node being *entered* - exactly the current heap + # Dijkstra's cost model. (directed=False would symmetrize each edge + # to l / (max(d_u, d_v)**2 + b) and pick slightly different paths) + distances, predecessors = dijkstra( + cavity_graph, directed=True, indices=start_local, + return_predecessors=True) + parent_to_children = defaultdict(list) + + for node, parent in enumerate(predecessors): + if parent >= 0: + parent_to_children[parent].append(node) + + paths = {} + stack = [(start_local, [start_local])] + while stack: + node, path = stack.pop() + paths[node] = path + for child in parent_to_children.get(node, []): + stack.append((child, path + [child])) + + # A channel ends where it first touches the surface, i.e. at whichever + # mouth absorbed it - so emit a candidate for every *reachable* mouth, + # not for a pre-sampled subset of them. Sampling the targets before the + # search (the old `end_tetrahedra`, thinned by `sparsity`) can pick a + # target that sits behind another mouth: the path is absorbed at that + # nearer mouth and can go no further, the sampled target is never + # reached, and because only sampled targets emit channels the tunnel is + # reported nowhere at all. Which mouth happens to shadow which target is + # a tessellation accident, so real tunnels vanished at some meshes and + # not others. Every mouth is a legitimate terminus, so let every + # reachable one produce a candidate and leave exit identity to the + # dedup, where `sparsity` merges the mouths that share one opening. + for exit_local in terminals_local: + if exit_local == start_local: continue - - for neighbor in neighbors[current_tetra]: - if neighbor in tetrahedra_set: - weight = calculate_weight(current_tetra, neighbor) - distance = current_distance + weight - if distance < distances.get(neighbor, float('inf')): - distances[neighbor] = distance - previous[neighbor] = current_tetra - heapq.heappush(pq, (distance, neighbor)) - - return None - - tetrahedra_set = set(cavity.tetrahedra) - for exit_tetrahedron in cavity.end_tetrahedra: - for starting_tetrahedron in cavity.starting_tetrahedron: - if exit_tetrahedron != starting_tetrahedron: - path = dijkstra_algorithm(starting_tetrahedron, exit_tetrahedron, tetrahedra_set) - if path: - path_tetrahedra = np.array(path) - channel = Channel(path_tetrahedra, *self.process_channel(path_tetrahedra, vertices, points, vdw_radii, simplices)) - cavity.add_channel(channel) - - def calculate_max_radius(self, vertice, points, vdw_radii, simp): + if np.isinf(distances[exit_local]): + continue + + path_local = paths.get(exit_local) + if path_local is None: + continue + + path_global = cavity_tetra[path_local] + channel = Channel(path_global, *self.processChannel( + path_global, vertices, points, vdw_radii, simplices), + cost=float(distances[path_local[-1]])) + # the Dijkstra cost at every node, so that a path cut short at a + # reported exit can be re-costed at the node it was cut at + node_costs = np.asarray(distances)[np.asarray(path_local)] + candidates.append((channel, node_costs)) + + if truncate_at_surface: + self._addDedupedChannels(cavity, candidates, similarity, vertices, + points, vdw_radii, simplices) + else: + for channel, _costs in candidates: + cavity.addChannel(channel) + + def _addDedupedChannels(self, cavity, candidates, similarity, vertices, + points, vdw_radii, simplices): + # Cheapest first, and each candidate is judged only against the channels + # already kept - so the kept channel is always the cheapest of its group + # and the result does not depend on the order candidates arrive in. + # + # Step 1, CUT. A reported channel's exit sphere (centred on its exit + # vertex, radius the clearance there) is the volume of that opening. If a + # candidate's route passes through it, the candidate has left the protein + # at that opening: whatever it does afterwards is a hop across the outside, + # not part of a tunnel. So cut it there. This is what stops a path from + # slipping past a mouth through a twin tetrahedron and claiming some far + # exit on the other side. Note the cut is made only against the handful of + # *already reported* exits, never against all mouths - truncating against + # every mouth is what used to demolish real tunnels. + # + # Step 2, COMPARE. Two channels are the same tunnel when they leave by the + # same opening AND take the same corridor to get there. Both halves are + # needed: a corridor that forks near the surface and exits twice through + # one opening is one tunnel counted twice (merge), but two genuinely + # different corridors that happen to surface at the same opening are two + # tunnels (keep), and one corridor reaching two separate openings is also + # two tunnels (keep). A candidate that was cut in step 1 is compared on its + # cut route, which is the only part of it that is really a tunnel. + # + # Route identity is measured GEOMETRICALLY - how much of one centerline + # runs alongside the other - not as a shared prefix of tetrahedra. A prefix + # is the wrong instrument twice over: it is blind to rejoining (two routes + # that split near the seed and then run together to the same exit share + # almost no prefix, yet are plainly one tunnel) and it is fooled by + # containment (a short route that is a prefix of a long one scores ~1.0 and + # deletes the long one, which is the channel carrying the distinctive + # route). Comparing the curves is immune to both, and to the tessellation: + # a tetrahedron count is not mesh-invariant, so the same physical fork + # scores differently at different max_deviation. + prepared = sorted(candidates, key=lambda c: c[0].cost) + if not prepared: + return + + kept = [] # (channel, pts, exit_xyz, opening_radius) + for channel, node_costs in prepared: + tetra = np.asarray(channel.tetrahedra) + pts = vertices[tetra] + + # step 1: cut at the first reported opening this route enters + cut, cutter = None, None + for i in range(1, len(pts)): + for kxyz, kr in ((k[2], k[3]) for k in kept): + if np.linalg.norm(pts[i] - kxyz) < kr: + cut = i + break + if cut is not None: + cutter = next(k for k in kept + if np.linalg.norm(pts[cut] - k[2]) < k[3]) + break + if cut is not None: + tetra = tetra[:cut + 1] + pts = pts[:cut + 1] + channel = Channel(tetra, *self.processChannel( + tetra, vertices, points, vdw_radii, simplices), + cost=float(node_costs[cut])) + + # step 2: same opening AND same corridor -> the same tunnel. + # The corridor is compared OUTSIDE the shared opening. Inside it the + # routes are already through the mouth and merely fanning out across + # it, and that fan is not evidence of a different corridor: the + # Voronoi network splays where a tunnel widens into its opening, so + # sibling paths peel off in the last few Angstrom and end on + # neighbouring exit tetrahedra. Counting that splay as divergence is + # what used to report one tunnel as a bundle of near-copies. + duplicate = False + for _kc, kpts, kxyz, kr in kept: + if np.linalg.norm(pts[-1] - kxyz) >= kr: + continue # a different opening + if self._routeCoverage(pts, kpts, center=kxyz, + radius=kr) >= similarity: + duplicate = True + break + if not duplicate: + if cut is not None: + # A cut channel stops inside an opening that is already + # reported, so it INHERITS that opening rather than declaring + # its own. Its last tetrahedron is an interior one that merely + # happens to lie in the exit volume, and its inscribed sphere + # is not a mouth - promoting it to a cutting surface would let + # an interior sphere start truncating other candidates. (It + # survives to here only when it reached that opening by a + # genuinely different corridor, which is a distinct tunnel and + # must be kept. Note its cost, taken at the cut node, is + # necessarily below that of the channel that cut it, since the + # cut lies upstream of that channel's mouth - so cost orders + # the output but does not mean the cut channel is "better".) + kept.append((channel, pts, cutter[2], cutter[3])) + else: + # One radius stands for this opening everywhere: it cuts routes + # that pass through it, it decides which channels share it, and + # it is the region discounted when their corridors are compared. + # The clearance at the exit vertex measures the mouth, but on a + # coarse tessellation it is erratic and can collapse to almost + # nothing, fragmenting one physical mouth into several; the + # sparsity floor keeps it mesh-independent. + kept.append((channel, pts, pts[-1], + max(self.calculateMaxRadius( + pts[-1], points, vdw_radii, + simplices[tetra[-1]]), self.sparsity))) + for channel, _pts, _xyz, _r in kept: + cavity.addChannel(channel) + + def _routeCoverage(self, a, b, tol=None, center=None, radius=0.0): + """Fraction of the SHORTER centerline's arc length that runs within ``tol`` + Angstrom of the longer one. + + Answers "does the longer channel follow the shorter one's corridor?". + ``1.0`` means the shorter route lies wholly inside the longer one's + corridor, so they took the same way out - the longer one simply carried on + past the point where the shorter one surfaced. That continuation is *not* + counted as a difference, which is the point: two channels leaving through + one opening are one tunnel even if one of them runs on and exits a few + Angstrom further along. It is safe to ignore the continuation only because + this is gated on the two channels sharing an opening; without that gate, + scoring against the shorter route would delete long channels that head off + to a quite different exit. + + ``center`` and ``radius`` describe that shared opening, and the part of + either route lying inside it is discarded before the comparison. A tunnel + splays as it widens into its mouth, so sibling paths peel apart over the + last few Angstrom and land on neighbouring exit tetrahedra; that fan says + nothing about which corridor they took, and counting it makes one tunnel + look like several. + + Note this deliberately says nothing about *where* the routes differ, or how + sharply the uncovered part turns away - only how much of the shorter route + is shared. Where two corridors genuinely part company, they do so for a + large fraction of the route, and the score falls.""" + from scipy.spatial import cKDTree + + if tol is None: + tol = self.route_tolerance + if center is not None and radius > 0: + a = a[np.linalg.norm(a - center, axis=1) > radius] + b = b[np.linalg.norm(b - center, axis=1) > radius] + if len(a) < 2 or len(b) < 2: + # Nothing survives outside the opening, so all either route ever + # did was cross the mouth: there is no corridor to tell apart. + return 1.0 + if len(a) < 2 or len(b) < 2: + return 0.0 + + def arclen(p): + return float(np.linalg.norm(np.diff(p, axis=0), axis=1).sum()) + + long_p, short_p = (a, b) if arclen(a) >= arclen(b) else (b, a) + steps = np.linalg.norm(np.diff(short_p, axis=0), axis=1) + total = steps.sum() + if total <= 0: + return 0.0 + + # each node carries half of each adjacent segment, so its weight is the + # arc length it stands for + weight = np.zeros(len(short_p)) + weight[:-1] += steps / 2.0 + weight[1:] += steps / 2.0 + + near = cKDTree(long_p).query(short_p)[0] <= tol + return float(weight[near].sum() / total) + + def calculateMaxRadius(self, vertice, points, vdw_radii, simp): atom_positions = points[simp] radii = vdw_radii[simp] distances = np.linalg.norm(atom_positions - vertice, axis=1) - radii return np.min(distances) - def calculate_radius_spline(self, tetrahedra, voronoi_vertices, points, vdw_radii, simp): - vertices = voronoi_vertices[tetrahedra] - radii = np.array([self.calculate_max_radius(v, points, vdw_radii, s) for v, s in zip(vertices, simp[tetrahedra])]) - return radii, np.min(radii) - - def process_channel(self, tetrahedra, voronoi_vertices, points, vdw_radii, simp): + def _pathGates(self, tetrahedra, voronoi_vertices, points, vdw_radii, + simp, vertex_radii): + # Per-edge gate clearance along the path: the minimum clearance on each + # shared Delaunay face between consecutive circumcenters (the edge + # bottleneck radius), read from the cache and recomputed for any edge the + # map lacks. Each gate is <= the clearance at both its endpoints, so the + # path minimum is the reported bottleneck and the gates are where the + # tube pinches for the volume. Length len(tetrahedra) - 1; empty for a + # single-tetrahedron path. + n = len(tetrahedra) + if n < 2: + return np.empty(0) + eb = self._edge_bottleneck + gates = np.empty(n - 1) + for k in range(n - 1): + i, j = int(tetrahedra[k]), int(tetrahedra[k + 1]) + key = (i, j) if i < j else (j, i) + g = eb.get(key) if eb is not None else None + if g is None: + shared = np.intersect1d(simp[i], simp[j], assume_unique=True) + if len(shared) != 3: + # not face-adjacent (should not happen on a graph path); + # fall back to the tighter of the two endpoints + g = float(min(vertex_radii[k], vertex_radii[k + 1])) + else: + g = self._edgeBottleneck(voronoi_vertices[i], + voronoi_vertices[j], shared, + points, vdw_radii) + gates[k] = g + return gates + + def calculateRadiusSpline(self, tetrahedra, voronoi_vertices, points, + vdw_radii, simp): + tetrahedra = np.asarray(tetrahedra) + # The per-vertex clearance is the same min buildSparseGraph already took + # over every simplex; read it back instead of recomputing it per path. + if self._vertex_clearance is not None: + radii = self._vertex_clearance[tetrahedra] + else: + vertices = voronoi_vertices[tetrahedra] + radii = np.array([self.calculateMaxRadius(v, points, vdw_radii, s) + for v, s in zip(vertices, simp[tetrahedra])]) + gates = self._pathGates(tetrahedra, voronoi_vertices, points, + vdw_radii, simp, radii) + return radii, gates + + def processChannel(self, tetrahedra, voronoi_vertices, points, vdw_radii, + simp): from scipy.interpolate import CubicSpline centers = voronoi_vertices[tetrahedra] - radii, bottleneck = self.calculate_radius_spline(tetrahedra, voronoi_vertices, points, vdw_radii, simp) - + radii, gates = self.calculateRadiusSpline(tetrahedra, + voronoi_vertices, + points, vdw_radii, simp) + bottleneck = float(np.min(gates)) if len(gates) else float(np.min(radii)) + t = np.arange(len(centers)) centerline_spline = CubicSpline(t, centers, bc_type='natural') - radius_spline = CubicSpline(t, radii, bc_type='natural') - - length = self.calculate_channel_length(centerline_spline) - volume = self.calculate_channel_volume(centerline_spline, radius_spline) + # The tube pinches at the gates, not at the wide circumcenters, so give + # the radius profile a knot at each gate (midway between its two + # vertices) carrying the gate clearance. The centerline keeps only the + # vertex knots; both splines share the same t domain, so the volume + # integral samples them consistently and the endpoints (hence the cap + # radii) are unchanged. + if len(gates): + knot_t = np.empty(2 * len(centers) - 1) + knot_t[0::2] = t + knot_t[1::2] = t[:-1] + 0.5 + knot_r = np.empty_like(knot_t) + knot_r[0::2] = radii + knot_r[1::2] = gates + radius_spline = CubicSpline(knot_t, knot_r, bc_type='natural') + else: + radius_spline = CubicSpline(t, radii, bc_type='natural') + + length = self.calculateChannelLength(centerline_spline) + volume = self.calculateChannelVolume(centerline_spline, radius_spline) return centerline_spline, radius_spline, length, bottleneck, volume - def find_biggest_tetrahedron(self, tetrahedra, voronoi_vertices, points, vdw_radii, simp): - radii = np.array([self.calculate_max_radius(voronoi_vertices[tetra], points, vdw_radii, simp[tetra]) for tetra in tetrahedra]) + def findBiggestTetrahedron(self, tetrahedra, voronoi_vertices, points, + vdw_radii, simp): + radii = np.array([self.calculateMaxRadius(voronoi_vertices[tetra], points, vdw_radii, simp[tetra]) for tetra in tetrahedra]) max_radius_index = np.argmax(radii) return tetrahedra[max_radius_index] - def get_end_tetrahedra(self, tetrahedra, voronoi_vertices, points, vdw_radii, simp, sparsity): - end_tetrahedra = [] - current_tetrahedron = self.find_biggest_tetrahedron(tetrahedra, voronoi_vertices, points, vdw_radii, simp) - end_tetrahedra.append(current_tetrahedron) - end_tetrahedra_set = {current_tetrahedron} - + def getEndTetrahedra(self, tetrahedra, voronoi_vertices, points, vdw_radii, + simp, sparsity): + # Greedy sparse sampling of the mouth (exit) tetrahedra: seed with the + # widest tetrahedron, then repeatedly add the widest tetrahedron that is + # still at least `sparsity` away from every already-selected one, until + # none qualify. Vectorized rewrite of the former O(N_exit x M^2) double + # loop (which called np.linalg.norm once per candidate/selected pair and + # re-scanned radii via findBiggestTetrahedron every pass): + # * the "far enough from all selected" test is exactly "running min + # distance to the selected set >= sparsity", so keep one min_dist + # array and fold in each new pick with a single vectorized norm; + # * inscribed radii are geometry-only, so precompute them once instead + # of recomputing find_biggest over the shrinking candidate set. + # Selection order and argmax first-tie-break match the original, so the + # returned end tetrahedra are identical. + tetrahedra = np.asarray(tetrahedra) + n = len(tetrahedra) + if n == 0: + return tetrahedra + + verts = voronoi_vertices[tetrahedra] # (n, 3) circumcenters + radii = np.array([ + self.calculateMaxRadius(voronoi_vertices[tetra], points, + vdw_radii, simp[tetra]) + for tetra in tetrahedra]) + + min_dist = np.full(n, np.inf) + selected = np.zeros(n, dtype=bool) + order = [] + + current = int(np.argmax(radii)) # widest tetrahedron (seed) while True: - found_tetrahedra = [] - for tetra in tetrahedra: - if tetra in end_tetrahedra_set: - continue - - all_far_enough = True - for selected_tetra in end_tetrahedra: - distance = np.linalg.norm(voronoi_vertices[selected_tetra] - voronoi_vertices[tetra]) - if distance < sparsity: - all_far_enough = False - break - - if all_far_enough: - found_tetrahedra.append(tetra) + order.append(current) + selected[current] = True + min_dist = np.minimum(min_dist, np.linalg.norm(verts - verts[current], axis=1)) - if not found_tetrahedra: + feasible = (min_dist >= sparsity) & ~selected # >= sparsity from every pick + if not feasible.any(): break - - biggest_tetrahedron = self.find_biggest_tetrahedron(found_tetrahedra, voronoi_vertices, points, vdw_radii, simp) - end_tetrahedra.append(biggest_tetrahedron) - end_tetrahedra_set.add(biggest_tetrahedron) + # widest feasible tetrahedron; np.argmax breaks ties toward the + # lowest index, matching the original input-order scan. + current = int(np.argmax(np.where(feasible, radii, -np.inf))) - return np.array(end_tetrahedra) + return tetrahedra[order] - def filter_cavities(self, cavities, min_depth): + def filterCavities(self, cavities, min_depth): return [cavity for cavity in cavities if cavity.depth >= min_depth] - def filter_channels_by_bottleneck(self, cavities, bottleneck): + def filterChannelsByBottleneck(self, cavities, bottleneck): for cavity in cavities: cavity.channels = [channel for channel in cavity.channels if channel.bottleneck >= bottleneck] - def filter_channels_by_volume(self, cavities, min_volume=None, max_volume=None): + def filterChannelsByVolume(self, cavities, min_volume=None, max_volume=None): """Filter channels by volume.""" for cavity in cavities: @@ -2634,7 +4440,8 @@ def filter_channels_by_volume(self, cavities, min_volume=None, max_volume=None): filtered_channels.append(channel) cavity.channels = filtered_channels - def filter_cavities_by_tetrahedra(self, cavities, min_tetrahedra=None, max_tetrahedra=None): + def filterCavitiesByTetrahedra(self, cavities, min_tetrahedra=None, + max_tetrahedra=None): """Filter cavities by cavity volume.""" filtered = [] @@ -2647,7 +4454,7 @@ def filter_cavities_by_tetrahedra(self, cavities, min_tetrahedra=None, max_tetra filtered.append(cavity) return filtered - def calculate_tetrahedron_volume(self, a, b, c, d): + def calculateTetrahedronVolume(self, a, b, c, d): return abs(np.dot(a - d, np.cross(b - d, c - d))) / 6.0 def calculate_cavity_volumes(self, cavities, simplices, coords): @@ -2658,10 +4465,10 @@ def calculate_cavity_volumes(self, cavities, simplices, coords): for tetra in cavity.tetrahedra: atom_ids = simplices[tetra] a, b, c, d = coords[atom_ids] - volume += self.calculate_tetrahedron_volume(a, b, c, d) + volume += self.calculateTetrahedronVolume(a, b, c, d) cavity.volume = volume - def filter_cavities_by_volume(self, cavities, min_volume=None, max_volume=None): + def filterCavitiesByVolume(self, cavities, min_volume=None, max_volume=None): """Filter cavities by approximate volume.""" filtered_cavities = [] @@ -2673,59 +4480,77 @@ def filter_cavities_by_volume(self, cavities, min_volume=None, max_volume=None): filtered_cavities.append(cavity) return filtered_cavities - def save_channels_to_pdb(self, cavities, filename, separate=False, num_samples=5): + def saveChannelsToPdb(self, channels, filename, separate=False, num_samples=5): + # ``channels`` is a flat list, already ordered by cost - that order is + # the order they are written and numbered here. Each channel is preceded + # by a REMARK reporting its length, bottleneck radius, curvature and cost. filename = str(filename) - - # All channels will be provided always when PDB/PQR will be created + + # All channels in a single file, one after another in list (cost) order. with open(filename, 'w') as pqr_file: atom_index = 1 - for cavity in cavities: - for channel in cavity.channels: - centerline_spline, radius_spline = channel.get_splines() + for channel_index, channel in enumerate(channels): + centerline_spline, radius_spline = channel.getSplines() + samples = len(channel.tetrahedra) * num_samples + t = np.linspace(centerline_spline.x[0], centerline_spline.x[-1], samples) + centers = centerline_spline(t) + radii = radius_spline(t) + + pqr_file.write(self._channelRemark(channel_index, channel)) + pdb_lines = [] + # Each channel gets its own residue number so the channels stay + # separable at the record level, matching saveCavitiesToPdb. + for i, (x, y, z, radius) in enumerate(zip(centers[:, 0], centers[:, 1], centers[:, 2], radii), start=atom_index): + pdb_lines.append("ATOM %5d H FIL T%4d %8.3f%8.3f%8.3f%6.2f%6.2f\n" % (i, channel_index + 1, x, y, z, 1.00, radius)) + + # Bond consecutive samples of THIS channel only, using the global + # atom serial numbers (start=atom_index). No CONECT spans two + # channels, so each one is a separate strand in the viewer. + for i in range(atom_index, atom_index + samples - 1): + pdb_lines.append("CONECT%5d%5d\n" % (i, i + 1)) + + pqr_file.writelines(pdb_lines) + pqr_file.write("\n") + atom_index += samples + + # When separate is set to True also separate PDB/PQR files will be + # created, one per channel, numbered by the same cost order. + if separate: + for channel_index, channel in enumerate(channels): + # TODO channel suffix is rather long, making it hard to read in pymol, shorten or user defined? + channel_filename = filename.replace('.pqr', '_chl{0}.pqr'.format(channel_index)) + channel_filename = channel_filename.replace('.pdb', '_chl{0}.pdb'.format(channel_index)) + + with open(channel_filename, 'w') as pqr_file: + atom_index = 1 + centerline_spline, radius_spline = channel.getSplines() samples = len(channel.tetrahedra) * num_samples t = np.linspace(centerline_spline.x[0], centerline_spline.x[-1], samples) centers = centerline_spline(t) radii = radius_spline(t) + pqr_file.write(self._channelRemark(channel_index, channel)) pdb_lines = [] for i, (x, y, z, radius) in enumerate(zip(centers[:, 0], centers[:, 1], centers[:, 2], radii), start=atom_index): - pdb_lines.append("ATOM %5d H FIL T 1 %8.3f%8.3f%8.3f%6.2f%6.2f\n" % (i, x, y, z, 1.00, radius)) + pdb_lines.append("ATOM %5d H FIL T%4d %8.3f%8.3f%8.3f%6.2f%6.2f\n" % (i, channel_index + 1, x, y, z, 1.00, radius)) - for i in range(1, samples): + for i in range(atom_index, atom_index + samples - 1): pdb_lines.append("CONECT%5d%5d\n" % (i, i + 1)) - + pqr_file.writelines(pdb_lines) - pqr_file.write("\n") - atom_index += samples - - # When separate is set to True also separate PDB/PQR files will be created - if separate: - channel_index = 0 - for cavity in cavities: - for channel in cavity.channels: - channel_filename = filename.replace('.pqr', '_channel{0}.pqr'.format(channel_index)) - - with open(channel_filename, 'w') as pqr_file: - atom_index = 1 - centerline_spline, radius_spline = channel.get_splines() - samples = len(channel.tetrahedra) * num_samples - t = np.linspace(centerline_spline.x[0], centerline_spline.x[-1], samples) - centers = centerline_spline(t) - radii = radius_spline(t) - - pdb_lines = [] - for i, (x, y, z, radius) in enumerate(zip(centers[:, 0], centers[:, 1], centers[:, 2], radii), start=atom_index): - pdb_lines.append("ATOM %5d H FIL T 1 %8.3f%8.3f%8.3f%6.2f%6.2f\n" % (i, x, y, z, 1.00, radius)) - - for i in range(1, samples): - pdb_lines.append("CONECT%5d%5d\n" % (i, i + 1)) - - pqr_file.writelines(pdb_lines) - - channel_index += 1 - - - def save_cavities_to_pdb(self, cavities, vertices, filename, separate=False): + + @staticmethod + def _channelRemark(channel_index, channel): + """One-line PQR/PDB REMARK with a channel's basic geometry: + length, bottleneck radius, curvature and Dijkstra cost.""" + curv = 'n/a' if np.isnan(channel.curvature) else "%.3f" % channel.curvature + cost = 'n/a' if channel.cost is None else "%.4g" % channel.cost + return ("REMARK Channel %d length=%.3f A bottleneck=%.3f A " + "curvature=%s cost=%s\n" % ( + channel_index, channel.length, channel.bottleneck, curv, cost)) + + + def saveCavitiesToPdb(self, cavities, vertices, filename, separate=False): """Save surface cavities to a PDB/PQR file as dummy atoms.""" filename = str(filename) @@ -2768,66 +4593,217 @@ def save_cavities_to_pdb(self, cavities, vertices, filename, separate=False): cavity_index += 1 - def calculate_channel_length(self, centerline_spline): + def calculateChannelLength(self, centerline_spline): t_values = np.linspace(centerline_spline.x[0], centerline_spline.x[-1], len(centerline_spline.x) * 10) points = centerline_spline(t_values) diffs = np.diff(points, axis=0) lengths = np.linalg.norm(diffs, axis=1) return np.sum(lengths) - def calculate_channel_volume(self, centerline_spline, radius_spline): - import warnings - from scipy.integrate import quad, IntegrationWarning - - warnings.filterwarnings("ignore", category=IntegrationWarning) - + def calculateChannelVolume(self, centerline_spline, radius_spline): + # Tube volume V = \int pi r(t)^2 |x'(t)| dt evaluated by vectorized + # composite Simpson instead of an adaptive scipy.quad that called the + # integrand thousands of times per channel. The centerline/radius are + # piecewise cubic, so a uniform grid scaled to the number of spline + # segments (~128 points/segment) reaches ~1e-8 relative accuracy - better + # than quad's default tolerance - at a fraction of the cost. + from scipy.integrate import simpson + t_min = centerline_spline.x[0] t_max = centerline_spline.x[-1] - - def differential_volume(t): - r = radius_spline(t) - area = np.pi * r**2 - dx_dt = centerline_spline(t, 1) - centerline_derivative = np.linalg.norm(dx_dt) - return area * centerline_derivative - - volume, error = quad(differential_volume, t_min, t_max) - + n_segments = max(1, len(centerline_spline.x) - 1) + n = n_segments * 128 + 1 + + t = np.linspace(t_min, t_max, n) + r = radius_spline(t) + speed = np.linalg.norm(centerline_spline(t, 1), axis=1) + volume = simpson(np.pi * r ** 2 * speed, x=t) + r_start = radius_spline(t_min) r_end = radius_spline(t_max) - + hemisphere_volume_start = (2/3) * np.pi * r_start**3 hemisphere_volume_end = (2/3) * np.pi * r_end**3 - + total_volume = volume + hemisphere_volume_start + hemisphere_volume_end - + return total_volume - def set_starting_tetrahedra_from_point(self, cavities, vertices, start_point): + def selectSeedTetrahedron(self, cavity, vertices, points, vdw_radii, simp, + neighbors, sp, search_radius): + '''Map `sp` to the seed tetrahedron of one cavity. + + The tetrahedron nearest `sp` (the anchor) is frequently a tight one, and every + channel of the cavity leaves through it: its inscribed radius then caps all of + their bottlenecks, and the shared first links show up as one common bottleneck + at the joint beginning of the bundle. So the anchor only says where to look. + The seed is the widest (largest inscribed radius) tetrahedron of this cavity + that lies within `search_radius` of `sp`, is no shallower than the anchor, and + is reachable from the anchor through the tetrahedra within `search_radius`. + Reachability is over the adjacency of the cleared tetrahedra, which is free + space, so the seed can only move through the void the start point sits in and + never hops across a wall into a lobe that merely passes nearby; the depth floor + keeps it from sliding outward towards the mouth, where tetrahedra are wide but + no longer inside the site. Note that the floor filters the seed, not the walk: + a marginally shallower cell in between must not wall off the wider region + behind it. + + `search_radius` <= 0 restores the plain nearest-vertex seed. + + :returns: dict of the seed and anchor properties (`seed`, `anchor`, and their + `_vertex`, `_distance` from `sp`, inscribed `_radius` and `_depth`), plus + the number of tetrahedra `searched` and how many of them were `eligible`''' + + from collections import deque + + depths = cavity.tetrahedra_depths + tet = np.asarray(cavity.tetrahedra) + d2 = np.sum((vertices[tet] - sp) ** 2, axis=1) + anchor = int(tet[int(np.argmin(d2))]) + + def properties(tetra): + return dict( + vertex=vertices[tetra], + distance=float(np.linalg.norm(vertices[tetra] - sp)), + radius=float(self.calculateMaxRadius( + vertices[tetra], points, vdw_radii, simp[tetra])), + depth=float(depths.get(tetra, 0.0))) + + def report(seed, searched, eligible): + info = {'seed': seed, 'anchor': anchor, + 'searched': searched, 'eligible': eligible} + for name, tetra in (('seed', seed), ('anchor', anchor)): + for key, value in properties(tetra).items(): + info['{0}_{1}'.format(name, key)] = value + return info + + if not search_radius or search_radius <= 0: + return report(anchor, 1, 1) + + near = set(int(t) for t, close in zip(tet, d2 <= search_radius ** 2) if close) + + # BFS from the anchor, staying inside the sphere and inside the cavity. + reachable = [anchor] + seen = {anchor} + queue = deque([anchor]) + while queue: + current = queue.popleft() + for neighbor in neighbors[current]: + neighbor = int(neighbor) + if neighbor in near and neighbor not in seen: + seen.add(neighbor) + reachable.append(neighbor) + queue.append(neighbor) + + anchor_depth = depths.get(anchor, 0) + eligible = [t for t in reachable if depths.get(t, 0) >= anchor_depth] + + reach = np.array(eligible, dtype=np.intp) + atoms = simp[reach] + clearance = (np.linalg.norm(points[atoms] - vertices[reach][:, None, :], axis=2) + - vdw_radii[atoms]) + radii = clearance.min(axis=1) + # The anchor is eligible and comes first (BFS order), so argmax ties to it. + best = int(np.argmax(radii)) + + return report(int(reach[best]), len(reachable), len(eligible)) + + def setStartingTetrahedraFromPoint(self, cavities, vertices, start_point, + points, vdw_radii, simp, neighbors, + restrict=False, search_radius=5.0): '''Set starting tetrahedra using a user-defined 3D point. - The starting tetrahedron is selected as the one whose Voronoi vertex is closest - to `start_point` (Euclidean distance). - + The starting tetrahedron of a cavity is the widest one `selectSeedTetrahedron` + finds in the neighbourhood of `start_point`; with ``search_radius=0`` it is + simply the one whose Voronoi vertex is closest to `start_point`. + :arg cavities: list of cavity objects :arg vertices: Voronoi vertices (array of shape (n, 3)) - :arg start_point: point [x, y, z] in Å (list/tuple/ndarray of length 3)''' - + :arg start_point: point [x, y, z] in Å (list/tuple/ndarray of length 3) + :arg points: atom coordinates (array of shape (n_atoms, 3)), used to compute + the inscribed radius of candidate tetrahedra + :arg vdw_radii: per-atom van der Waals radii (array of shape (n_atoms,)) + :arg simp: simplices (tetrahedron -> its 4 atom indices) + :arg neighbors: tetrahedron adjacency (tetrahedron -> its 4 neighbours, -1 none) + :arg restrict: if True, only the single cavity whose closest tetrahedron is + globally nearest to `start_point` is seeded and returned, so channels are + computed exclusively for the region around `start_point`. If False (default), + every cavity is seeded with its own seed tetrahedron and all cavities are + returned unchanged. + :type restrict: bool + :arg search_radius: radius, in Angstrom, of the neighbourhood of `start_point` + searched for a wider seed. 0 disables the search. + :type search_radius: float + :returns: list of cavities to search: all cavities when `restrict` is False, the + single selected cavity when `restrict` is True, or an empty list if no cavity + has any tetrahedra''' + sp = np.asarray(start_point, dtype=float).reshape(3,) - for cavity in cavities: + best_cavity = None + best_info = None + + for i, cavity in enumerate(cavities): tet = cavity.tetrahedra if tet is None or len(tet) == 0: continue - # Voronoi vertex per tetrahedron: vertices[tetra_id] -> (x,y,z) - v = vertices[tet] - d2 = np.sum((v - sp) ** 2, axis=1) - chosen = tet[int(np.argmin(d2))] - - cavity.set_starting_tetrahedron(np.array([chosen])) - - - def trim_cavities_by_depth(self, cavities, max_depth): + info = self.selectSeedTetrahedron( + cavity, vertices, points, vdw_radii, simp, neighbors, sp, search_radius) + + if not restrict: + cavity.setStartingTetrahedron(np.array([info['seed']])) + self.reportSeedTetrahedron(info, search_radius, cavity_index=i) + + # The cavity is still chosen by proximity to start_point: widening moves the + # seed inside a cavity, it must never decide between cavities. + if best_info is None or info['anchor_distance'] < best_info['anchor_distance']: + best_info = info + best_cavity = cavity + + if not restrict: + return cavities + + if best_cavity is None: + _warn("start_point was provided but no cavity contains any " + "tetrahedron; no channels will be computed.") + return [] + + best_cavity.setStartingTetrahedron(np.array([best_info['seed']])) + self.reportSeedTetrahedron(best_info, search_radius) + LOGGER.info(" restricting the channel search to the cavity that contains it " + "({0} tetrahedra, depth {1:.1f} A).".format(len(best_cavity.tetrahedra), + float(best_cavity.depth))) + + return [best_cavity] + + def reportSeedTetrahedron(self, info, search_radius, cavity_index=None): + '''Log the seed tetrahedron `selectSeedTetrahedron` picked, and, when it is not + the one nearest the start point, the anchor it replaced -- the two radii are what + tell the user whether the seed was capping the bottlenecks of the cavity.''' + + where = '' if cavity_index is None else ' of cavity {0}'.format(cavity_index) + LOGGER.info("start_point seeded at tetrahedron {0}{1} (Voronoi vertex at " + "[{2:.3f}, {3:.3f}, {4:.3f}], {5:.3f} A from start_point, inscribed radius " + "{6:.3f} A, depth {7:.1f} A)." + .format(info['seed'], where, info['seed_vertex'][0], info['seed_vertex'][1], + info['seed_vertex'][2], info['seed_distance'], info['seed_radius'], + info['seed_depth'])) + + if info['seed'] != info['anchor']: + LOGGER.info(" widened from the nearest tetrahedron {0} ({1:.3f} A away, " + "inscribed radius {2:.3f} A, depth {3:.1f} A), the widest of the {4} tetrahedra " + "no shallower than it among the {5} reachable within {6:.1f} A; seeding " + "the narrow one would have capped every channel here at its radius." + .format(info['anchor'], info['anchor_distance'], info['anchor_radius'], + info['anchor_depth'], info['eligible'], info['searched'], + float(search_radius))) + elif search_radius and search_radius > 0: + LOGGER.info(" already the widest of the {0} tetrahedra no shallower than " + "it among the {1} reachable within {2:.1f} A." + .format(info['eligible'], info['searched'], float(search_radius))) + + + def trimCavitiesByDepth(self, cavities, max_depth): """Filtering cavities by max_depth.""" for cavity in cavities: