diff --git a/prody/proteins/channels.py b/prody/proteins/channels.py index 60b95cdba..993b0d33e 100644 --- a/prody/proteins/channels.py +++ b/prody/proteins/channels.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- """This module is called CaviTracer and defines functions for calculating -channels, tunnels, and surface cavities within protein structure. +channels, tunnels, pores, and surface cavities within protein structure. """ __author__ = 'Karolina Mikulska-Ruminska', 'Jan Brezovsky', 'Eryk Trzcinski' @@ -31,7 +31,10 @@ 'getSurfaceCavityResidueNamesMultipleFrames', 'getSurfaceCavityParametersMultipleFrames', 'getChannelParametersMultipleFrames', '_reportAtomsInputComposition', - 'getChannelResidueNamesMultipleFrames', 'scanChannelParameters'] + 'getChannelResidueNamesMultipleFrames', 'calcPoresFromChannels', + 'showPores', 'getPoreParameters', 'getPoreResidueNames', + 'calcPoresFromChannelsMultipleFrames', 'getPoreParametersMultipleFrames', + 'getPoreResidueNamesMultipleFrames', 'scanChannelParameters'] # Sampling of the enclosure test used to strip the moat (see # ChannelCalculator.calcEnclosure). These are constants, not knobs: the enclosure @@ -160,6 +163,27 @@ def _surfaceFromPqrWorker(args): return surface +def _calcChannelsMultipleFramesWorker(args): + """Compute channels. Supporting function for muliprocessing in :func:`calcChannelsMultipleFrames`.""" + frame_nr, atoms, frame_coords, frame_output_path, separate, start_point, return_details, kwargs = args + + LOGGER.info("Frame/model: {0}".format(frame_nr)) + atoms_copy = atoms.copy() + atoms_copy.setCoords(frame_coords) + + return calcChannels(atoms_copy, output_path=frame_output_path, separate=separate, + start_point=start_point, return_details=return_details, **kwargs) + + +def _calcPoresFromChannelsWorker(args): + """Reconstruct pores from channels. Supporting function for multiprocessing + in :func:`calcPoresFromChannelsMultipleFrames`.""" + frame_nr, channels, details, output_path, separate, kwargs = args + LOGGER.info("Frame/model: {0}".format(frame_nr)) + return calcPoresFromChannels(channels, details, output_path=output_path, + separate=separate, **kwargs) + + def _reportAtomsInputComposition(atoms): """Report the composition of atoms supplied for channel analysis. @@ -366,7 +390,7 @@ def getVmdModel(vmd_path, atoms, representation='NewCartoon'): def showChannels(channels, model=None, surface=None): - """Visualizes the channels, and optionally, the molecular model and + """Visualizes the channels or pores, and optionally, the molecular model and surface, using Open3D. This function renders a 3D visualization of molecular channels based on @@ -486,6 +510,8 @@ def create_mesh_from_spline(centerline_spline, radius_spline, n=5): else: LOGGER.info("Nothing to visualize.") +showPores = showChannels + def showCavities(surface, show_surface=False): """Visualizes the cavities within a molecular surface using Open3D. @@ -823,7 +849,8 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, 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): + weighted_cache=True, weighted_mouth_depth=2.5, edge_cost=None, + return_details=False): """Computes and identifies channels within a molecular structure using Voronoi and Delaunay tessellations. @@ -1100,6 +1127,12 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, unaffected by this choice. :type edge_cost: str or None + :arg return_details: If True return an additional dictionary containing + internal calculation data, including the channel calculator, simplices, + neighboring tetrahedra, Voronoi vertices, atomic coordinates, and van der + Waals radii. Default is False. + :type return_details: bool + :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. @@ -1401,7 +1434,7 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, c_filtered_cavities, min_volume, max_volume) merged_cavities = calculator.mergeCavities(c_filtered_cavities, s_clr.simp) - + # Early-return for the calcSurfaceCavities function: if cavities_only: LOGGER.info("Returning surface cavities") @@ -1424,7 +1457,7 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, 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), @@ -1472,11 +1505,227 @@ def calcChannels(atoms, output_path=None, separate=False, start_point=None, LOGGER.info("No output path given.") LOGGER.report('Channel calculation completed in %.2fs.', '_prody_calcChannels') + + # Additional information can be obtained + if return_details: + details = {'calculator': calculator, + 'simplices': s_clr.simp, + 'neighbors': s_clr.neigh, + 'vertices': s_clr.verti, + 'coords': coords, + 'vdw_radii': vdw_radii} + + return channels, [coords, s_srf.simp, merged_cavities, s_clr.simp], details + return channels, [coords, s_srf.simp, merged_cavities, s_clr.simp] + +def calcPoresFromChannels(channels, details, min_end_to_end=None, max_end_to_end=None, + min_bottleneck=None, max_bottleneck=None, min_length=None, max_length=None, + min_volume=None, max_volume=None, output_path=None, separate=False): + """Construct potential pores from previously identified channels using + :func:`calcChannels`. This function performs a post-processing analysis of + channels and requires ``return_details`` set to ``True`` in :func:`calcChannels`. + The `separate` parameter controls whether each pore is additionally saved to a + separate file. + + The pore-construction procedure consists of the following steps: + + 1. Group channels according to their starting tetrahedron. + 2. Generate all unique pairs of channels within each group. + 3. Identify the common initial segment and the last tetrahedron shared by + each pair of channel paths. + 4. Join the non-overlapping parts of the two channels at their branching + tetrahedron to obtain a surface-to-surface path. + 5. Reject paths containing loops or discontinuities between neighboring + tetrahedra. + 6. Remove identical paths and paths differing only in direction. + 7. Recalculate the centerline spline, radius profile, length, bottleneck, + and volume of each resulting pore using approach implemented for channels + identification and visualization. + 8. Pores are filtered based on the given criteria (``min_end_to_end``, + ``max_end_to_end``, ``min_bottleneck``, ``max_bottleneck``, ``min_length``, + ``max_length``, ``min_volume``, ``max_volume``). + + :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 details: Additional calculation data returned by + :func:`calcChannels` with ``return_details=True``. The dictionary must + contain ``calculator``, ``simplices``, ``neighbors``, ``vertices``, + ``coords``, and ``vdw_radii``. + :type details: dict + + :arg min_end_to_end: Minimum allowed distance between the two pore + openings. Pores with a smaller end-to-end distance will be excluded. + Default is None. + :type min_end_to_end: int, float + + :arg max_end_to_end: Maximum allowed distance between the two pore + openings. Pores with a larger end-to-end distance will be excluded. + Default is None. + :type max_end_to_end: int, float + + :arg min_bottleneck: Minimum allowed bottleneck radius of the pore. + Pores with a smaller bottleneck will be excluded. + Default is None. + :type min_bottleneck: int, float + + :arg max_bottleneck: Maximum allowed bottleneck radius of the pore. + Pores with a larger bottleneck will be excluded. + Default is None. + :type max_bottleneck: int, float + + :arg min_length: Minimum allowed length of the pore. Pores shorter than this + value will be excluded. Default is None. + :type min_length: int, float + + :arg max_length: Maximum allowed length of the pore. Pores longer than this + value will be excluded. Default is None. + :type max_length: int, float + + :arg min_volume: Minimum allowed volume of the pore. Pores with a smaller volume + will be excluded. The value is given in cubic Angstroms. Default is None. + :type min_volume: int, float + + :arg max_volume: Maximum allowed volume of the pore. Pores with a larger volume + will be excluded. The value is given in cubic Angstroms. Default is None. + :type max_volume: int, float + + :arg output_path: Optional path to save the resulting pores and + associated data in PQR (or PDB) format. If None, results are not saved. + Default is None. + :type output_path: str or None + + :returns: Potential pores constructed from compatible channel pairs. + :rtype: list of Channel + + Usage: + channels, surface, details = calcChannels(protein, return_details=True) + pores = calcPoresFromChannels(channels, details, output_path='pores', separate=True) + """ + + if PY3K: + from pathlib import Path + else: + from pathlib2 import Path + + calculator = details['calculator'] + simplices = details['simplices'] + vertices = details['vertices'] + coords = details['coords'] + vdw_radii = details['vdw_radii'] + neighbors = details['neighbors'] + + pores = [] + pore_paths = [] + seen_paths = set() + channel_groups = {} + + # Group channels using their starting tetrahedron + for channel_index, channel in enumerate(channels): + path = np.asarray(channel.tetrahedra, dtype=np.intp) + # If channel is smaller than two tetrahedra (probably very rare) + # Those channels should be excluded because they can not be connected with others + if len(path) < 2: + continue + + start_tetrahedron = int(path[0]) + channel_groups.setdefault(start_tetrahedron, []).append((channel_index, channel, path)) + + from itertools import combinations + # Generate all channel pairs within each group + for start_tetrahedron, group in channel_groups.items(): + if len(group) < 2: + continue + + for (channel1_index, channel1, path1), (channel2_index, channel2, path2) in combinations(group, 2): + common_length = 0 + + for tetrahedron1, tetrahedron2 in zip(path1, path2): + if tetrahedron1 != tetrahedron2: + break + common_length += 1 + + if common_length == 0: + continue + + # If we have for example: path1: start → A → B → C → mouth 1 and path2: start → A → B → D → mouth 2 + # it will create mouth 1 → C → B → D → mouth 2 + branch_index = common_length - 1 + pore_path = np.concatenate((path1[branch_index:][::-1], path2[branch_index + 1:])) + # Reject paths containing loops + if len(np.unique(pore_path)) != len(pore_path): + continue + + # Continulity check of the pores (neighbours) + is_continuous = True + for tetrahedron1, tetrahedron2 in zip(pore_path[:-1], pore_path[1:]): + if tetrahedron2 not in neighbors[tetrahedron1]: + is_continuous = False + break + if not is_continuous: + continue + + # Remove identical paths + path_key = tuple(int(tetrahedron) for tetrahedron in pore_path) + canonical_key = min(path_key, path_key[::-1]) + + if canonical_key in seen_paths: + continue + + seen_paths.add(canonical_key) + pore_paths.append(pore_path) + + # Pores reconstruction + for pore_path in pore_paths: + # Filters - Distance between two ends of the pore + end_to_end = np.linalg.norm(vertices[pore_path[0]] - vertices[pore_path[-1]]) + + if min_end_to_end is not None and end_to_end < min_end_to_end: + continue + if max_end_to_end is not None and end_to_end > max_end_to_end: + continue + + centerline_spline, radius_spline, length, bottleneck, volume = calculator.processChannel( + pore_path, vertices, coords, vdw_radii, simplices) + + # Filters - bottleneck, length, volume + if min_bottleneck is not None and bottleneck < min_bottleneck: + continue + if max_bottleneck is not None and bottleneck > max_bottleneck: + continue + + if min_length is not None and length < min_length: + continue + if max_length is not None and length > max_length: + continue + + if min_volume is not None and volume < min_volume: + continue + if max_volume is not None and volume > max_volume: + continue + + pore = Channel(pore_path, centerline_spline, radius_spline, length, bottleneck, volume, 0.0) + pores.append(pore) + + if output_path: + output_path = Path(output_path) + if output_path.is_dir(): + output_path = output_path / "pores.pqr" + elif output_path.suffix not in (".pdb", ".pqr"): + output_path = output_path.with_suffix(".pqr") + + calculator.saveChannelsToPdb(pores, output_path, separate=separate) + + return pores + + def calcChannelsMultipleFrames(atoms, trajectory=None, output_path=None, - separate=False, start_point=None, **kwargs): + separate=False, start_point=None, max_proc=2, mp_context=None, **kwargs): """Compute channels for each frame in a given trajectory or multi-model PDB file. @@ -1511,6 +1760,21 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, output_path=None, seed selection based on the deepest tetrahedron). Coordinates must be given in Å. :type start_point: list, tuple, or ndarray (length 3), or None + :arg max_proc: Maximum number of parallel processes used for calculation. + If 1, files are processed serially. If None, all available CPU + cores are used. Default is 2. + :type max_proc: int or None + + :arg mp_context: Multiprocessing start method used for parallel pore + calculations. If `None`, the default method for the operating system + is used. Windows and macOS use the ``'spawn'`` method by default, + whereas Linux typically uses ``'fork'``. Setting + ``mp_context='spawn'`` can be potentially used on Linux, but might + be slower. Available values may include ``'spawn'``, ``'fork'``, + and ``'forkserver'``, depending on the operating system. + Default is `None`. + :type mp_context: str or None + :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. @@ -1529,7 +1793,6 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, output_path=None, channels_all, surfaces_all = calcChannelsMultipleFrames(atoms, trajectory=traj, output_path="channels.pdb", separate=False, start_point=[-10.353, -0.133, 5.608]) """ - if PY3K: if not checkAndImport('pathlib'): @@ -1554,6 +1817,10 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, output_path=None, channels_all = [] surfaces_all = [] + details_all = [] + tasks = [] + + return_details = kwargs.pop('return_details', False) start_frame = kwargs.pop('start_frame', 0) stop_frame = kwargs.pop('stop_frame', -1) @@ -1573,33 +1840,72 @@ def calcChannelsMultipleFrames(atoms, trajectory=None, output_path=None, 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) + frame_output_path = str(output_path) + "{0}.pqr".format(j0) else: - channels, surfaces = calcChannels(atoms_copy, start_point=start_point, **kwargs) - channels_all.append(channels) - surfaces_all.append(surfaces) + frame_output_path = None + + tasks.append((j0, atoms_copy, np.array(frame0.getCoords(), copy=True), + frame_output_path, separate, start_point, return_details, kwargs)) trajectory._nfi = nfi else: if atoms.numCoordsets() > 1: + coordsets = atoms.getCoordsets() for i in range(len(atoms.getCoordsets()[start_frame:stop_frame])): - LOGGER.info("Model: {0}".format(i+start_frame)) - atoms.setACSIndex(i+start_frame) + model_nr = 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) + frame_output_path = str(output_path) + "{0}.pqr".format(model_nr) else: - channels, surfaces = calcChannels(atoms, start_point=start_point, **kwargs) - channels_all.append(channels) - surfaces_all.append(surfaces) + frame_output_path = None + + tasks.append((model_nr, atoms, np.array(coordsets[model_nr], copy=True), + frame_output_path, separate, start_point, return_details, kwargs)) + else: LOGGER.info("Include trajectory or use multi-model PDB file.") + + import multiprocessing + + if len(tasks) == 0: + if return_details: + return channels_all, surfaces_all, details_all + return channels_all, surfaces_all + + if max_proc is None: + max_proc = multiprocessing.cpu_count() + + max_proc = max(1, min(int(max_proc), len(tasks))) + + if max_proc == 1: + results = [_calcChannelsMultipleFramesWorker(task) for task in tasks] + else: + if mp_context is None: + ctx = multiprocessing.get_context() + else: + ctx = multiprocessing.get_context(mp_context) + + with ctx.Pool(processes=max_proc) as pool: + results = pool.map(_calcChannelsMultipleFramesWorker, tasks) + + for result in results: + if return_details: + channels, surfaces, details = result + details_all.append(details) + else: + channels, surfaces = result + + channels_all.append(channels) + surfaces_all.append(surfaces) + + if return_details: + return channels_all, surfaces_all, details_all + return channels_all, surfaces_all @@ -1750,6 +2056,106 @@ def calcSurfaceCavitiesMultipleFrames(atoms, trajectory=None, output_path=None, return cavities_all, surfaces_all +def calcPoresFromChannelsMultipleFrames(channels_all, details_all, output_path=None, + separate=False, max_proc=2, mp_context=None, **kwargs): + """Construct pores for multiple trajectory frames or multi-model PDBs from + channels previously calculated with :func:`calcChannelsMultipleFrames`. + + This function applies :func:`calcPoresFromChannels` independently to each + frame or model. The channel list and calculation details at the same index + must correspond to the same frame/model. + + When using parallel calculations on Windows or macOS, the call to this + function should be placed inside an ``if __name__ == '__main__':`` block + to prevent child processes from executing the main script again. + + :arg channels_all: Lists of channels returned by :func:`calcChannelsMultipleFrames`. + Each element contains channels calculated for one trajectory frame or model. + :type channels_all: list of lists + + :arg details_all: Calculation details returned by :func:`calcChannelsMultipleFrames` + with ``return_details=True``. Each dictionary must contain ``calculator``, + ``simplices``, ``neighbors``, ``vertices``, ``coords``, and ``vdw_radii`` + for the corresponding frame/model. + :type details_all: list of dict + + :arg max_proc: Maximum number of parallel processes used for calculation. + If 1, files are processed serially. If None, all available CPU + cores are used. Default is 2. + :type max_proc: int or None + + :arg mp_context: Multiprocessing start method used for parallel pore + calculations. If `None`, the default method for the operating system + is used. Windows and macOS use the ``'spawn'`` method by default, + whereas Linux typically uses ``'fork'``. Setting + ``mp_context='spawn'`` can be potentially used on Linux, but might + be slower. Available values may include ``'spawn'``, ``'fork'``, + and ``'forkserver'``, depending on the operating system. + Default is `None`. + :type mp_context: str or None + + :arg kwargs: Pore-filtering parameters passed to + :func:`calcPoresFromChannels`, including ``min_end_to_end``, + ``max_end_to_end``, ``min_bottleneck``, ``max_bottleneck``, + ``min_length``, ``max_length``, ``min_volume``, and ``max_volume``. + :type kwargs: dict + + :returns: A list containing pores constructed for each frame/model. + Each element corresponds to the channels and calculation details at + the same index in ``channels_all`` and ``details_all``. + :rtype: list of lists + + Example usage: + channels_all, surfaces_all, details_all = calcChannelsMultipleFrames( + protein, trajectory=dcd, return_details=True) + + pores_all = calcPoresFromChannelsMultipleFrames(channels_all, details_all, + min_end_to_end=40, min_bottleneck=0.7, output_path='poresALL_', separate=True)""" + + if PY3K: + from pathlib import Path + else: + from pathlib2 import Path + + import multiprocessing + + if len(channels_all) != len(details_all): + raise ValueError("channels_all and details_all must contain the same number of frames") + + if output_path is not None: + output_path = Path(output_path) + if output_path.suffix not in ('.pqr', '.pdb'): + output_path = output_path.with_suffix('.pqr') + + tasks = [] + for frame_nr, (channels, details) in enumerate(zip(channels_all, details_all)): + if output_path is not None: + frame_output_path = output_path.with_name( + "{0}_frame{1}{2}".format(output_path.stem, frame_nr, output_path.suffix)) + else: + frame_output_path = None + + tasks.append((frame_nr, channels, details, frame_output_path, separate, kwargs)) + + if max_proc is None: + max_proc = multiprocessing.cpu_count() + + max_proc = max(1, min(int(max_proc), len(tasks))) + + if max_proc == 1: + pores_all = [_calcPoresFromChannelsWorker(task) for task in tasks] + else: + if mp_context is None: + ctx = multiprocessing.get_context() + else: + ctx = multiprocessing.get_context(mp_context) + + with ctx.Pool(processes=max_proc) as pool: + pores_all = pool.map(_calcPoresFromChannelsWorker, tasks) + + return pores_all + + def parseParameters(channels, **kwargs): """Extracts and returns the lengths, bottlenecks, and volumes of each channel in a given list of channels. """ @@ -1829,6 +2235,62 @@ def getChannelParameters(channels, **kwargs): return multi_model_param +def getPoreParameters(pores, **kwargs): + """Extracts and returns the lengths, bottlenecks, and volumes of each + pore in a given list of pores identified using :func:`calcPoresFromChannels`. + + This function iterates through a list of pore objects, extracting the + length, bottleneck, and volume of each pore. These values are collected + into separate lists, which are returned as a tuple for further use. + + :arg pores: A list of pores objects, where each pore has attributes + `length`, `bottleneck`,and `volume`. These attributes represent the + length of the pore, the minimum radius (bottleneck) along its path, + and the total volume of the pore, respectively. + :type pores: 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_pores.txt' suffix. + :type param_file_name: str + + :returns: Three lists containing the lengths, bottlenecks, and volumes of + the pores. + :rtype: tuple (list, list, list) + + Example usage: + lengths, bottlenecks, volumes = getPoreParameters(pores) """ + + multi_model_param = [] + param_file_name = kwargs.get('param_file_name', None) + + try: + results_L_B_V = parseParameters(pores, **kwargs) + lengths, bottlenecks, volumes = results_L_B_V + LOGGER.info("Pore {0}: \t{1} \t{2} \t{3}".format('ID', 'Volume [ų]', + 'Length [Å]', + 'Bottleneck [Å]')) + for i in range(len(lengths)): + LOGGER.info("pore {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 + + except: + for nr_i,i in enumerate(pores): + safe_param_file_name = param_file_name if param_file_name is not None else "" + results = parseParameters(pores[nr_i], param_file_name=safe_param_file_name + str(nr_i)) + multi_model_param.append(results) + + LOGGER.info("Pore {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)) + for i in range(len(lengths)): + LOGGER.info("pore {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 + + def getChannelParametersMultipleFrames(channels_all, **kwargs): """Extract channel parameters for multiple frames or models. @@ -1863,6 +2325,40 @@ def getChannelParametersMultipleFrames(channels_all, **kwargs): return parameters_all +def getPoreParametersMultipleFrames(pores_all, **kwargs): + """Extract pore parameters for multiple frames or models. + + This function is a multi-frame wrapper for :func:`getPoreParameters`. + It extracts pore parameters for each model or trajectory frame separately. + Each element of ``pores_all`` is treated as the list of pores calculated + for one frame/model. + + This function should be used with pores returned by + :func:`calcPoresMultipleFrames`. + + :arg pores_all: list of pore lists returned by + :func:`calcChannelsMultipleFrames`. Each element corresponds to one + model or trajectory frame. + :type pores_all: list + + :arg param_file_name: base name for the output parameter files. If provided, + one file will be written for each model/frame with the frame/model index + added to the file name. + :type param_file_name: str + + :returns: A list of parameter tuples for each model/frame. Each tuple contains + pore lengths, bottlenecks, and volumes. + :rtype: list """ + + parameters_all = [] + for frame_nr, pores in enumerate(pores_all): + LOGGER.info("Frame/model: {0}".format(frame_nr)) + params = getPoreParameters(pores, **kwargs) + parameters_all.append(params) + + return parameters_all + + def parseSurfaceCavityParameters(cavities, **kwargs): """Extract depths, volumes, and tetrahedra counts for surface cavities.""" @@ -2068,19 +2564,23 @@ def convert_lines_to_atomic(atom_lines): return channels_atomic -def getChannelResidueNames(atoms, channels, **kwargs): - '''Provides the resnames and resid of residues that are forming the channel(s). +def getObjectResidueNames(atoms, objects, object_type='channel', **kwargs): + '''Provides the resnames and resid of residues that are forming the object(s). Residues are extracted based on distA which is the distance between FIL atoms - (channel atoms) and protein residues. + (object atoms) and protein residues. Results could be save as txt file by providing the `residues_file_name` parameter. :arg atoms: an Atomic object from which residues are selected - :type atoms: :class:`.Atomic`, :class:`.LigandInteractionsTrajectory` + :type atoms: :class:`.Atomic` - :arg channels: A list of channel objects. Each channel has a method + :arg objects: A list of objects. Each object has a method `getSplines()` that returns the centerline spline and radius spline of - the channel. - :type channels: list + the object. + :type objects: list + + :arg object_type: Type of the object; "channel" or "pore". + Default is "channel". + :type object_type: str :arg distA: Residues will be provided based on this value. default is 4 [Ang] @@ -2088,7 +2588,7 @@ def getChannelResidueNames(atoms, channels, **kwargs): :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 + '_Residues_All_{object_type}.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 @@ -2107,6 +2607,9 @@ def getChannelResidueNames(atoms, channels, **kwargs): raise TypeError('coords must be an object ' 'with `getCoords` method') + if object_type not in ('channel', 'pore'): + raise ValueError("object_type must be 'channel' or 'pore'") + distA = kwargs.pop('distA', 4) residues_file_name = kwargs.pop('residues_file_name', None) @@ -2114,38 +2617,43 @@ def getChannelResidueNames(atoms, channels, **kwargs): if one_letter_aa == True: from prody.atomic.atomic import AAMAP - if isinstance(channels, list): - # Multiple channels + if isinstance(objects, list): + # Multiple objects selected_residues_ch = [] - for i, channel in enumerate(channels): - atoms_protein = getChannelAtoms(channel, atoms) + for i, object in enumerate(objects): + atoms_protein = getChannelAtoms(object, atoms) residues = atoms_protein.select('same residue as exwithin '+str(distA)+' of resname FIL') if residues is not None: resnames = residues.select('name CA').getResnames() if one_letter_aa == True: - resnames_1letter = [AAMAP["HIS"] if aa in ("HSD", "HSP") else AAMAP[aa] for aa in resnames] + resnames_1letter = [AAMAP["HIS"] if aa in ("HSD", "HSP", "HSE", "HID", "HIE", "HIP") + else AAMAP[aa] for aa in resnames] resnames = resnames_1letter resnums = residues.select('name CA').getResnums() residues_info = ["{}{}".format(resname, resnum) for resname, resnum in zip(resnames, resnums)] residues_list = ", ".join(residues_info) - residues_list = 'channel'+str(i)+': '+residues_list + if object_type == "channel": + residues_list = 'channel'+str(i)+': '+residues_list + elif object_type == "pore": + residues_list = "pore"+str(i)+': '+residues_list selected_residues_ch.append(residues_list) else: residues_list = "None" else: - # Single channel analysis in case someone provide channels[0] - atoms_protein = getChannelAtoms(channels, atoms) + # Single object analysis in case someone provide objects[0] + atoms_protein = getChannelAtoms(objects, atoms) residues = atoms_protein.select('same residue as exwithin '+str(distA)+' of resname FIL') selected_residues_ch = [] if residues is not None: resnames = residues.select('name CA').getResnames() if one_letter_aa == True: - resnames_1letter = [AAMAP["HIS"] if aa in ("HSD", "HSP") else AAMAP[aa] for aa in resnames] + resnames_1letter = [AAMAP["HIS"] if aa in ("HSD", "HSP", "HSE", "HID", "HIE", "HIP") + else AAMAP[aa] for aa in resnames] resnames = resnames_1letter resnums = residues.select('name CA').getResnums() @@ -2153,67 +2661,76 @@ def getChannelResidueNames(atoms, channels, **kwargs): residues_list = ", ".join(residues_info) selected_residues_ch.append(residues_list) else: - residues_list = "None" + selected_residues_ch.append("None") if residues_file_name is not None: - output_file = residues_file_name + '_Residues_All_channels.txt' + if object_type == "channel": + output_file = residues_file_name + '_Residues_All_channels.txt' + elif object_type == "pore": + output_file = residues_file_name + '_Residues_All_pores.txt' + with open(output_file, "a") as f_res: for k in selected_residues_ch: f_res.write(("{0}_{1}\n".format(residues_file_name, k))) - LOGGER.info("Channel residues were saved to: {0}".format(output_file)) + if object_type == "channel": + LOGGER.info("Channel residues were saved to: {0}".format(output_file)) + elif object_type == "pore": + LOGGER.info("Pore residues were saved to: {0}".format(output_file)) return selected_residues_ch -def getChannelResidueNamesMultipleFrames(atoms, channels_all, trajectory=None, **kwargs): - """Provides residue names for channels calculated for multiple frames/models. - - This function is a multi-frame wrapper for :func:`getChannelResidueNames`. - For each model/frame, the atomic coordinates are matched with the - corresponding channel prediction. - - :arg atoms: an Atomic object from which residues are selected +def getObjectResidueNamesMultipleFrames(atoms, objects_all, trajectory=None, object_type='channel', **kwargs): + '''Provides the resnames and resid of residues that are forming the object(s) in + multiple frames/models. + Residues are extracted based on distA which is the distance between FIL atoms + (object atoms) and protein residues. + Results could be save as txt file by providing the `residues_file_name` parameter. + + :arg atoms: an Atomic object from which residues are selected :type atoms: :class:`.Atomic` - :arg channels_all: list of channel lists returned by :func:`calcChannelsMultipleFrames`. - :type channels_all: list - + :arg objects_all: A list of objects. Each object has a method + `getSplines()` that returns the centerline spline and radius spline of + the object. + :type objects_all: list + :arg trajectory: optional trajectory object. If provided, coordinates are taken from trajectory frames. If None, a multi-model PDB is assumed and models are selected using ``setACSIndex``. :type trajectory: :class:`.Trajectory` or None - :arg start_frame: first frame/model index. Default is 0. - :type start_frame: int - - :arg stop_frame: last frame/model index. Default is -1, meaning all available - frames/models in ``channels_all``. - :type stop_frame: int - - :arg residues_file_name: base name for output residue files. If provided, - one file will be written for each frame/model. - :type residues_file_name: str - - :arg distA: maximal distance between channel FIL atoms and protein residues. - Default is 4 Å. - :type distA: int, float + :arg object_type: Type of the object; "channel" or "pore". + Default is "channel". + :type object_type: str - :arg one_letter_aa: whether to apply one-letter code to residue names. - Default is False. - :type one_letter_aa: bool - - :returns: A list of residue-name lists for each frame/model. - :rtype: 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_{object_type}.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 + + :arg one_letter_aa: Whether to apply 1-latter code to residue name + by defult is False + :type one_letter_aa: bool ''' start_frame = kwargs.pop('start_frame', 0) stop_frame = kwargs.pop('stop_frame', -1) residues_file_name = kwargs.pop('residues_file_name', None) selected_residues_all = [] + if object_type not in ('channel', 'pore'): + raise ValueError("object_type must be 'channel' or 'pore'") + if trajectory is None: # multi-model PDB - for frame_pos, channels in enumerate(channels_all): + for frame_pos, objects in enumerate(objects_all): model_index = start_frame + frame_pos if stop_frame != -1 and model_index > stop_frame: @@ -2226,9 +2743,13 @@ def getChannelResidueNamesMultipleFrames(atoms, channels_all, trajectory=None, * frame_residues_file_name = residues_file_name + "_model{}".format(model_index) else: frame_residues_file_name = None - - residues = getChannelResidueNames(atoms, channels, - residues_file_name=frame_residues_file_name, **kwargs) + + if object_type == "channel": + residues = getChannelResidueNames(atoms, objects, + residues_file_name=frame_residues_file_name, **kwargs) + elif object_type == "pore": + residues = getPoreResidueNames(atoms, objects, + residues_file_name=frame_residues_file_name, **kwargs) selected_residues_all.append(residues) @@ -2248,7 +2769,7 @@ def getChannelResidueNamesMultipleFrames(atoms, channels_all, trajectory=None, * for frame_pos, frame in enumerate(traj): frame_index = start_frame + frame_pos - if frame_pos >= len(channels_all): + if frame_pos >= len(objects_all): break LOGGER.info("Frame: {0}".format(frame_index)) @@ -2259,9 +2780,13 @@ def getChannelResidueNamesMultipleFrames(atoms, channels_all, trajectory=None, * else: frame_residues_file_name = None - residues = getChannelResidueNames(atoms_copy, channels_all[frame_pos], - residues_file_name=frame_residues_file_name, **kwargs) - + if object_type == "channel": + residues = getChannelResidueNames(atoms_copy, objects_all[frame_pos], + residues_file_name=frame_residues_file_name, **kwargs) + elif object_type == "pore": + residues = getPoreResidueNames(atoms_copy, objects_all[frame_pos], + residues_file_name=frame_residues_file_name, **kwargs) + selected_residues_all.append(residues) if nfi is not None: @@ -2270,6 +2795,136 @@ def getChannelResidueNamesMultipleFrames(atoms, channels_all, trajectory=None, * return selected_residues_all +def getChannelResidueNames(atoms, channels, **kwargs): + '''Provides the resnames and resid of residues that are forming the channel(s). + Residues are extracted based on distA which is the distance between FIL atoms + (channel atoms) and protein residues. + Results could be save as txt file by providing the `residues_file_name` parameter. + + :arg atoms: an Atomic object from which residues are selected + :type atoms: :class:`.Atomic` + + :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 + residues_file_name is PDB+chain for example: '1bbhA'. + :type residues_file_name: str + + :arg one_letter_aa: Whether to apply 1-latter code to residue name + by defult is False + :type one_letter_aa: bool ''' + + return getObjectResidueNames(atoms, channels, object_type='channel', **kwargs) + + +def getPoreResidueNames(atoms, pores, **kwargs): + '''Provides the resnames and resid of residues that are forming the pore(s). + Residues are extracted based on distA which is the distance between FIL atoms + (pore atoms) and protein residues. + Results could be save as txt file by providing the `residues_file_name` parameter. + + :arg atoms: an Atomic object from which residues are selected + :type atoms: :class:`.Atomic` + + :arg pores: A list of pore objects. Each pore has a method + `getSplines()` that returns the centerline spline and radius spline of + the pore. + :type pores: 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_pores.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 + + :arg one_letter_aa: Whether to apply 1-latter code to residue name + by defult is False + :type one_letter_aa: bool ''' + + return getObjectResidueNames(atoms, pores, object_type='pore', **kwargs) + + +def getChannelResidueNamesMultipleFrames(atoms, channels, trajectory=None, **kwargs): + '''Provides the resnames and resid of residues that are forming the channel(s). + Residues are extracted based on distA which is the distance between FIL atoms + (channel atoms) and protein residues. + Results could be save as txt file by providing the `residues_file_name` parameter. + + :arg atoms: an Atomic object from which residues are selected + :type atoms: :class:`.Atomic` + + :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 + residues_file_name is PDB+chain for example: '1bbhA'. + :type residues_file_name: str + + :arg one_letter_aa: Whether to apply 1-latter code to residue name + by defult is False + :type one_letter_aa: bool ''' + + return getObjectResidueNamesMultipleFrames(atoms, channels, trajectory=trajectory, + object_type='channel', **kwargs) + + +def getPoreResidueNamesMultipleFrames(atoms, pores, trajectory=None, **kwargs): + '''Provides the resnames and resid of residues that are forming the pore(s). + Residues are extracted based on distA which is the distance between FIL atoms + (pore atoms) and protein residues. + Results could be save as txt file by providing the `residues_file_name` parameter. + + :arg atoms: an Atomic object from which residues are selected + :type atoms: :class:`.Atomic` + + :arg pores: A list of pore objects. Each pore has a method + `getSplines()` that returns the centerline spline and radius spline of + the pore. + :type pores: 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_pores.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 + + :arg one_letter_aa: Whether to apply 1-latter code to residue name + by defult is False + :type one_letter_aa: bool ''' + + return getObjectResidueNamesMultipleFrames(atoms, pores, trajectory=trajectory, + object_type='pore', **kwargs) + + def getSurfaceCavityResidueNames(atoms, cavities, surface, **kwargs): '''Provides the resnames and resid of residues that form surface cavities. @@ -2662,6 +3317,16 @@ def calcChannelSurfaceOverlaps(**kwargs): 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 mp_context: Multiprocessing start method used for parallel pore + calculations. If `None`, the default method for the operating system + is used. Windows and macOS use the ``'spawn'`` method by default, + whereas Linux typically uses ``'fork'``. Setting + ``mp_context='spawn'`` can be potentially used on Linux, but might + be slower. Available values may include ``'spawn'``, ``'fork'``, + and ``'forkserver'``, depending on the operating system. + Default is `None`. + :type mp_context: str or None :arg output_file_name: The name of the PDB file with overlapping surfaces. :type output_file_name: str @@ -2699,6 +3364,7 @@ def calcChannelSurfaceOverlaps(**kwargs): resolution = kwargs.pop('resolution', 0.5) max_proc = kwargs.pop('max_proc', 2) + mp_context = kwargs.pop('mp_context', None) pqr_files = kwargs.pop('pqr_files', False) if pqr_files == False or pqr_files is None: @@ -2715,7 +3381,18 @@ 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') - + + # PRQ files might be empty + valid_pqr_files = [] + for pqr_file in pqr_files: + if not os.path.isfile(pqr_file) or os.path.getsize(pqr_file) == 0: + LOGGER.warn("Skipping empty PQR file: {0}".format(pqr_file)) + continue + + valid_pqr_files.append(pqr_file) + + pqr_files = valid_pqr_files + if len(pqr_files) == 0: LOGGER.info("No PQR files found.") return None @@ -2742,7 +3419,12 @@ def calcChannelSurfaceOverlaps(**kwargs): 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: + if mp_context is None: + ctx = multiprocessing.get_context() + else: + ctx = multiprocessing.get_context(mp_context) + + with ctx.Pool(processes=max_proc) as pool: for surface in pool.imap_unordered(_surfaceFromPqrWorker, tasks, chunksize=chunksize): merged_surface.update(surface)