CaviTracer: Pores reconstruction

prody/proteins/channels.py

@@ -31,7 +31,8 @@
            'getSurfaceCavityResidueNamesMultipleFrames',
            'getSurfaceCavityParametersMultipleFrames', 
            'getChannelParametersMultipleFrames', '_reportAtomsInputComposition',
-           'getChannelResidueNamesMultipleFrames']
+           'getChannelResidueNamesMultipleFrames', 'calcPoresFromChannels',
+           'showPores']
 
 # Sampling of the enclosure test used to strip the moat (see
 # ChannelCalculator.calcEnclosure). These are constants, not knobs: the enclosure
@@ -366,7 +367,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 +487,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 +826,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 +1104,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 +1411,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 +1434,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,9 +1482,137 @@ 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):
+    """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 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.
+
+    :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
+
+    :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)    """
+
+    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:
+        centerline_spline, radius_spline, length, bottleneck, volume = calculator.processChannel(
+                                                pore_path, vertices, coords, vdw_radii, simplices)
+        pore = Channel(pore_path, centerline_spline, radius_spline, length, bottleneck, volume, 0.0)
+        pores.append(pore)
+
+    return pores
+
+
 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