From d6c025bd32b82bbd83154bc3509847557a24e707 Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 12:21:58 -0400 Subject: [PATCH 01/16] changes sum to the builtin sum to prevent issue with numpy depreciation of np.sum(generator) --- prody/proteins/channels.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index cbc182f40..dd9849cc2 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -27,6 +27,8 @@ from .compare import * from prody.measure import calcTransformation, calcDistance, calcRMSD, superpose +import builtins +sum = builtins.sum __all__ = ['getVmdModel', 'calcChannels', 'calcChannelsMultipleFrames', 'getChannelParameters', 'getChannelAtoms', 'showChannels', 'showCavities', From f02cb20de121d73fca5850eab48f7f8b892bf208 Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 12:29:34 -0400 Subject: [PATCH 02/16] changed calcChannels to match new calculator functions --- prody/proteins/channels.py | 68 +++++++++++++++++++++----------------- 1 file changed, 37 insertions(+), 31 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index dd9849cc2..521e7cbb9 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -377,7 +377,7 @@ def showCavities(surface, show_surface=False): vis.destroy_window() -def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3, r2=1.25, min_depth=10, bottleneck=1, sparsity=15): +def calcChannels(atoms, r1=3, r2=1.25, min_depth=10, bottleneck=1, sparsity=15, start_point=None, **kwargs): """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 @@ -403,11 +403,6 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 are saved in a single PDB 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 Å. - :type start_point: list, tuple, or ndarray (length 3), 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. :type r1: float @@ -446,8 +441,6 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 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]) - 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) """ @@ -474,7 +467,7 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 from pathlib import Path else: from pathlib2 import Path - + if start_point is not None: if not isListLike(start_point): raise TypeError("start_point must be a list/tuple/ndarray with three numeric values") @@ -487,19 +480,24 @@ 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])) - - + + output_path = kwargs.pop('output_path',None) + separate = kwargs.pop('separate',False) + filename = kwargs.pop('filename',None) + # Initialize the channel detection helper with the provided parameters. calculator = ChannelCalculator(atoms, r1, r2, min_depth, bottleneck, sparsity) + # Filter the input atoms to exclude heteroatoms and hydrogens, since channels are defined by the protein backbone. atoms = atoms.select('not hetero and noh') # Excluding hydrogens coords = atoms.getCoords() vdw_radii = calculator.get_vdw_radii(atoms.getElements()) - + + # Compute Delaunay triangulation and Voronoi tessellation from the atomic coordinates. dela = Delaunay(coords) - voro = Voronoi(coords) - - s_prt = State(dela.simplices, dela.neighbors, voro.vertices) + verts = calculator.calc_circumcenters(dela) + # The state object holds the current simplices, neighbor relations, and Voronoi vertices. + s_prt = State(dela.simplices, dela.neighbors, verts) #voro.vertices) if PY3K: s_tmp = State(*s_prt.get_state()) @@ -507,12 +505,14 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 else: s_tmp = apply(State, s_prt.get_state()) s_prv = State(None, None, None) - + + # Iteratively delete simplices until the surface state converges. + while True: s_prv.set_state(*s_tmp.get_state()) if PY3K: - #s_tmp.set_state(*calculator.delete_simplices3d(coords, *(s_tmp.get_state() + [vdw_radii, r1, True]))) + # Remove simplices that are outside the molecule surface for the outer radius threshold. s_tmp.set_state(*calculator.delete_simplices3d(coords, *(s_tmp.get_state() + tuple([vdw_radii, r1, True])))) else: tmp_state = calculator.delete_simplices3d(coords, *(s_tmp.get_state() + [vdw_radii, r1, True])) @@ -520,35 +520,41 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 if s_tmp == s_prv: break - + # State after surface deletion is the outer surface state. s_srf = State(*s_tmp.get_state()) - #s_inr = State(*calculator.delete_simplices3d(coords, *(s_srf.get_state() + [vdw_radii, r2, False]))) + # Perform a second deletion with the inner radius threshold to identify inner surface simplices. s_inr = State(*calculator.delete_simplices3d(coords, *(s_srf.get_state() + tuple([vdw_radii, r2, False])))) + # Separate the first surface layer and the second interior layer for cavity detection. 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())) - + # Find connected regions of tetrahedra, then identify cavities connected to the surface. 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) - + # Determine the deepest entry points and filter by minimum depth. calculator.find_deepest_tetrahedra(c_surface_cavities, s_clr.neigh) - 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) + if not c_filtered_cavities: + LOGGER.warning("No cavities found.") + return [], [] merged_cavities = calculator.merge_cavities(c_filtered_cavities, s_clr.simp) - + + # For each cavity, run Dijkstra search to find probable channel paths. + simplices, neighbors, vertices = s_clr.get_state() + #build graph for dijkstra + graph = calculator.build_sparse_graph(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.dijkstra(cavity, graph, simplices, neighbors, vertices, coords, vdw_radii) + # Remove channels narrower than the specified bottleneck radius. calculator.filter_channels_by_bottleneck(c_filtered_cavities, bottleneck) channels = [channel for cavity in c_filtered_cavities for channel in cavity.channels] - + no_of_channels = len(channels) LOGGER.info("Detected " + str(no_of_channels) + " channels.") - + if output_path: + if filename: + output_path = output_path + filename output_path = Path(output_path) if output_path.is_dir(): @@ -564,7 +570,7 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, r1=3 calculator.save_channels_to_pdb(c_filtered_cavities, output_path, separate) else: LOGGER.info("No output path given.") - + # Return the detected channels and surface state data for further analysis or visualization. return channels, [coords, s_srf.simp, merged_cavities, s_clr.simp] From 0a88ea13bf02dd18fa7e0e3a5aa14c21abd7f3ae Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 12:33:02 -0400 Subject: [PATCH 03/16] added calcChannelsMultipleAtomGroups which runs calcChannels for multiple unrelated atomgroups --- prody/proteins/channels.py | 40 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 521e7cbb9..bf74ca7e5 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -573,6 +573,46 @@ def calcChannels(atoms, r1=3, r2=1.25, min_depth=10, bottleneck=1, sparsity=15, # Return the detected channels and surface state data for further analysis or visualization. return channels, [coords, s_srf.simp, merged_cavities, s_clr.simp] +def _ag_worker(args): + atoms,kwargs,i,filename= args + # LOGGER.info("Processing atom group " + str(filename)) + start = time.perf_counter() + try: + channels, surfaces = calcChannels(atoms,filename = filename, **kwargs) + if not channels: + return i,[],[],0 + total = time.perf_counter() - start + return i, channels, surfaces, total + except: + return i, [], [], 0 + + +def calcChannelsMultipleAtomGroups(atomgroups, **kwargs): + from multiprocessing import Pool, cpu_count + filenames = kwargs.pop('filenames',[None]*len(atomgroups)) + max_proc = kwargs.pop('max_proc',None) + if max_proc is None: + max_proc = max(1, cpu_count()//2) + + if max_proc == 1: + results = [_ag_worker((ag,kwargs,i,filenames[i])) for i, ag in enumerate(atomgroups)] + else: + tasks = ((ag,kwargs,i,filenames[i]) for i, ag in enumerate(atomgroups)) + with Pool(processes=max_proc) as pool: + results = pool.map(_ag_worker, tasks,chunksize=builtins.max(1, len(atomgroups)//(max_proc*4))) + + results.sort(key=lambda x: x[0]) + + channels_all = [r[1] for r in results] + for channels in channels_all: + for channel in channels: + channel.build_splines() + surfaces_all = [r[2] for r in results] + times_all = [r[3] for r in results] + failed = [filenames[r[0]] for r in results if (not r[1] and not r[2])] + if failed: + LOGGER.warning(f"WARNING: {len(failed)} proteins failed or No Channels Detected: {', '.join(str(f) for f in failed)}") + return channels_all, surfaces_all, times_all 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. From 2665d4a002d55c1e540e30aefe83aeb433790102 Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 12:34:30 -0400 Subject: [PATCH 04/16] updated calcChannelsMultipleFrames to use multiprocessing POOL --- prody/proteins/channels.py | 205 +++++++++++++++++++++++-------------- 1 file changed, 128 insertions(+), 77 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index bf74ca7e5..d1a193c7a 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -614,121 +614,172 @@ def calcChannelsMultipleAtomGroups(atomgroups, **kwargs): LOGGER.warning(f"WARNING: {len(failed)} proteins failed or No Channels Detected: {', '.join(str(f) for f in failed)}") return channels_all, surfaces_all, times_all -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. +def _frame_worker(args): + atoms, coords, kwargs, frame_id = args - 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. + atoms_local = atoms.copy() + atoms_local.setCoords(coords) - :param atoms: Atomic data or object containing atomic coordinates and methods for accessing them. - :type atoms: object + try: + channels, surfaces = calcChannels(atoms_local, **kwargs) + if not channels: + return frame_id,[],[] + except: + return frame_id, [], [] + return frame_id, channels, surfaces + +def _model_worker(args): + """Process a single model from multi-model PDB""" + atoms, kwargs, frame_id = args + start_frame = kwargs.pop(0,"start_frame") + frame_idx = frame_id + start_frame + LOGGER.info("Model: {0}".format(frame_idx)) + # Make a copy to avoid state conflicts + atoms_copy = atoms.copy() + atoms_copy.setACSIndex(frame_idx) + + try: + channels, surfaces = calcChannels(atoms_copy, **kwargs) + if not channels: + return frame_id,[],[] + except: + return frame_id, [], [] + + return frame_id, channels, surfaces + +def calcChannelsMultipleFrames(atoms, trajectory, **kwargs): + """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. - :param trajectory: Trajectory object containing multiple frames or a multi-model PDB file. - :type trajectory: Atomic or Ensemble object + 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. + :type atoms: `Atoms` 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. + :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. :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. + :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. :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 Å. - :type start_point: list, tuple, or ndarray (length 3), 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. + :type r1: float - :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. - See the available parameters in calcChannels(). - :type kwargs: dict + :param r2: The second radius threshold used to define the inner surface of the channels. Default is 1.25. + :type r2: float - :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 + :param min_depth: The minimum depth a cavity must have to be considered as a channel. Default is 10. + :type min_depth: float - 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, start_point=[-10.353, -0.133, 5.608]) """ + :param bottleneck: The minimum allowed bottleneck size (narrowest point) for the channels. Default is 1. + :type bottleneck: 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 + + :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. + :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. + + Example usage: + channels, surface = calcChannels(atoms, output_path="channels", separate=True) - if PY3K: - if not checkAndImport('pathlib'): - errorMsg = 'To run showChannels, please install open3d.' - raise ImportError(errorMsg) - - from pathlib import Path - else: - if not checkAndImport('pathlib2'): - errorMsg = 'To run showChannels, please install pathlib2 for Python 2.7.' - raise ImportError(errorMsg) - - from pathlib2 import Path - + 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) """ + from multiprocessing import Pool, cpu_count + try: coords = getCoords(atoms) except AttributeError: try: checkCoords(coords) except TypeError: - raise TypeError('coords must be an object with `getCoords` method') + raise TypeError('coords must be an object with `getCoords` method') + + max_proc = kwargs.pop('max_proc',None) + if max_proc is None: + max_proc = max(1, cpu_count()//2) - channels_all = [] - surfaces_all = [] start_frame = kwargs.pop('start_frame', 0) stop_frame = kwargs.pop('stop_frame', -1) - - if output_path: - output_path = Path(output_path) - if output_path.suffix == ".pqr": - output_path = output_path.with_suffix('') if trajectory is not None: if isinstance(trajectory, Atomic): trajectory = Ensemble(trajectory) - nfi = trajectory._nfi - trajectory.reset() - if stop_frame == -1: traj = trajectory[start_frame:] else: traj = trajectory[start_frame:stop_frame+1] - atoms_copy = atoms.copy() - for j0, frame0 in enumerate(traj, start=start_frame): - LOGGER.info("Frame: {0}".format(j0)) - atoms_copy.setCoords(frame0.getCoords()) - if output_path: - channels, surfaces = calcChannels(atoms_copy, str(output_path) + "{0}.pqr".format(j0), separate, start_point=start_point, **kwargs) - else: - channels, surfaces = calcChannels(atoms_copy, start_point=start_point, **kwargs) - channels_all.append(channels) - surfaces_all.append(surfaces) - trajectory._nfi = nfi + if max_proc == 1: + results = [_frame_worker((atoms,frame.getCoords(),kwargs,i)) for i, frame in enumerate(traj)] + else: + tasks = ((atoms,frame.getCoords(),kwargs,i) for i, frame in enumerate(trajectory)) + with Pool(processes=max_proc) as pool: + results = pool.map(_frame_worker, tasks,chunksize=builtins.max(1, len(traj)//(max_proc*4))) + results.sort(key=lambda x: x[0]) + + channels_all = [r[1] for r in results] + surfaces_all = [r[2] for r in results] + times_all = [r[3] for r in results] + + return channels_all, surfaces_all, times_all else: - if atoms.numCoordsets() > 1: - for i in range(len(atoms.getCoordsets()[start_frame:stop_frame])): - LOGGER.info("Model: {0}".format(i+start_frame)) - atoms.setACSIndex(i+start_frame) - if output_path: - channels, surfaces = calcChannels(atoms, str(output_path) + "{0}.pqr".format(i+start_frame), separate, start_point=start_point, **kwargs) - else: - channels, surfaces = calcChannels(atoms, start_point=start_point, **kwargs) - channels_all.append(channels) - surfaces_all.append(surfaces) + if atoms.numCoordsets() > 1: + if stop_frame == -1: + num_models = len(atoms.getCoordsets()[start_frame:]) + else: + num_models = len(atoms.getCoordsets()[start_frame:stop_frame+1]) + + if max_proc == 1: + results = [_model_worker((atoms,kwargs,i)) for i in range(num_models)] + else: + tasks = ((atoms,kwargs,i) for i in range(num_models)) + + with Pool(processes=max_proc) as pool: + results = pool.map(_model_worker, tasks,chunksize=builtins.max(1, len(traj)//(max_proc*4))) + + results.sort(key=lambda x: x[0]) + + channels_all = [r[1] for r in results] + surfaces_all = [r[2] for r in results] + times_all = [r[3] for r in results] + + return channels_all, surfaces_all, times_all + else: LOGGER.info("Include trajectory or use multi-model PDB file.") - return channels_all, surfaces_all - def parseParameters(channels, **kwargs): """Extracts and returns the lengths, bottlenecks, and volumes of each channel in a given list of channels. """ From b64bfcc8c64eb09bb51a8e92effe813476dd4bdf Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 12:42:10 -0400 Subject: [PATCH 05/16] updated Channel class to be compatible with functions using multiprocressing Poolsuch as _ag_worker and _frame_worker --- prody/proteins/channels.py | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index d1a193c7a..b94baa6a8 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -1250,17 +1250,35 @@ def write_merge_surf_pdb(merged_surface, filename, nr_pdbs): class Channel: - def __init__(self, tetrahedra, centerline_spline, radius_spline, length, bottleneck, volume): + def __init__(self, tetrahedra, centerline_spline, radius_spline, length, bottleneck, volume, centers, radii): + # A detected channel path through the interior of the protein. self.tetrahedra = tetrahedra self.centerline_spline = centerline_spline self.radius_spline = radius_spline self.length = length self.bottleneck = bottleneck self.volume = volume - + self.centers = centers + self.radii = radii + def get_splines(self): + # Return the parametric centerline and radius curves for visualization or analysis. return self.centerline_spline, self.radius_spline + def __getstate__(self): + state = self.__dict__.copy() + state["centerline_spline"] = None + state["radius_spline"] = None + return state + + def build_splines(self): + from scipy.interpolate import CubicSpline + centers = self.centers + radii = self.radii + t = np.arange(len(centers)) + self.centerline_spline = CubicSpline(t, centers, bc_type='natural') + self.radius_spline = CubicSpline(t, radii, bc_type='natural') + class State: def __init__(self, simplices, neighbors, vertices): From ec6abb88ab4c13b4420d2aa443526ace3ce5167b Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 12:44:25 -0400 Subject: [PATCH 06/16] vectorized sphere_fit and delete_simplices3d so that spherefit is calculated for all relevant tetrahedra in the vector, rather than a for loop --- prody/proteins/channels.py | 52 ++++++++++++++++++-------------------- 1 file changed, 24 insertions(+), 28 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index b94baa6a8..881352a10 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -1335,38 +1335,34 @@ def __init__(self, atoms, r1=3, r2=1.25, min_depth=10, bottleneck=1, sparsity=15 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) - + def sphere_fit(self,points, vertices, simplices, vdw_radii, r): + d = np.linalg.norm(points[simplices]-vertices[:,None,:],axis=-1) + d_sum = np.sum(d, axis=1) + r_sum = 4*r + np.sum(vdw_radii[simplices], axis=1) + return d_sum >= r_sum + 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 + + simplices = np.ascontiguousarray(simplices) + vertices = np.ascontiguousarray(vertices) + n = len(simplices) + if surface: + boundary = np.any(neighbors == -1, axis=1) + keep_mask = np.ones(n, dtype=bool) + keep_mask[boundary] = ~self.sphere_fit(points,vertices[boundary],simplices[boundary],vdw_radii,r) + else: + keep_mask = self.sphere_fit(points,vertices,simplices,vdw_radii,r) + simp = simplices[keep_mask] + neigh = neighbors[keep_mask] + verti = vertices[keep_mask] + mapping = -np.ones(n, dtype=int) + mapping[np.where(keep_mask)[0]] = np.arange(np.sum(keep_mask)) + mask = neigh != -1 + neigh[mask] = mapping[neigh[mask]] + return simp, neigh, verti def delete_section(self, simplices_subset, simplices, neighbors, vertices, reverse=False): From 3928489aeb16fdd3d96aac9d79ba5d9e3ea5217d Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 13:03:38 -0400 Subject: [PATCH 07/16] updated delete_section for faster mask creation and filtering --- prody/proteins/channels.py | 37 +++++++++++++++---------------------- 1 file changed, 15 insertions(+), 22 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 881352a10..23f07e75e 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -1366,33 +1366,26 @@ def delete_simplices3d(self, points, simplices, neighbors, vertices, vdw_radii, 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]) + matches = (simplices[:, None,:] == simplices_subset[None,:,:]).all(axis=2).any(axis=1) + if reverse: + keep_mask = matches + else: + keep_mask = ~matches + + delete_mask = ~keep_mask - simp, neigh, verti = map(np.array, [simp, neigh, verti]) - deleted = np.array(deleted) + simp = simplices[keep_mask] + neigh = neighbors[keep_mask] + verti = vertices[keep_mask] + + deleted = np.flatnonzero(delete_mask) mask = np.isin(neigh, deleted) neigh[mask] = -1 - for i in reversed(deleted): - neigh = np.where((neigh > i) & (neigh != -1), neigh - 1, neigh) + valid = neigh >= 0 + + neigh[valid] -= np.searchsorted(deleted, neigh[valid]) return simp, neigh, verti From 236b013153e8a6774dafc4307b5359ee60a6b9a1 Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 13:06:57 -0400 Subject: [PATCH 08/16] updated surface_layer to take advantage of contiguousarrays and views --- prody/proteins/channels.py | 72 +++++++++++++++++++++----------------- 1 file changed, 39 insertions(+), 33 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 23f07e75e..b5c399104 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -1406,40 +1406,46 @@ 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) - ] - - 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) - ] - - 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) - ] + + surface_mask = np.any(shape_neighbors == -1, axis=1) + + surface_simplices = shape_simplices[surface_mask] + surface_neighbors = shape_neighbors[surface_mask] + interior_simplices = shape_simplices[~surface_mask] + + row_dtype = np.dtype( + (np.void,shape_simplices.dtype.itemsize * shape_simplices.shape[1])) + + filtered_view = ( + np.ascontiguousarray(filtered_simplices).view(row_dtype).ravel()) + + surface_view = ( + np.ascontiguousarray(surface_simplices).view(row_dtype).ravel()) + + surface_filter_mask = np.isin(surface_view, filtered_view) + + filtered_surface_simplices = surface_simplices[surface_filter_mask] + filtered_surface_neighbors = surface_neighbors[surface_filter_mask] + filtered_surface_neighbors = np.unique(filtered_surface_neighbors) + filtered_surface_neighbors = filtered_surface_neighbors[filtered_surface_neighbors >= 0] + + interior_view = (np.ascontiguousarray(interior_simplices).view(row_dtype).ravel()) + + interior_filter_mask = np.isin(interior_view, filtered_view) + + filtered_interior_simplices = interior_simplices[interior_filter_mask] + + surface_layer_neighbor_simplices = (shape_simplices[filtered_surface_neighbors]) + + neighbor_view = (np.ascontiguousarray(surface_layer_neighbor_simplices).view(row_dtype).ravel()) + + filtered_interior_view = (np.ascontiguousarray(filtered_interior_simplices).view(row_dtype).ravel()) + + second_layer_mask = np.isin(filtered_interior_view,neighbor_view) + + second_layer = filtered_interior_simplices[second_layer_mask] + return filtered_surface_simplices, second_layer From eb72e3b54365e8998d1b84addbcd0c3a7df88935 Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 13:09:30 -0400 Subject: [PATCH 09/16] updated dijkstra to use scipy algorithms and to build paths as nodes and branches --- prody/proteins/channels.py | 132 +++++++++++++++++++++++-------------- 1 file changed, 82 insertions(+), 50 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index b5c399104..1b0fadfcb 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -1521,61 +1521,93 @@ def find_deepest_tetrahedra(self, cavities, neighbors): cavity.set_starting_tetrahedron(np.array([deepest_tetrahedron])) cavity.set_depth(max_depth) - 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) + def build_sparse_graph(self, simplices, neighbors, vertices, points, vdw_radii): + + import numpy as np + from scipy.sparse import csr_matrix + + tetra_points = points[simplices] + distances = np.linalg.norm( + tetra_points - vertices[:, None, :], + axis=2) + bottleneck = np.min(distances - vdw_radii[simplices],axis=1) + + rows = [] + cols = [] + data = [] + + b = 1e-3 + for tetra, neighs in enumerate(neighbors): + for neigh in neighs: + if neigh == -1: + continue + l = np.linalg.norm( + vertices[tetra] - vertices[neigh]) + d = bottleneck[neigh] + weight = l / (d * d + b) + rows.append(tetra) + cols.append(neigh) + data.append(weight) + graph = csr_matrix((data, (rows, cols)),shape=(len(simplices), len(simplices))) + return graph - 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 + def dijkstra(self, cavity, graph, simplices, neighbors, vertices, points, vdw_radii): + + import numpy as np + 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)] + + for start_global in cavity.starting_tetrahedron: + if start_global not in global_to_local: + continue + start_local = global_to_local[start_global] + distances, predecessors = dijkstra(cavity_graph,directed=False,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) - b = 1e-3 - return l / (d**2 + b) - - def dijkstra_algorithm(start, goal, tetrahedra_set): - pq = [(0, start)] - distances = {start: 0} - previous = {start: None} + paths = {} + stack = [(start_local, [start_local])] - 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]: + while stack: + node, path = stack.pop() + paths[node] = path + + for child in parent_to_children.get(node, []): + stack.append((child, path + [child])) + + for exit_global in cavity.end_tetrahedra: + if exit_global == start_global: + continue + if exit_global not in global_to_local: + continue + exit_local = global_to_local[exit_global] + if np.isinf(distances[exit_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) + path_local = paths.get(exit_local) + if path_local is None: + continue + + path_global = cavity_tetra[path_local] + + node = exit_local + + channel = Channel(path_global,*self.process_channel( + path_global, + vertices, + points, + vdw_radii, + simplices)) + cavity.add_channel(channel) def calculate_max_radius(self, vertice, points, vdw_radii, simp): atom_positions = points[simp] From 4b82727714a718031c4612afd9ca6e1d1895f9a1 Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 13:11:08 -0400 Subject: [PATCH 10/16] added centers and radii to process_channel outputs to match changes to Channel class --- prody/proteins/channels.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 1b0fadfcb..b2b185d50 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -1633,7 +1633,7 @@ def process_channel(self, tetrahedra, voronoi_vertices, points, vdw_radii, simp) length = self.calculate_channel_length(centerline_spline) volume = self.calculate_channel_volume(centerline_spline, radius_spline) - return centerline_spline, radius_spline, length, bottleneck, volume + return centerline_spline, radius_spline, length, bottleneck, volume,centers, radii 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]) From e52a47e4ba5934298514358c31ed01957ea9ef6b Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 13:12:20 -0400 Subject: [PATCH 11/16] added calc_circumcenters function to ChannelCalculator class --- prody/proteins/channels.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index b2b185d50..9bb2a3090 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -1789,4 +1789,9 @@ def set_starting_tetrahedra_from_point(self, cavities, vertices, start_point): cavity.set_starting_tetrahedron(np.array([chosen])) + def calc_circumcenters(self, dela): + eq = dela.equations + scale = dela.paraboloid_scale + centers = -eq[:, :-2] / (2 * scale * eq[:, -2][:, None]) + return centers From 8893430e3e50eaa25f4c1eb0fe2ff20b7fb60c71 Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 13:46:56 -0400 Subject: [PATCH 12/16] added normalModeCavityAnalysis --- prody/proteins/channels.py | 144 +++++++++++++++++++++++++++++++++---- 1 file changed, 130 insertions(+), 14 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 9bb2a3090..eab09a488 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -16,6 +16,8 @@ from prody.utilities import importLA, checkCoords, showFigure, getCoords, isListLike from prody.measure import calcDistance, calcAngle, calcCenter from prody.measure.contacts import findNeighbors +from prody.dynamics.sampling import traverseMode +from prody.dynamics.anm import ANM from prody.proteins import writePDB, parsePDB, parsePQR from collections import Counter @@ -574,6 +576,7 @@ def calcChannels(atoms, r1=3, r2=1.25, min_depth=10, bottleneck=1, sparsity=15, return channels, [coords, s_srf.simp, merged_cavities, s_clr.simp] def _ag_worker(args): + """Process a single atomgroup for calcChannelsMultipleAtomsGroups""" atoms,kwargs,i,filename= args # LOGGER.info("Processing atom group " + str(filename)) start = time.perf_counter() @@ -588,6 +591,76 @@ def _ag_worker(args): def calcChannelsMultipleAtomGroups(atomgroups, **kwargs): + """"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. + + 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. + :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. + :type output_path: str or None + + :param filenames: A list of filenames for the output files. If None, the default filenames will be used. + :type filenames: list of 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. + :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 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. + :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: float + + :param bottleneck: The minimum allowed bottleneck size (narrowest point) for the channels. Default is 1. + :type bottleneck: 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 + + :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. + :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. + + Example usage: + channels, surface = calcChannels(atoms, output_path="channels", separate=True) + + 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) """ + from multiprocessing import Pool, cpu_count filenames = kwargs.pop('filenames',[None]*len(atomgroups)) max_proc = kwargs.pop('max_proc',None) @@ -615,13 +688,14 @@ def calcChannelsMultipleAtomGroups(atomgroups, **kwargs): return channels_all, surfaces_all, times_all def _frame_worker(args): - atoms, coords, kwargs, frame_id = args - + """Process a single frame for calcChannelsMultipleFrames""" + atoms, coords, kwargs, frame_id, filename = args + atoms_local = atoms.copy() atoms_local.setCoords(coords) try: - channels, surfaces = calcChannels(atoms_local, **kwargs) + channels, surfaces = calcChannels(atoms_local, filename=filename, **kwargs) if not channels: return frame_id,[],[] except: @@ -629,8 +703,8 @@ def _frame_worker(args): return frame_id, channels, surfaces def _model_worker(args): - """Process a single model from multi-model PDB""" - atoms, kwargs, frame_id = args + """Process a single model for calcChannelsMultipleFrames""" + atoms, kwargs, frame_id, filename = args start_frame = kwargs.pop(0,"start_frame") frame_idx = frame_id + start_frame LOGGER.info("Model: {0}".format(frame_idx)) @@ -639,7 +713,7 @@ def _model_worker(args): atoms_copy.setACSIndex(frame_idx) try: - channels, surfaces = calcChannels(atoms_copy, **kwargs) + channels, surfaces = calcChannels(atoms_copy, filename=filename, **kwargs) if not channels: return frame_id,[],[] except: @@ -669,6 +743,9 @@ def calcChannelsMultipleFrames(atoms, trajectory, **kwargs): If None, results are not saved. Default is None. :type output_path: str or None + :param filenames: A list of filenames for the output files. If None, the default filenames will be used. + :type filenames: list of 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. :type separate: bool @@ -727,7 +804,7 @@ def calcChannelsMultipleFrames(atoms, trajectory, **kwargs): max_proc = kwargs.pop('max_proc',None) if max_proc is None: max_proc = max(1, cpu_count()//2) - + start_frame = kwargs.pop('start_frame', 0) stop_frame = kwargs.pop('stop_frame', -1) @@ -739,11 +816,11 @@ def calcChannelsMultipleFrames(atoms, trajectory, **kwargs): traj = trajectory[start_frame:] else: traj = trajectory[start_frame:stop_frame+1] - + filenames = kwargs.pop('filenames',[None]*len(traj)) if max_proc == 1: - results = [_frame_worker((atoms,frame.getCoords(),kwargs,i)) for i, frame in enumerate(traj)] + results = [_frame_worker((atoms,frame.getCoords(),kwargs,i,filenames[i])) for i, frame in enumerate(traj)] else: - tasks = ((atoms,frame.getCoords(),kwargs,i) for i, frame in enumerate(trajectory)) + tasks = ((atoms,frame.getCoords(),kwargs,i,filenames[i]) for i, frame in enumerate(trajectory)) with Pool(processes=max_proc) as pool: results = pool.map(_frame_worker, tasks,chunksize=builtins.max(1, len(traj)//(max_proc*4))) @@ -760,11 +837,11 @@ def calcChannelsMultipleFrames(atoms, trajectory, **kwargs): num_models = len(atoms.getCoordsets()[start_frame:]) else: num_models = len(atoms.getCoordsets()[start_frame:stop_frame+1]) - + filenames = kwargs.pop('filenames',[None]*len(num_models)) if max_proc == 1: - results = [_model_worker((atoms,kwargs,i)) for i in range(num_models)] + results = [_model_worker((atoms,kwargs,i,filenames[i])) for i in range(num_models)] else: - tasks = ((atoms,kwargs,i) for i in range(num_models)) + tasks = ((atoms,kwargs,i,filenames[i]) for i in range(num_models)) with Pool(processes=max_proc) as pool: results = pool.map(_model_worker, tasks,chunksize=builtins.max(1, len(traj)//(max_proc*4))) @@ -1248,6 +1325,46 @@ def write_merge_surf_pdb(merged_surface, filename, nr_pdbs): merged_surface = merge_surfaces(surfaces) write_merge_surf_pdb(merged_surface, output_file_name, nr_pdbs) +def normalModeCavityAnalysis(structure,num_modes=20,scale=2,n_steps=20,**kwargs): + + if isinstance(structure, str) or hasattr(structure, '__fspath__'): + ag = parsePDB(str(structure)) + elif isinstance(structure, AtomGroup): + ag = structure.copy() + else: + raise ValueError("structure must be a file path or AtomGroup") + title = kwargs.pop('title',ag.getTitle()) + if (nma := kwargs.pop('nma', None)): + if not isinstance(nma, NMA): + raise ValueError("nma must be an instance of prody.NMA") + else: + nma = ANM() + if (selection := kwargs.pop('selection', 'calpha')): + ag_cg = ag.select(selection) + nma.buildHessian(ag) + nma.calcModes(n_modes=num_modes) + if (extend := kwargs.pop('extend', 'all')): + if extend == 'all': + ext_nma, ext_ag = nma.extendModel(nma,ag_cg,ag) + elif extend: + ext_nma, ext_ag = nma.extendModel(nma,ag_cg,ag.select(extend)) + else: + ext_nma, ext_ag = nma, ag_cg + + agc = ext_ag.copy() + calcMap=[f"{ag.getTitle()}"] + for mode in range(num_modes): + trav_ens = traverseMode(ext_nma[mode],ext_ag,n_steps,rmsd=scale) + agc.addCoordset(trav_ens[:n_steps//2]) + agc.addCoordset(trav_ens[(n_steps//2)+1:]) + for n in range(n_steps*2+1): + if n == n_steps: + continue + else: + output = f"{title}_mode{mode+1}_step{n-n_steps}" + calcMap.append(f"{output}") + channels, surfaces = calcChannelsMultipleFrames(agc,output_path=output,separate=kwargs.get('separate', False),filenames=calcMap) + return channels,surfaces class Channel: def __init__(self, tetrahedra, centerline_spline, radius_spline, length, bottleneck, volume, centers, radii): @@ -1522,7 +1639,6 @@ def find_deepest_tetrahedra(self, cavities, neighbors): cavity.set_depth(max_depth) def build_sparse_graph(self, simplices, neighbors, vertices, points, vdw_radii): - import numpy as np from scipy.sparse import csr_matrix From 0ba87b1ba950ed5262fef1512a2f9b3d750c210e Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 13:53:44 -0400 Subject: [PATCH 13/16] removed normalModeCavityAnalysis, needed updating --- prody/proteins/channels.py | 42 ++------------------------------------ 1 file changed, 2 insertions(+), 40 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index eab09a488..638622764 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -13,6 +13,7 @@ from prody import LOGGER, SETTINGS, PY3K from prody.atomic import AtomGroup, Atom, Atomic, Selection, Select from prody.atomic import flags, sliceAtomicData +from prody.dynamics.nma import NMA from prody.utilities import importLA, checkCoords, showFigure, getCoords, isListLike from prody.measure import calcDistance, calcAngle, calcCenter from prody.measure.contacts import findNeighbors @@ -1325,46 +1326,7 @@ def write_merge_surf_pdb(merged_surface, filename, nr_pdbs): merged_surface = merge_surfaces(surfaces) write_merge_surf_pdb(merged_surface, output_file_name, nr_pdbs) -def normalModeCavityAnalysis(structure,num_modes=20,scale=2,n_steps=20,**kwargs): - - if isinstance(structure, str) or hasattr(structure, '__fspath__'): - ag = parsePDB(str(structure)) - elif isinstance(structure, AtomGroup): - ag = structure.copy() - else: - raise ValueError("structure must be a file path or AtomGroup") - title = kwargs.pop('title',ag.getTitle()) - if (nma := kwargs.pop('nma', None)): - if not isinstance(nma, NMA): - raise ValueError("nma must be an instance of prody.NMA") - else: - nma = ANM() - if (selection := kwargs.pop('selection', 'calpha')): - ag_cg = ag.select(selection) - nma.buildHessian(ag) - nma.calcModes(n_modes=num_modes) - if (extend := kwargs.pop('extend', 'all')): - if extend == 'all': - ext_nma, ext_ag = nma.extendModel(nma,ag_cg,ag) - elif extend: - ext_nma, ext_ag = nma.extendModel(nma,ag_cg,ag.select(extend)) - else: - ext_nma, ext_ag = nma, ag_cg - - agc = ext_ag.copy() - calcMap=[f"{ag.getTitle()}"] - for mode in range(num_modes): - trav_ens = traverseMode(ext_nma[mode],ext_ag,n_steps,rmsd=scale) - agc.addCoordset(trav_ens[:n_steps//2]) - agc.addCoordset(trav_ens[(n_steps//2)+1:]) - for n in range(n_steps*2+1): - if n == n_steps: - continue - else: - output = f"{title}_mode{mode+1}_step{n-n_steps}" - calcMap.append(f"{output}") - channels, surfaces = calcChannelsMultipleFrames(agc,output_path=output,separate=kwargs.get('separate', False),filenames=calcMap) - return channels,surfaces + class Channel: def __init__(self, tetrahedra, centerline_spline, radius_spline, length, bottleneck, volume, centers, radii): From 6632c29f517dd2692945423d7c56a596062d5e6e Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 26 Jun 2026 15:23:23 -0400 Subject: [PATCH 14/16] added additional error handling for multiprocessing --- prody/proteins/channels.py | 84 ++++++++++++++++++++++++++------------ 1 file changed, 58 insertions(+), 26 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 638622764..ee065bf32 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -580,17 +580,25 @@ def _ag_worker(args): """Process a single atomgroup for calcChannelsMultipleAtomsGroups""" atoms,kwargs,i,filename= args # LOGGER.info("Processing atom group " + str(filename)) - start = time.perf_counter() try: channels, surfaces = calcChannels(atoms,filename = filename, **kwargs) if not channels: return i,[],[],0 - total = time.perf_counter() - start - return i, channels, surfaces, total + return i, channels, surfaces except: return i, [], [], 0 - +def _run_with_pool(_worker,tasks, max_proc,chunksize,start_method=None): + from multiprocessing import Pool + if start_method: + ctx = multiprocessing.get_context(start_method) + PoolClass = ctx.Pool + else: + PoolClass = Pool + + with PoolClass(processes=max_proc) as pool: + return pool.map(_worker, tasks,chunksize=chunksize) + def calcChannelsMultipleAtomGroups(atomgroups, **kwargs): """"Computes and identifies channels within a molecular structure using Voronoi and Delaunay tessellations. @@ -666,14 +674,24 @@ def calcChannelsMultipleAtomGroups(atomgroups, **kwargs): filenames = kwargs.pop('filenames',[None]*len(atomgroups)) max_proc = kwargs.pop('max_proc',None) if max_proc is None: - max_proc = max(1, cpu_count()//2) + max_proc = builtins.max(1, cpu_count()//2) if max_proc == 1: results = [_ag_worker((ag,kwargs,i,filenames[i])) for i, ag in enumerate(atomgroups)] else: tasks = ((ag,kwargs,i,filenames[i]) for i, ag in enumerate(atomgroups)) - with Pool(processes=max_proc) as pool: - results = pool.map(_ag_worker, tasks,chunksize=builtins.max(1, len(atomgroups)//(max_proc*4))) + chunksize = builtins.max(1, len(atomgroups)//(max_proc*4)) + try: + results = _run_with_pool(_ag_worker,tasks,max_proc,chunksize) + except (OSError,EOFError): + try: + results = _run_with_pool(_ag_worker,tasks,max_proc,chunksize,'fork') + except (OSError, EOFError): + try: + results = _run_with_pool(_ag_worker,tasks,max_proc,chunksize,'spawn') + except Exception as e: + raise RuntimeError("Multiprocessing failed with all start methods. " + "Try running with max_proc=1.") from e results.sort(key=lambda x: x[0]) @@ -682,11 +700,10 @@ def calcChannelsMultipleAtomGroups(atomgroups, **kwargs): for channel in channels: channel.build_splines() surfaces_all = [r[2] for r in results] - times_all = [r[3] for r in results] failed = [filenames[r[0]] for r in results if (not r[1] and not r[2])] if failed: LOGGER.warning(f"WARNING: {len(failed)} proteins failed or No Channels Detected: {', '.join(str(f) for f in failed)}") - return channels_all, surfaces_all, times_all + return channels_all, surfaces_all def _frame_worker(args): """Process a single frame for calcChannelsMultipleFrames""" @@ -706,7 +723,7 @@ def _frame_worker(args): def _model_worker(args): """Process a single model for calcChannelsMultipleFrames""" atoms, kwargs, frame_id, filename = args - start_frame = kwargs.pop(0,"start_frame") + start_frame = kwargs.pop("start_frame",0) frame_idx = frame_id + start_frame LOGGER.info("Model: {0}".format(frame_idx)) # Make a copy to avoid state conflicts @@ -722,7 +739,7 @@ def _model_worker(args): return frame_id, channels, surfaces -def calcChannelsMultipleFrames(atoms, trajectory, **kwargs): +def calcChannelsMultipleFrames(atoms, trajectory=None, **kwargs): """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 @@ -804,7 +821,7 @@ def calcChannelsMultipleFrames(atoms, trajectory, **kwargs): max_proc = kwargs.pop('max_proc',None) if max_proc is None: - max_proc = max(1, cpu_count()//2) + max_proc = builtins.max(1, cpu_count()//2) start_frame = kwargs.pop('start_frame', 0) stop_frame = kwargs.pop('stop_frame', -1) @@ -812,7 +829,7 @@ def calcChannelsMultipleFrames(atoms, trajectory, **kwargs): if trajectory is not None: if isinstance(trajectory, Atomic): trajectory = Ensemble(trajectory) - + traj=None if stop_frame == -1: traj = trajectory[start_frame:] else: @@ -821,39 +838,56 @@ def calcChannelsMultipleFrames(atoms, trajectory, **kwargs): if max_proc == 1: results = [_frame_worker((atoms,frame.getCoords(),kwargs,i,filenames[i])) for i, frame in enumerate(traj)] else: - tasks = ((atoms,frame.getCoords(),kwargs,i,filenames[i]) for i, frame in enumerate(trajectory)) - with Pool(processes=max_proc) as pool: - results = pool.map(_frame_worker, tasks,chunksize=builtins.max(1, len(traj)//(max_proc*4))) + tasks = ((atoms,frame.getCoords(),kwargs,i,filenames[i]) for i, frame in enumerate(traj)) + chunksize = builtins.max(1, len(traj)//(max_proc*4)) + try: + results = _run_with_pool(_frame_worker,tasks,max_proc,chunksize) + except (OSError,EOFError): + try: + results = _run_with_pool(_frame_worker,tasks,max_proc,chunksize,'fork') + except (OSError, EOFError): + try: + results = _run_with_pool(_frame_worker,tasks,max_proc,chunksize,'spawn') + except Exception as e: + raise RuntimeError("Multiprocessing failed with all start methods. " + "Try running with max_proc=1.") from e results.sort(key=lambda x: x[0]) channels_all = [r[1] for r in results] surfaces_all = [r[2] for r in results] - times_all = [r[3] for r in results] - return channels_all, surfaces_all, times_all + return channels_all, surfaces_all else: if atoms.numCoordsets() > 1: if stop_frame == -1: num_models = len(atoms.getCoordsets()[start_frame:]) else: num_models = len(atoms.getCoordsets()[start_frame:stop_frame+1]) - filenames = kwargs.pop('filenames',[None]*len(num_models)) + filenames = kwargs.pop('filenames',[None]*num_models) if max_proc == 1: results = [_model_worker((atoms,kwargs,i,filenames[i])) for i in range(num_models)] else: tasks = ((atoms,kwargs,i,filenames[i]) for i in range(num_models)) - - with Pool(processes=max_proc) as pool: - results = pool.map(_model_worker, tasks,chunksize=builtins.max(1, len(traj)//(max_proc*4))) + chunksize = builtins.max(1, num_models//(max_proc*4)) + try: + results = _run_with_pool(_model_worker,tasks,max_proc,chunksize) + except (OSError,EOFError): + try: + results = _run_with_pool(_model_worker,tasks,max_proc,chunksize,'fork') + except (OSError, EOFError): + try: + results = _run_with_pool(_model_worker,tasks,max_proc,chunksize,'spawn') + except Exception as e: + raise RuntimeError("Multiprocessing failed with all start methods. " + "Try running with max_proc=1.") from e results.sort(key=lambda x: x[0]) channels_all = [r[1] for r in results] surfaces_all = [r[2] for r in results] - times_all = [r[3] for r in results] - return channels_all, surfaces_all, times_all + return channels_all, surfaces_all else: LOGGER.info("Include trajectory or use multi-model PDB file.") @@ -1325,8 +1359,6 @@ def write_merge_surf_pdb(merged_surface, filename, nr_pdbs): nr_pdbs = nr_pdbs+1 merged_surface = merge_surfaces(surfaces) write_merge_surf_pdb(merged_surface, output_file_name, nr_pdbs) - - class Channel: def __init__(self, tetrahedra, centerline_spline, radius_spline, length, bottleneck, volume, centers, radii): From e3618d57020d8fd7beec2e8ca105945d2e94fedc Mon Sep 17 00:00:00 2001 From: MatthewLicht Date: Sat, 4 Jul 2026 21:12:39 -0400 Subject: [PATCH 15/16] removed unnecessary import and altered formatting for legibility --- prody/proteins/channels.py | 96 ++++++++++++++++++++++++++++++-------- 1 file changed, 77 insertions(+), 19 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index ee065bf32..2e6976da3 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -13,7 +13,6 @@ from prody import LOGGER, SETTINGS, PY3K from prody.atomic import AtomGroup, Atom, Atomic, Selection, Select from prody.atomic import flags, sliceAtomicData -from prody.dynamics.nma import NMA from prody.utilities import importLA, checkCoords, showFigure, getCoords, isListLike from prody.measure import calcDistance, calcAngle, calcCenter from prody.measure.contacts import findNeighbors @@ -677,18 +676,37 @@ def calcChannelsMultipleAtomGroups(atomgroups, **kwargs): max_proc = builtins.max(1, cpu_count()//2) if max_proc == 1: - results = [_ag_worker((ag,kwargs,i,filenames[i])) for i, ag in enumerate(atomgroups)] + results = [_ag_worker((ag,kwargs,i,filenames[i])) + for i, ag in enumerate(atomgroups)] else: - tasks = ((ag,kwargs,i,filenames[i]) for i, ag in enumerate(atomgroups)) + tasks = ((ag,kwargs,i,filenames[i]) + for i, ag in enumerate(atomgroups)) chunksize = builtins.max(1, len(atomgroups)//(max_proc*4)) try: - results = _run_with_pool(_ag_worker,tasks,max_proc,chunksize) + results = _run_with_pool( + _ag_worker, + tasks + ,max_proc, + chunksize + ) except (OSError,EOFError): try: - results = _run_with_pool(_ag_worker,tasks,max_proc,chunksize,'fork') + results = _run_with_pool( + _ag_worker, + tasks, + max_proc, + chunksize, + 'fork' + ) except (OSError, EOFError): try: - results = _run_with_pool(_ag_worker,tasks,max_proc,chunksize,'spawn') + results = _run_with_pool( + _ag_worker, + tasks, + max_proc, + chunksize, + 'spawn' + ) except Exception as e: raise RuntimeError("Multiprocessing failed with all start methods. " "Try running with max_proc=1.") from e @@ -700,7 +718,8 @@ def calcChannelsMultipleAtomGroups(atomgroups, **kwargs): for channel in channels: channel.build_splines() surfaces_all = [r[2] for r in results] - failed = [filenames[r[0]] for r in results if (not r[1] and not r[2])] + failed = [filenames[r[0]] + for r in results if (not r[1] and not r[2])] if failed: LOGGER.warning(f"WARNING: {len(failed)} proteins failed or No Channels Detected: {', '.join(str(f) for f in failed)}") return channels_all, surfaces_all @@ -812,12 +831,14 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, **kwargs): from multiprocessing import Pool, cpu_count try: - coords = getCoords(atoms) + coords = (coords._getCoords() if hasattr(coords, '_getCoords') else + coords.getCoords()) except AttributeError: try: checkCoords(coords) except TypeError: - raise TypeError('coords must be an object with `getCoords` method') + raise TypeError('coords must be a Numpy array or an object ' + 'with `getCoords` method') max_proc = kwargs.pop('max_proc',None) if max_proc is None: @@ -836,9 +857,24 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, **kwargs): traj = trajectory[start_frame:stop_frame+1] filenames = kwargs.pop('filenames',[None]*len(traj)) if max_proc == 1: - results = [_frame_worker((atoms,frame.getCoords(),kwargs,i,filenames[i])) for i, frame in enumerate(traj)] + results = [_frame_worker( + (atoms, + frame.getCoords(), + kwargs, + i, + filenames[i]) + ) for i, frame in enumerate(traj)] + else: - tasks = ((atoms,frame.getCoords(),kwargs,i,filenames[i]) for i, frame in enumerate(traj)) + tasks = ( + (atoms, + frame.getCoords(), + kwargs, + i, + filenames[i]) + for i, frame in enumerate(traj) + ) + chunksize = builtins.max(1, len(traj)//(max_proc*4)) try: results = _run_with_pool(_frame_worker,tasks,max_proc,chunksize) @@ -866,18 +902,37 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, **kwargs): num_models = len(atoms.getCoordsets()[start_frame:stop_frame+1]) filenames = kwargs.pop('filenames',[None]*num_models) if max_proc == 1: - results = [_model_worker((atoms,kwargs,i,filenames[i])) for i in range(num_models)] + results = [_model_worker((atoms,kwargs,i,filenames[i])) + for i in range(num_models)] else: - tasks = ((atoms,kwargs,i,filenames[i]) for i in range(num_models)) + tasks = ((atoms,kwargs,i,filenames[i]) + for i in range(num_models)) chunksize = builtins.max(1, num_models//(max_proc*4)) try: - results = _run_with_pool(_model_worker,tasks,max_proc,chunksize) + results = _run_with_pool( + _model_worker, + tasks, + max_proc, + chunksize + ) except (OSError,EOFError): try: - results = _run_with_pool(_model_worker,tasks,max_proc,chunksize,'fork') + results = _run_with_pool( + _model_worker, + tasks, + max_proc, + chunksize, + 'fork' + ) except (OSError, EOFError): try: - results = _run_with_pool(_model_worker,tasks,max_proc,chunksize,'spawn') + results = _run_with_pool( + _model_worker, + tasks, + max_proc, + chunksize, + 'spawn' + ) except Exception as e: raise RuntimeError("Multiprocessing failed with all start methods. " "Try running with max_proc=1.") from e @@ -959,7 +1014,7 @@ def getChannelParameters(channels, **kwargs): lengths, bottlenecks, volumes = frame LOGGER.info("Frame {0}".format(frame_nr)) 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))) + 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 multi_model_param @@ -1573,7 +1628,9 @@ def dfs(tetra_index): if not visited[index]: visited[index] = True current_group.append(index) - stack.extend(neighbor for neighbor in neigh[index] if neighbor != -1 and not visited[neighbor]) + stack.extend(neighbor for + neighbor in neigh[index] + if neighbor != -1 and not visited[neighbor]) return np.array(current_group) for i in range(x): @@ -1604,7 +1661,8 @@ def get_surface_cavities(self, cavities, interior_simplices, second_layer, state def merge_cavities(self, cavities, simplices): - merged_tetrahedra = np.concatenate([cavity.tetrahedra for cavity in cavities]) + merged_tetrahedra = np.concatenate([cavity.tetrahedra + for cavity in cavities]) return simplices[merged_tetrahedra] def find_deepest_tetrahedra(self, cavities, neighbors): From d9b95d0acdcf73ec0ad2fae9ba9f01fda22c5e54 Mon Sep 17 00:00:00 2001 From: MatthewLicht Date: Mon, 6 Jul 2026 15:31:22 -0400 Subject: [PATCH 16/16] resolved formatting issues --- prody/proteins/channels.py | 75 ++++++++++++++++++++++++-------------- 1 file changed, 48 insertions(+), 27 deletions(-) diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 2e6976da3..a66f15da7 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -877,13 +877,29 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, **kwargs): chunksize = builtins.max(1, len(traj)//(max_proc*4)) try: - results = _run_with_pool(_frame_worker,tasks,max_proc,chunksize) + results = _run_with_pool( + _frame_worker, + tasks, + max_proc, + chunksize + ) except (OSError,EOFError): try: - results = _run_with_pool(_frame_worker,tasks,max_proc,chunksize,'fork') + results = _run_with_pool( + _frame_worker, + tasks, + max_proc, + chunksize, + 'fork') except (OSError, EOFError): try: - results = _run_with_pool(_frame_worker,tasks,max_proc,chunksize,'spawn') + results = _run_with_pool( + _frame_worker, + tasks, + max_proc, + chunksize, + 'spawn' + ) except Exception as e: raise RuntimeError("Multiprocessing failed with all start methods. " "Try running with max_proc=1.") from e @@ -902,37 +918,39 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, **kwargs): num_models = len(atoms.getCoordsets()[start_frame:stop_frame+1]) filenames = kwargs.pop('filenames',[None]*num_models) if max_proc == 1: - results = [_model_worker((atoms,kwargs,i,filenames[i])) - for i in range(num_models)] + results = [ + _model_worker((atoms,kwargs,i,filenames[i])) + for i in range(num_models) + ] else: tasks = ((atoms,kwargs,i,filenames[i]) for i in range(num_models)) chunksize = builtins.max(1, num_models//(max_proc*4)) try: results = _run_with_pool( - _model_worker, - tasks, - max_proc, - chunksize - ) + _model_worker, + tasks, + max_proc, + chunksize + ) except (OSError,EOFError): try: results = _run_with_pool( - _model_worker, - tasks, - max_proc, - chunksize, - 'fork' - ) + _model_worker, + tasks, + max_proc, + chunksize, + 'fork' + ) except (OSError, EOFError): try: results = _run_with_pool( - _model_worker, - tasks, - max_proc, - chunksize, - 'spawn' - ) + _model_worker, + tasks, + max_proc, + chunksize, + 'spawn' + ) except Exception as e: raise RuntimeError("Multiprocessing failed with all start methods. " "Try running with max_proc=1.") from e @@ -941,6 +959,9 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, **kwargs): channels_all = [r[1] for r in results] surfaces_all = [r[2] for r in results] + for channels in channels_all: + for channel in channels: + channel.build_splines() return channels_all, surfaces_all @@ -1014,7 +1035,7 @@ def getChannelParameters(channels, **kwargs): lengths, bottlenecks, volumes = frame LOGGER.info("Frame {0}".format(frame_nr)) 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))) + 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 multi_model_param @@ -1131,8 +1152,8 @@ def getChannelResidueNames(atoms, channels, **kwargs): try: checkCoords(coords) except TypeError: - raise TypeError('coords must be an object ' - 'with `getCoords` method') + raise TypeError('coords must be a Numpy array or an object ' + 'with `getCoords` method') distA = kwargs.pop('distA', 4) residues_file_name = kwargs.pop('residues_file_name', None) @@ -1231,8 +1252,8 @@ def selectChannelBySelection(atoms, residue_sele, **kwargs): try: checkCoords(coords) except TypeError: - raise TypeError('coords must be an object ' - 'with `getCoords` method') + raise TypeError('coords must be a Numpy array or an object ' + 'with `getCoords` method') import os, shutil import numpy as np