From b08824a04ed3dec5834c6d87836a61c8a4a8be5d Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Tue, 4 Aug 2026 17:55:54 +0200 Subject: [PATCH 01/10] New clustering functions - Hierarchical, KMedoids --- prody/proteins/interactions.py | 1462 +++++++++++++++++++++++++++++++- 1 file changed, 1461 insertions(+), 1 deletion(-) diff --git a/prody/proteins/interactions.py b/prody/proteins/interactions.py index dfca7a1d3..01aace7c0 100644 --- a/prody/proteins/interactions.py +++ b/prody/proteins/interactions.py @@ -53,7 +53,13 @@ 'calcSminaBindingAffinity', 'calcSminaPerAtomInteractions', 'calcSminaTermValues', 'showSminaTermValues', 'showPairEnergy', 'checkNonstandardResidues', 'saveInteractionsAsDummyAtoms', 'createFoldseekAlignment', 'runFoldseek', 'runDali', - 'runBLAST', 'extractMultiModelPDB', 'calcSignatureInteractions'] + 'runBLAST', 'extractMultiModelPDB', 'calcSignatureInteractions', + 'alignTrajectory', 'calcRMSDfromReference', 'showRMSDfromReference', + 'calcPairwiseRMSD', 'showPairwiseRMSDHeatmap', 'showRMSDHistogram', + 'calcClusterPopulations', 'getCluster', 'getClusterMedoid', + 'calcClusterStatistics', 'calcAllClusterStatistics', 'showClusterStatisticsTable', + 'showClusterRMSDComparison', 'clusterHierarchical', 'showDendrogram', + 'clusterKMedoids', 'writeClusters'] def cleanNumbers(listContacts): @@ -4237,7 +4243,1461 @@ def calcSignatureInteractions(PDB_folder, **kwargs): # Proceed with plotting plot_barh(result, bond_type, n_per_plot=n_per_plot, min_height=min_height) + + + + + +def alignTrajectory(atoms, trajectory, align = 'protein and backbone', select = 'all'): + """ + Aligns each trajectory frame to the reference structure and returns a tuple of the + reference coordinates and the aligned coordinates of the selected atoms. + + The trajectory frames are aligned to the reference structure using the atoms specified by + ``align``. After alignment, the coordinates of the atoms specified by ``select`` are + extracted for every frame. + + + :arg atoms: reference structure used for the alignment. + :type atoms: :class:`prody.Atomic` + + :arg trajectory: trajectory containing the coordinate sets to align + :type trajectory: :class:`prody.Trajectory` + + :arg align: atom selection used to calculate the alignment transformation. + Must be a valid ProDy selection string. + Default is ``"protein and backbone"`` + :type align: str + + :arg select: atom selection whose coordinates are returned. + Must be a valid ProDy selection string. + Default is ``"all"`` + :type select: str + + :returns: a tuple containing: + * ref_coords (numpy.ndarray): coordinates of the selected atoms in the reference structure. + * aligned_coords (numpy.ndarray): aligned coordinates of the selected atoms for every + trajectory frame with shape ``(n_frames, n_atoms, 3)``. + :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) + + Example usage: + >>> pdb = prody.parsePDB("structure.pdb") + >>> dcd = prody.Trajectory("trajectory.dcd") + >>> ref_coords, aligned_coords = prody.alignTrajectory(pdb, dcd, select='resname IOA') + """ + + if trajectory.numAtoms() != atoms.numAtoms(): + raise ValueError("Trajectory atoms count does not match structure atoms count.") + + # Save original coordinates to restore in the end + orig_coords = atoms.getCoords().copy() + + atom_align = atoms.select(align) + if atom_align is None: + raise ValueError(f"No atoms match '{align}' in the structure.") + ref_align = atom_align.copy() + + atom_select = atoms.select(select) + if atom_select is None: + raise ValueError(f"No atoms match '{select}' in the structure.") + ref_coords = atom_select.getCoords().copy() + + trajectory.link(atoms) # linking trajectory to update coordinates frame-by-frame + + n_frames = trajectory.numFrames() + n_atoms = atom_select.numAtoms() + + trajectory.reset() + aligned_coords = np.zeros((n_frames, n_atoms, 3)) + + try: + for i, frame in enumerate(trajectory): + trans = calcTransformation(atom_align, ref_align) + trans.apply(atom_select) + aligned_coords[i] = atom_select.getCoords() + finally: + atoms.setCoords(orig_coords) + trajectory.reset() + + return ref_coords, aligned_coords + + +def calcRMSDfromReference(reference_coords, aligned_coords): + """ + Calculates the RMSD of each aligned trajectory frame from the reference coordinates. + Uses a vectorized approach for better efficiency. + + + :arg reference_coords: reference coordinates with shape ``(n_atoms, 3)``. + :type reference_coords: :class:`numpy.ndarray` + + :arg aligned_coords: aligned coordinates with shape ``(n_frames, n_atoms, 3)``. + Recommended to generate them using :func:`alignTrajectory`. + :type aligned_coords: :class:`numpy.ndarray` + + :returns: an array containing the RMSD of each frame from the reference + :rtype: :class:`numpy.ndarray` + + Example usage: + >>> pdb = prody.parsePDB("structure.pdb") + >>> dcd = prody.Trajectory("trajectory.dcd") + >>> ref_coords, aligned_coords = prody.alignTrajectory(pdb, dcd, select='resname IOA') + >>> rmsd = prody.calcRMSDfromReference(ref_coords, aligned_coords) + """ + + aligned_coords = np.asarray(aligned_coords) + reference_coords = np.asarray(reference_coords) + + if reference_coords.ndim != 2: + raise ValueError(f"reference_coords must be a 2D array of shape (n_atoms, 3), but got {reference_coords.shape}") + + if aligned_coords.ndim != 3: + raise ValueError(f"aligned_coords must have shape (n_frames, n_atoms, 3), but got {aligned_coords.shape}.") + + if reference_coords.shape != aligned_coords.shape[1:]: + raise ValueError(f"Incompatible shapes: reference is {reference_coords.shape}, but aligned frames have {aligned_coords.shape[1:]}.") + + # Vectorized RMSD calculation + sd = (aligned_coords - reference_coords)**2 + msd = np.mean(sd.sum(axis=2), axis=1) + rmsd = np.sqrt(msd) + + return rmsd + + +def showRMSDfromReference(rmsd_array, title = "Frame-to-Reference RMSD", + color = "#36454F", lw = 0.75, label = None, ax = None): + """ + Plots the RMSD of an aligned trajectory from the reference on the current axis. + + + :arg rmsd_array: one-dimensional array containing the RMSD values (in Å) of each + trajectory frame relative to the reference structure. + This array is typically generated using :func:`calcRMSDfromReference`. + :type rmsd_array: :class:`numpy.ndarray` + + :arg title: the title of the generated plot. + Default is ``"Frame-to-Reference RMSD"`` + :type title: str + + :arg color: the color of the plot. + Default is ``"#36454F"`` + :type color: str + + :arg lw: the width of the line. + Default is ``0.75`` + :type lw: float + + :arg label: the label for the plot. + Default is ``None`` + :type label: str + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :returns: the Matplotlib axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes`, None + + Example usage: + >>> import matplotlib.pyplot as plt + >>> rmsd_run1 = prody.calcRMSDfromReference(ref_coords, aligned_coords1) + >>> rmsd_run2 = prody.calcRMSDfromReference(ref_coords, aligned_coords2) + >>> plt.figure(figsize=(8, 6)) + >>> prody.showRMSDfromReference(rmsd_run1, label="Run 1") + >>> prody.showRMSDfromReference(rmsd_run2, color="red", label="Run 2") + >>> plt.show() + """ + + import matplotlib.pyplot as plt + + + if ax is None: + ax = plt.gca() + + ax.plot(rmsd_array, color = color, lw = lw, label = label) + ax.set_xlabel("# Frame") + ax.set_ylabel("RMSD from Reference [Å]") + ax.set_title(title) + ax.grid(alpha=0.3) + + if label is not None: + ax.legend() + + return ax + + +def calcPairwiseRMSD(aligned_coords): + """ + Calculates the frame-to-frame pairwise RMSD matrix using aligned structures + and returns a symmetric distance matrix. + Uses a vectorized approach for better efficiency. + + + :arg aligned_coords: aligned coordinates with shape ``(n_frames, n_atoms, 3)``. + Recommended to generate them using :func:`alignTrajectory`. + :type aligned_coords: :class:`numpy.ndarray` + + :returns: symmetric matrix containing the pairwise RMSD between all frames with shape + ``(n_frames, n_frames)``. + :rtype: :class:`numpy.ndarray` + + Example usage: + >>> ref_coords, aligned_coords = prody.alignTrajectory(pdb, dcd, select='resname IOA') + >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) + """ + + from scipy.spatial.distance import cdist + + + aligned_coords = np.asarray(aligned_coords) + + if aligned_coords.ndim != 3 or aligned_coords.shape[2] != 3: + raise ValueError(f"aligned_coords must have shape (n_frames, n_atoms, 3), but got {aligned_coords.shape}.") + + if aligned_coords.shape[0] == 0: + raise ValueError("aligned_coords contains no frames.") + + if aligned_coords.shape[1] == 0: + raise ValueError("aligned_coords contains no atoms.") + + + n_frames, n_atoms, _ = aligned_coords.shape + + # Flatten frame coordinates (n_frames, n_atoms, 3) -> (n_frames, n_atoms * 3) + coords_flat = aligned_coords.reshape(n_frames, 3 * n_atoms) + + euclidean_dists = cdist(coords_flat, coords_flat, metric = 'euclidean') + distance_matrix = euclidean_dists / np.sqrt(n_atoms) + + return distance_matrix + + +def showPairwiseRMSDHeatmap(distance_matrix, title = "Pairwise RMSD Distance Matrix", + label = "RMSD [Å]", cmap = "viridis", ax = None): + """ + Plots the heatmap of the frame-to-frame pairwise RMSDs using the distance matrix + on the current axis. + + + :arg distance_matrix: two-dimensional array constituting the pairwise distance matrix. + Typically generated using :func:`calcPairwiseRMSD`. + :type distance_matrix: :class:`numpy.ndarray` + + :arg title: the title of the generated plot. + Default is ``"Pairwise RMSD Distance Matrix"`` + :type title: str + + :arg label: the label for the heatmap legend. + Default is ``"RMSD [Å]"`` + :type label: str + + :arg cmap: the colormap for the heatmap. + Default is ``"viridis"`` + :type cmap: str + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :returns: the Matplotlib axes containing the heatmap. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) + >>> plt.figure(figsize=(8, 6)) + >>> prody.showPairwiseRMSDHeatmap(distance_matrix) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + + + if distance_matrix.ndim != 2: + raise ValueError(f"distance_matrix must be a 2D array, but got shape {distance_matrix.shape}.") + if distance_matrix.shape[0] != distance_matrix.shape[1]: + raise ValueError("distance_matrix must be square.") + + if ax is None: + ax = plt.gca() + + im = ax.imshow(distance_matrix, cmap = cmap) + ax.figure.colorbar(im, ax = ax, label = label) + + ax.set_xlabel("# Frame") + ax.set_ylabel("# Frame") + ax.set_title(title) + + return ax + + +def showRMSDHistogram(rmsd_data, bins = 50, title = 'Distribution of RMSDs', + xlabel = 'RMSD [Å]', ylabel = 'Frequency', label = None, + kde = False, element = 'bars', stat = 'count', + color = 'teal', edgecolor = 'black', lw = 0.8, alpha = 0.5, ax = None): + """ + Plots the distribution of RMSD values on the current axis. + The input may be either a one-dimensional RMSD array (e.g., frame-to-reference or intra-cluster RMSDs) + or a 2D pairwise RMSD distance matrix. + + + :arg rmsd_data: 1D array of RMSDs or 2D pairwise distance matrix. + If 2D, the upper triangle (excluding diagonal) is automatically extracted. + Recommended to generate them using :func:`calcRMSDfromReference` or + :func:`calcPairwiseRMSD` respectively. + :type rmsd_data: :class:`numpy.ndarray` + + :arg bins: number of histogram bins. + Default is ``50`` + :type bins: int + + :arg title: the title of the generated plot. + Default is ``'Distribution of RMSDs'`` + :type title: str + + :arg xlabel: the label for the X axis. + Default is ``'RMSD [Å]'`` + :type xlabel: str + + :arg ylabel: the label for the Y axis. + Default is ``'Frequency'`` + :type ylabel: str + + :arg label: the label for the plot legend. + Default is ``None`` + :type label: str + + :arg kde: whether to plot a kernel density estimate. + Default is ``False`` + :type kde: bool + + :arg element: visual representation of the histogram bins ('bars', 'step', or 'poly'). + Default is ``'bars'`` + :type element: str + + :arg stat: aggregate statistic to compute in each bin ('count', 'frequency', 'probability', 'percent', 'density'). + Default is ``'count'`` + :type stat: str + + :arg color: the color of the plot. + Default is ``'teal'`` + :type color: str + + :arg edgecolor: the color of the bin edges. + Default is ``'black'`` + :type edgecolor: str + + :arg lw: the line width of the bin edges. + Default is ``0.8`` + :type lw: float + + :arg alpha: the transparency of the bins. + Default is ``0.5`` + :type alpha: float + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :returns: the Matplotlib axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) + >>> plt.figure(figsize = (8, 6)) + >>> prody.showRMSDHistogram(distance_matrix, bins=70, kde=True) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + try: + import seaborn as sns + except ImportError: + raise ImportError("The 'seaborn' package is required to display the Histogram." + "\nPlease install it using 'pip install seaborn'." + "\nAlternatively, use standard matplotlib.pyplot.hist() for basic plots.") + + + data_array = np.asarray(rmsd_data) + + if data_array.ndim == 2: + if data_array.shape[0] != data_array.shape[1]: + raise ValueError(f"Pairwise RMSD matrix must be square, but got shape {data_array.shape}") + rmsd_values = data_array[np.triu_indices_from(data_array, k=1)] + elif data_array.ndim == 1: + rmsd_values = data_array + else: + raise ValueError(f"Expected 1D or 2D array, but got shape {data_array.shape}") + + if ax is None: + ax = plt.gca() + + sns.histplot(rmsd_values, bins = bins, element = element, stat = stat, kde = kde, + alpha = alpha, color = color, edgecolor = edgecolor, linewidth = lw, + label = label, ax = ax) + + ax.set_title(title) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + + if label is not None: + ax.legend() + + ax.grid(axis = 'y', alpha = 0.3) + + return ax + + +def calcClusterPopulations(cluster_ids): + """ + Uses the cluster IDs array, assigning each frame to a cluster, to calculate the population + of each cluster and its corresponding percentage. + + + :arg cluster_ids: a one-dimensional array matching each frame to a cluster. + :type cluster_ids: :class:`numpy.ndarray` + + :returns: a dictionary mapping each cluster ID to its population statistics. + Each value is a dictionary with the keys ``"count"`` and ``"pct"``. + :rtype: dict + + Example usage: + >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) + >>> populations = prody.calcClusterPopulations(cluster_ids) + >>> print(f"In cluster 1: {populations[1]['pct']:.2f}% of the objects") + """ + + cluster_ids = np.asarray(cluster_ids) + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}") + if cluster_ids.size == 0: + raise ValueError("cluster_ids is empty.") + + total_frames = len(cluster_ids) + clusters, frequencies = np.unique(cluster_ids, return_counts=True) + + population_data = {} + for c, n in zip(clusters, frequencies): + population_data[int(c)] = {'count': int(n), 'pct': (n / total_frames) * 100} + + return population_data + + +def getCluster(cluster_ids, cluster_number): + """ + Returns the indices of the frames grouped in a specific cluster. + + + :arg cluster_ids: a one-dimensional array assigning each frame to a cluster. + :type cluster_ids: :class:`numpy.ndarray` + + :arg cluster_number: the ID of the cluster whose members we want to return. + :type cluster_number: int + + :returns: an array of the trajectory indices assigned to the specified cluster. + :rtype: :class:`numpy.ndarray` + + Example usage: + >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) + >>> cluster1_indices = prody.getCluster(cluster_ids, 1) + """ + + cluster_ids = np.asarray(cluster_ids) + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") + if cluster_ids.size == 0: + raise ValueError("cluster_ids is empty.") + + if isinstance(cluster_number, bool) or not isinstance(cluster_number, (int, np.integer)): + raise TypeError(f"cluster_number must be an integer, but got {type(cluster_number).__name__}.") + + cluster_indices = np.where(cluster_ids == cluster_number)[0] + + if cluster_indices.size == 0: + raise ValueError(f'No frames belong to Cluster {cluster_number}') + + return cluster_indices + + +def getClusterMedoid(distance_matrix, cluster_indices = None, cluster_ids = None, cluster_number = None): + """ + Determines the medoid (representative) of the cluster from the minimum total pairwise + distance to all other cluster members. + Returns a dictionary containing the medoid index within the cluster ("local") and in + the full trajectory ("global"). + + The cluster can be specified either by: + 1. cluster_indices (recommended and can be generated by :func:`getCluster`) or + 2. cluster_ids and cluster_number + + + :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg cluster_indices: a one-dimensional array of the indices of the cluster members + Default is ``None`` + :type cluster_indices: :class:`numpy.ndarray` + + :arg cluster_ids: a one-dimensional array of the cluster IDs per element + Default is ``None`` + :type cluster_ids: :class:`numpy.ndarray` + + :arg cluster_number: the ID (number) of the cluster we want to investigate + Default is ``None`` + :type cluster_number: int + + :returns: a dictionary of the medoid index, with keys "local" for the index within the + cluster and "global" for the index within the full trajectory + :rtype: dict + + Example usage: + >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) + >>> cluster1_indices = prody.getCluster(cluster_ids, 1) + >>> cluster1_medoid = prody.getClusterMedoid(distance_matrix, cluster1_indices) + >>> print(f"Medoid of cluster 1 is # {cluster1_medoid['local']} in the cluster.") + """ + + cluster_indices = _resolveClusterIndices(cluster_indices=cluster_indices, + cluster_ids=cluster_ids, + cluster_number=cluster_number) + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + n_frames = distance_matrix.shape[0] + if np.any((cluster_indices < 0) | (cluster_indices >= n_frames)): + raise ValueError("cluster_indices contain indices outside the bounds of the distance matrix.") + + cluster_distance_matrix = distance_matrix[np.ix_(cluster_indices, cluster_indices)] + + sum_dist = cluster_distance_matrix.sum(axis=1) + medoid_local = int(np.argmin(sum_dist)) + medoid_global = int(cluster_indices[medoid_local]) + + return {"local": medoid_local, "global": medoid_global} + + +def calcClusterStatistics(distance_matrix, cluster_indices = None, cluster_medoid = None, + cluster_ids = None, cluster_number = None): + """ + Calculates descriptive statistics for a single cluster. + Returns a dictionary with the statistics quantity as a key. + + The cluster can be specified either by: + 1. cluster_indices (recommended and can be generated by :func:`getCluster`) or + 2. cluster_ids and cluster_number + + + :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg cluster_indices: a one-dimensional array of the indices of the cluster members + Default is ``None`` + :type cluster_indices: :class:`numpy.ndarray` + + :arg cluster_medoid: the dictionary of the cluster medoid. + Expects the 'cluster_medoid' dictionary generated by the function + :func:`getClusterMedoid`. + Default is ``None`` + :type cluster_medoid: dict + + :arg cluster_ids: a one-dimensional array of the cluster IDs per element + Default is ``None`` + :type cluster_ids: :class:`numpy.ndarray` + + :arg cluster_number: the ID (number) of the cluster we want to investigate + Default is ``None`` + :type cluster_number: int + + :returns: a dictionary of statistics quantities with their names as keys + :rtype: dict + + Example usage: + >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) + >>> cluster1_indices = prody.getCluster(cluster_ids, 1) + >>> cluster_stats1 = prody.calcClusterStatistics(distance_matrix, cluster1_indices) + >>> print(f"Mean distance from cluster Medoid: {cluster_stats1['mean']:.2f} ± {cluster_stats1['std']:.2f} [Å]") + """ + + cluster_indices = _resolveClusterIndices(cluster_indices, cluster_ids, cluster_number) + population = len(cluster_indices) + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + if cluster_medoid is None: + cluster_medoid = getClusterMedoid(distance_matrix, cluster_indices=cluster_indices) + else: + if not isinstance(cluster_medoid, dict): + raise TypeError("cluster_medoid must be a dictionary with keys 'global' and 'local'.") + + required_keys = {"global", "local"} + if not required_keys.issubset(cluster_medoid.keys()): + missing = required_keys - cluster_medoid.keys() + raise ValueError(f"Expected keys are 'global' and 'local'.\ncluster_medoid is missing required key(s): {missing}.") + + if cluster_medoid["global"] not in cluster_indices: + raise ValueError("Index not found.\ncluster_medoid does not belong to the specified cluster.") + if not (0 <= cluster_medoid["local"] < population): + raise ValueError("Local index is out of bounds.") + + if cluster_indices[cluster_medoid["local"]] != cluster_medoid["global"]: + raise ValueError("cluster_medoid['local'] and cluster_medoid['global'] are inconsistent.") + + distance_to_medoid = distance_matrix[cluster_indices, cluster_medoid["global"]] + total_frames = distance_matrix.shape[0] + + stats = { + "population" : population, + "pct" : (population / total_frames) * 100, + "medoid_global" : cluster_medoid["global"], + "medoid_local" : cluster_medoid["local"], + "distances" : distance_to_medoid, + "mean" : np.mean(distance_to_medoid), + "std" : np.std(distance_to_medoid), + "median" : np.median(distance_to_medoid), + "iqr" : np.percentile(distance_to_medoid, 75) - np.percentile(distance_to_medoid, 25), + "p95" : np.percentile(distance_to_medoid, 95), + "max" : np.max(distance_to_medoid) + } + + if cluster_number is not None: + stats["cluster"] = cluster_number + + return stats + + +def calcAllClusterStatistics(distance_matrix, cluster_ids): + """ + Calculates descriptive statistics for all clusters. + Returns a list of dictionaries, each corresponding to a cluster with the statistic quantities + as keys. + + + :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg cluster_ids: a one-dimensional array of the cluster IDs per element + :type cluster_ids: :class:`numpy.ndarray` + + :returns: a list of dictionaries, one dictionary for each cluster with the statistic + quantities as keys. + :rtype: list of dict + + Example usage: + >>> cluster_ids, _= prody.clusterHierarchical(distance_matrix, 4) + >>> all_stats = prody.calcAllClusterStatistics(distance_matrix, cluster_ids) + """ + + cluster_ids = np.asarray(cluster_ids) + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + if len(cluster_ids) != distance_matrix.shape[0]: + raise ValueError(f"Dimension mismatch: cluster_ids length ({len(cluster_ids)}) " + f"does not match distance_matrix frames ({distance_matrix.shape[0]}).") + + clusters = np.unique(cluster_ids) + all_stats = [] + + for c in clusters: + cluster_stats = calcClusterStatistics(distance_matrix, cluster_ids = cluster_ids, cluster_number = int(c)) + all_stats.append(cluster_stats) + + return all_stats + + +def _resolveClusterIndices(cluster_indices, cluster_ids, cluster_number): + """ + Checks whether the input is cluster_indices or cluster_ids and cluster_number, and for the latter returns + the cluster indices. + """ + + has_indices = cluster_indices is not None + has_ids_info = (cluster_ids is not None) and (cluster_number is not None) + + if has_indices == has_ids_info: + raise ValueError("Provide exactly one of 'cluster_indices' or both 'cluster_ids' and 'cluster_number'.") + + if has_indices: + cluster_indices = np.asarray(cluster_indices) + + if cluster_indices.ndim != 1: + raise ValueError(f"cluster_indices must be a 1D array, but shape {cluster_indices.shape} was given.") + if cluster_indices.size == 0: + raise ValueError("cluster_indices must not be empty.") + else: + cluster_ids = np.asarray(cluster_ids) + cluster_indices = getCluster(cluster_ids, cluster_number) + + return cluster_indices + + +def _validateDistanceMatrix(distance_matrix): + """ + Validates that a distance matrix is a non-empty square 2D NumPy array. + """ + + distance_matrix = np.asarray(distance_matrix) + if distance_matrix.ndim != 2: + raise ValueError(f"The distance matrix must be a 2D array, but got shape {distance_matrix.shape}.") + if distance_matrix.shape[0] != distance_matrix.shape[1]: + raise ValueError(f"distance_matrix must be square, but got shape {distance_matrix.shape}.") + if distance_matrix.size == 0: + raise ValueError("distance_matrix cannot be empty.") + + return distance_matrix + + +def showClusterStatisticsTable(all_stats, dissimilarity = "RMSD", units = "Å", show = True): + """ + Prints the cluster statistics in a table, where each column corresponds to a cluster. + Expects either the 'stats' dictionary generated by :func:`calcClusterStatistics` or + the 'all_stats' list of dictionaries generated by :func:`calcAllClusterStatistics`. + + + :arg all_stats: a dictionary, or list of dictionaries, with the cluster descriptive statistics + Recommended to generate the dictionaries from the functions :func:`calcClusterStatistics` + of a single cluster or :func:`calcAllClusterStatistics`. + :type all_stats: dict, or list of dict + + :arg dissimilarity: the dissimilarity measure + Default is `"RMSD"` + :type dissimilarity: str + + :arg units: the units of the dissimilarity measure + Default is `"Å"` + :type units: str + + :arg show: whether to print the table + Default is `True` + :type show: bool + + :returns: a dictionary containing the table, the row labels and the headers, with their names as keys + :rtype: dict + + Example usage: + >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) + >>> all_stats = prody.calcAllClusterStatistics(distance_matrix, cluster_ids) + >>> prody.showClusterStatisticsTable(all_stats); + """ + + try: + from tabulate import tabulate + except ImportError: + raise ImportError("The 'tabulate' package is required to display the table. " + "Please install it using 'pip install tabulate'.") + + + if isinstance(all_stats, dict): + all_stats = [all_stats] + elif not isinstance(all_stats, list): + raise TypeError(f"all_stats must be a dict or list of dicts, but got {type(all_stats).__name__}.") + + if not all_stats: + raise ValueError("all_stats cannot be empty.") + + required = {"population", "pct", "medoid_global", "medoid_local", "mean", "std", + "median", "iqr", "p95", "max"} + + headers = [] + for i, stat in enumerate(all_stats): + if not isinstance(stat, dict): + raise TypeError(f"Expected a dictionary in all_stats, but got {type(stat).__name__} at index {i}.") + + cluster_id = stat.get('cluster', i + 1) + headers.append(f"Cluster {int(cluster_id)}") + + missing = required - stat.keys() + if missing: + raise ValueError(f"Statistics dictionary is missing required keys: {sorted(missing)}.") + + metrics = [ + ("population", "Total Frames"), + ("pct", "Population Percentage (%)"), + ("medoid_global", "Medoid Frame (Global)"), + ("medoid_local", "Medoid Frame (Within Cluster)"), + ("mean", f"Mean {dissimilarity} [{units}]"), + ("std", f"Std {dissimilarity} [{units}]"), + ("median", f"Median {dissimilarity} [{units}]"), + ("iqr", f"IQR [{units}]"), + ("p95", f"95th Percentile [{units}]"), + ("max", f"Max {dissimilarity} [{units}]") + ] + + table = [] + row_labels = [] + + for key, label in metrics: + row_labels.append(label) + table.append([stat[key] for stat in all_stats]) + + if not isinstance(show, bool): + raise TypeError("show must be a bool.") + + if show: + print(tabulate(table, headers=headers, showindex=row_labels, + tablefmt="fancy_grid", floatfmt=".4f", stralign='center')) + + return {"table": table, "row_labels": row_labels, "headers": headers} + + +def showClusterRMSDComparison(all_stats, bins = 50, title = 'RMSD Distributions', + xlabel = 'RMSD to Medoid [Å]', ylabel = 'Frequency', + kde = True, element = 'step', stat = 'count', + alpha = 0.5, lw = 1.5, ax = None): + """ + Overlays the internal RMSD distributions of all clusters onto a single plot + on the current axis. + + Expects the list of cluster statistics dictionaries generated by :func:`calcAllClusterStatistics`. + + + :arg all_stats: list of cluster statistics dictionaries. Each dictionary + must contain the keys ``"cluster"`` and ``"distances"``. + :type all_stats: list of dict + + :arg bins: number of histogram bins. + Default is ``50`` + :type bins: int + + :arg title: the title of the generated plot. + Default is ``'Intra-Cluster RMSD Distributions'`` + :type title: str + + :arg xlabel: the label for the X axis. + Default is ``'RMSD to Medoid [Å]'`` + :type xlabel: str + + :arg ylabel: the label for the Y axis. + Default is ``'Frequency'`` + :type ylabel: str + + :arg kde: whether to plot kernel density estimates. + Default is ``True`` + :type kde: bool + + :arg element: visual representation of the histogram bins ('step', 'bars', or 'poly'). + Default is ``'step'`` + :type element: str + + :arg stat: aggregate statistic to compute in each bin. + Default is ``'count'`` + :type stat: str + + :arg alpha: the transparency of the bins. + Default is ``0.5`` + :type alpha: float + + :arg lw: line width of the distribution lines. + Default is ``1.5`` + :type lw: float + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :returns: the Matplotlib axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> all_stats = prody.calcAllClusterStatistics(distance_matrix, clusterIDs_array) + >>> plt.figure(figsize = (8, 6)) + >>> prody.showClusterRMSDComparison(all_stats) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + try: + import seaborn as sns + except ImportError: + raise ImportError("The 'seaborn' package is required to display the Histogram." + "\nPlease install it using 'pip install seaborn'." + "\nAlternatively, use standard matplotlib.pyplot.hist() for basic plots.") + + + if len(all_stats) == 0: + raise ValueError("all_stats cannot be empty.") + + if ax is None: + ax = plt.gca() + + for cluster_stats in all_stats: + sns.histplot(cluster_stats["distances"], bins = bins, element = element, stat = stat, kde = kde, + alpha = alpha, linewidth = lw, label = f"Cluster {cluster_stats['cluster']}", + ax = ax) + + ax.set_title(title) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + + ax.legend(title = "Clusters") + ax.grid(axis = 'y', alpha = 0.3) + + return ax + + +def clusterHierarchical(distance_matrix, method='average', cutoff=None): + """ + Performs bottom-up hierarchical clustering from a pairwise distance matrix. + + Note that the resulting cluster IDs are 1-based (starting at 1, not 0). + + If cutoff is ``None``, it is automatically chosen as the midpoint of the largest + gap between consecutive linkage distances. + + Recommendation: Plot the dendrogram using showDendrogram() and choose the cutoff + manually whenever possible. + + + :arg distance_matrix: either a one-dimensional condensed distance matrix or + a two-dimensional pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg method: linkage criterion used when constructing the linkage matrix. + Default is ``'average'`` + :type method: str + + :arg cutoff: the cutoff distance that determines the number of clusters. + Default is ``None`` and the midpoint of the largest linkage gap + is used automatically. + :type cutoff: float + + :returns: + * a one-dimensional array containing the cluster ID for each object. + IDs are 1-indexed (1 to number of clusters for the cutoff) + * the hierarchical clustering linkage matrix + :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) + + Example usage: + >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) + >>> cluster_ids, linkage_matrix = prody.clusterHierarchical(distance_matrix) + """ + + from scipy.cluster.hierarchy import linkage, fcluster + + + condensed_distance_matrix = _condenseDistanceMatrix(distance_matrix) + linkage_matrix = linkage(condensed_distance_matrix, method=method) + + if cutoff is None: + cutoff = _calcAutoCutoff(linkage_matrix) + elif isinstance(cutoff, str): + raise ValueError(f"Invalid input for cutoff: '{cutoff}'. Use numeric value or None.") + + cluster_ids = fcluster(linkage_matrix, t=cutoff, criterion='distance') + + return cluster_ids, linkage_matrix + + +def showDendrogram(distance_matrix = None, linkage_matrix = None, method = 'average', cutoff = None, + truncate_mode = None, p = 30, title = "Hierarchical Clustering Dendrogram", + ylabel = "RMSD [Å]", ax = None): + """ + Plots a hierarchical clustering dendrogram using either a pairwise RMSD distance matrix + or a pre-computed linkage matrix. + + By default, the complete dendrogram is shown. To display a truncated dendrogram, + use ``truncate_mode`` together with ``p``. + + Exactly one of ``distance_matrix`` or ``linkage_matrix`` must be provided. + + + :arg distance_matrix: one-dimensional condensed or two-dimensional square pairwise + distance matrix. Typically generated using :func:`calcPairwiseRMSD`. + Default is ``None`` + :type distance_matrix: :class:`numpy.ndarray` + + :arg linkage_matrix: pre-computed hierarchical clustering linkage matrix. + Default is ``None`` + :type linkage_matrix: :class:`numpy.ndarray` + + :arg method: linkage criterion used when constructing the linkage matrix from + ``distance_matrix``. + Default is ``"average"`` + :type method: str + + :arg cutoff: dendrogram color threshold. If ``"auto"``, the cutoff is placed at the midpoint + of the largest linkage gap. + Default is ``None`` + :type cutoff: float, str + + :arg truncate_mode: dendrogram truncation mode passed to :func:`scipy.cluster.hierarchy.dendrogram`. + Default is ``None`` + :type truncate_mode: str, None + + :arg p: truncation parameter used together with ``truncate_mode``. + Default is ``30`` + :type p: int + + :arg title: title of the generated plot. + Default is ``"Hierarchical Clustering Dendrogram"`` + :type title: str + + :arg ylabel: Distance matric label for the y-axis + Default is ``"RMSD [Å]"`` + :type ylabel: str + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :returns: the matplotlib axes and the linkage matrix. + :rtype: tuple(:class:`matplotlib.axes.Axes`, :class:`numpy.ndarray`) + + Example usage: + >>> import matplotlib.pyplot as plt + >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) + >>> plt.figure(figsize=(12, 8)) + >>> ax, linkage = prody.showDendrogram(distance_matrix, cutoff='auto', truncate_mode='lastp', p=30) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + from scipy.cluster.hierarchy import dendrogram, linkage + + + if (distance_matrix is None) == (linkage_matrix is None): + raise ValueError("Provide exactly one of 'distance_matrix' or 'linkage_matrix'.") + + if linkage_matrix is None: + condensed_distance_matrix = _condenseDistanceMatrix(distance_matrix) + linkage_matrix = linkage(condensed_distance_matrix, method=method) + else: + linkage_matrix = np.asarray(linkage_matrix) + + if linkage_matrix.ndim != 2 or linkage_matrix.shape[1] != 4 or linkage_matrix.shape[0] == 0: + raise ValueError("linkage_matrix must be a non-empty array with shape (n_samples-1, 4).") + + if cutoff == "auto": + cutoff = _calcAutoCutoff(linkage_matrix) + elif isinstance(cutoff, str): + raise ValueError(f"Invalid string for cutoff: '{cutoff}'. Use 'auto' or a numeric value.") + + if ax is None: + ax = plt.gca() + + with plt.rc_context({"lines.linewidth": 0.6}): + dendrogram(linkage_matrix, truncate_mode = truncate_mode, p = p, color_threshold = cutoff, + no_labels = True, ax = ax) + + ax.set_xlabel("# Frame") + ax.set_ylabel(ylabel) + ax.set_title(title) + + if cutoff is not None: + ax.axhline(cutoff, color = "black", linestyle = "--", linewidth = 1.2, label = f"Cutoff ({cutoff:.2f} Å)") + ax.legend() + + return ax, linkage_matrix + + +def _calcAutoCutoff(linkage_matrix): + """ + Calculates the optimal distance cutoff based on the midpoint of the largest gap + in the linkage matrix. + """ + + linkage_matrix = np.asarray(linkage_matrix) + merge_dist = linkage_matrix[:, 2] # locations of the merges + + if len(merge_dist) == 0: + raise ValueError("At least two frames are required for clustering") + elif len(merge_dist) == 1: + return merge_dist[0] + 1.0 + + gaps = np.diff(merge_dist) + largest_gap_idx = np.argmax(gaps) + + return (merge_dist[largest_gap_idx] + merge_dist[largest_gap_idx + 1]) / 2.0 + + +def _condenseDistanceMatrix(distance_matrix): + """ + Converts a pairwise distance matrix into condensed form. + + Accepts either: + * a condensed one-dimensional distance matrix, or + * a square two-dimensional distance matrix. + + Returns the condensed distance matrix suitable for scipy.cluster.hierarchy.linkage. + """ + from scipy.spatial.distance import squareform + + + distance_matrix = np.asarray(distance_matrix) + + if distance_matrix.size == 0: + raise ValueError("distance_matrix cannot be empty.") + + if distance_matrix.ndim == 2: + if distance_matrix.shape[0] != distance_matrix.shape[1]: + raise ValueError(f"distance_matrix must be square, but got shape {distance_matrix.shape}.") + condensed_distance_matrix = squareform(distance_matrix) + elif distance_matrix.ndim == 1: + condensed_distance_matrix = distance_matrix + else: + raise ValueError("distance_matrix must be either a condensed 1D array or a square 2D array.") + + return condensed_distance_matrix + + +def clusterKMedoids(distance_matrix, k, method = 'alternate', method_sklearn = 'pam', + initial_medoids = None, seed = None, max_iter = 100, n_init = 10): + """ + Performs K-Medoids clustering using various algorithms. + + This function acts as a facade, routing the clustering task to the specified backend + ('alternate', 'pam', or 'sklearn'). 'alternate' is generally faster, while 'pam' typically + finds solutions with lower cost. 'sklearn' needs the sklearn_extra.cluster module to be installed. + + Note that the resulting cluster IDs are 1-based (starting at 1, not 0). + + + :arg distance_matrix: square, symmetric matrix of pairwise distances between all objects + :type distance_matrix: :class:`numpy.ndarray` + + :arg k: prespecified number of clusters to form + :type k: int + + :arg method: the clustering algorithm to use. + Options are 'alternate' (custom Alternating Medoids), 'pam' (custom PAM), 'sklearn' (sklearn_extra) + Default is 'alternate' + :type method: str + + :arg method_sklearn: the specific method to pass to sklearn_extra if method = 'sklearn' is chosen. + Options are 'pam', 'alternate' + Default is 'pam' + :type method_sklearn: str + + :arg initial_medoids: one-dimensional array of indices to use as starting medoids. + Default is ``None`` and the starting medoids are picked randomly. + It is not supported for the 'sklearn' method. + :type initial_medoids: :class:`numpy.ndarray` + + :arg seed: random seed for reproducibility + Default is ``None`` + :type seed: int + + :arg max_iter: maximum number of iterations per a single run + Default is ``100`` + :type max_iter: int + + :arg n_init: number of times the algorithm will be run with different initial medoids + and the best result with the lowest cost is returned + + :returns: a tuple of: + * an one-dimensional array containing the cluster ID for each object. + IDs are 1-indexed (1 to k) + * a one-dimensional array of shape (k,) containing the indices of the final cluster + medoids + * the final sum of distances from each point to its nearest medoid + :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`, float) + + Example usage: + >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) + >>> cluster_ids, medoids, _ = prody.clusterKMedoids(distance_matrix, 4, method = 'alternate', + seed = 42, n_init = 30) + """ + + distance_matrix = _validateDistanceMatrix(distance_matrix) + n = distance_matrix.shape[0] + + if not isinstance(k, (int, np.integer)): + raise TypeError(f"k must be an integer, got {type(k).__name__}") + if k <= 0 or k > n: + raise ValueError(f"k must be between 1 and {n} (the number of points).") + + if not isinstance(max_iter, (int, np.integer)) or max_iter <= 0: + raise ValueError("max_iter must be a positive integer.") + + if not isinstance(n_init, (int, np.integer)) or n_init <= 0: + raise ValueError("n_init must be a positive integer.") + + if seed is not None and not isinstance(seed, (int, np.integer)): + raise TypeError("seed must be an integer or None.") + + if initial_medoids is not None: + if method == 'sklearn': + raise ValueError("The 'sklearn' method does not support custom 'initial_medoids'. Use 'alternate' or 'pam' instead.") + + initial_medoids = np.asarray(initial_medoids, dtype=int) + if initial_medoids.ndim != 1: + raise ValueError("initial_medoids must be a 1D array.") + if len(initial_medoids) != k: + raise ValueError(f"Number of initial medoids ({len(initial_medoids)}) must equal k ({k}).") + if len(np.unique(initial_medoids)) != k: + raise ValueError("initial_medoids must be unique.") + if np.any((initial_medoids < 0) | (initial_medoids >= n)): + raise ValueError(f"initial_medoids contain indices out of bounds (must be 0 to {n-1}).") + + n_init = 1 + + if method == 'alternate': + return _clusterKMedoidsAlternating(distance_matrix, k, initial_medoids, seed, max_iter, n_init) + elif method == 'pam': + return _clusterKMedoidsPAM(distance_matrix, k, initial_medoids, seed, max_iter, n_init) + elif method == 'sklearn': + if method_sklearn not in ['alternate', 'pam']: + raise ValueError(f"method_sklearn '{method_sklearn}' is not valid. Options: 'alternate', 'pam'") + return _clusterKMedoidsSklearn(distance_matrix, k, method_sklearn, seed, max_iter, n_init) + else: + raise ValueError(f"Method '{method}' is not valid. Options: 'alternate', 'pam', 'sklearn'") + + +def _cost(distance_matrix, medoids): + """Calculates the sum of distances from all points to their nearest medoid.""" + distances = distance_matrix[:, medoids] + closest_distance = np.min(distances, axis=1) + return np.sum(closest_distance) + + +def _clusterKMedoidsAlternating(distance_matrix, k, initial_medoids, seed, max_iter, n_init): + """K-Medoids clustering with an Alternating Medoids algorithm""" + + n = distance_matrix.shape[0] + rng = np.random.default_rng(seed) + + best_cost = np.inf + best_medoids = None + + for run in range(n_init): + if initial_medoids is None: + medoids = np.asarray(rng.choice(n, size = k, replace = False)) + else: + medoids = initial_medoids.copy() + + for iteration in range(max_iter): + distances_to_medoids = distance_matrix[:, medoids] + cluster_ids = np.argmin(distances_to_medoids, axis = 1) + + new_medoids = np.zeros_like(medoids) + + for i in range(len(medoids)): + cluster_members = np.where(cluster_ids == i)[0] + + # Handles empty clusters by retaining the old medoid + if len(cluster_members) == 0: + new_medoids[i] = medoids[i] + continue + + intra_cluster_distances = distance_matrix[np.ix_(cluster_members, cluster_members)] + sum_distances = intra_cluster_distances.sum(axis=1) + best_medoid_idx_in_cluster = np.argmin(sum_distances) + new_medoids[i] = cluster_members[best_medoid_idx_in_cluster] + + current_cost = _cost(distance_matrix, new_medoids) + + if np.array_equal(medoids, new_medoids): + break + + medoids = new_medoids + + if current_cost < best_cost: + best_cost = current_cost + best_medoids = medoids.copy() + + distances = distance_matrix[:, best_medoids] + cluster_ids = np.argmin(distances, axis=1) + 1 + + return cluster_ids, best_medoids, best_cost + + +def _clusterKMedoidsPAM(distance_matrix, k, initial_medoids, seed, max_iter, n_init): + """K-Medoids clustering with PAM (Partitioning Around Medoids)""" + + n = distance_matrix.shape[0] + rng = np.random.default_rng(seed) + + best_cost = np.inf + best_medoids = None + + for run in range(n_init): + if initial_medoids is None: + medoids = np.asarray(rng.choice(n, size=k, replace=False)) + else: + medoids = initial_medoids.copy() + + # Using a set for removal/addition operations + non_medoids = set(i for i in range(n) if i not in medoids) + current_cost = _cost(distance_matrix, medoids) + + for iteration in range(max_iter): + best_cost_swap = current_cost + best_swap = None + + for medoid_idx, old_medoid in enumerate(medoids): + for new_medoid in non_medoids: + candidate = medoids.copy() + candidate[medoid_idx] = new_medoid + candidate_cost = _cost(distance_matrix, candidate) + + if candidate_cost < best_cost_swap: + best_cost_swap = candidate_cost + best_swap = (medoid_idx, old_medoid, new_medoid) + + if best_swap is None: + break + + idx, old_medoid, new_medoid = best_swap + medoids[idx] = new_medoid + non_medoids.remove(new_medoid) + non_medoids.add(old_medoid) + current_cost = best_cost_swap + + if current_cost < best_cost: + best_cost = current_cost + best_medoids = medoids.copy() + + distances = distance_matrix[:, best_medoids] + cluster_ids = np.argmin(distances, axis=1) + 1 + + return cluster_ids, best_medoids, best_cost + + +def _clusterKMedoidsSklearn(distance_matrix, k, method_sklearn, seed, max_iter, n_init): + """K-Medoids clustering with sklearn_extra""" + try: + from sklearn_extra.cluster import KMedoids + except ImportError: + raise ImportError("The 'sklearn_extra' package is required for this K-Medoids approach. " + "Please install it using 'pip install scikit-learn-extra'.") + + + best_cost = np.inf + best_labels = None + best_medoids = None + + rng = np.random.default_rng(seed) + + for run in range(n_init): + kmedoids = KMedoids(n_clusters=k, metric='precomputed', method=method_sklearn, init='random', + max_iter=max_iter, random_state=int(rng.integers(0, 1000000))) + + kmedoids.fit(distance_matrix) + medoids = kmedoids.medoid_indices_ + + cost = np.sum(np.min(distance_matrix[:, medoids], axis=1)) + + if cost < best_cost: + best_cost = cost + best_labels = kmedoids.labels_ + 1 + best_medoids = medoids.copy() + + return best_labels, best_medoids, best_cost + + +def writeClusters(atoms, trajectory, distance_matrix, cluster_ids, write_dcd = True, + align ="protein and backbone", system = "system", tag ="cluster"): + """ + Aligns the trajectory once, then exports representative medoid structures as PDB files and + cluster-specific DCD trajectories. + Returns the name of all exported .pdb and .dcd files. + + Note: This function loads all aligned coordinates into memory. It is highly optimized + for speed, provided the trajectory fits within available system RAM. + + :arg atoms: reference structure used for the alignment. + :type atoms: :class:`prody.Atomic` + + :arg trajectory: trajectory containing the coordinate sets to align + :type trajectory: :class:`prody.Trajectory` + + :arg distance_matrix: square, symmetric matrix of pairwise distances between all objects + :type distance_matrix: :class:`numpy.ndarray` + + :arg cluster_ids: a one-dimensional array of the cluster IDs per element + :type cluster_ids: :class:`numpy.ndarray` + + :arg write_dcd: determines whether to save a .dcd file of the cluster frames + Default is ``True`` + :type write_dcd: bool + + :arg align: atom selection used to calculate the alignment transformation. + Must be a valid ProDy selection string. + Default is ``"protein and backbone"`` + :type align: str + + :arg system: the name of the system under investigation + Default is ``"system"`` + :type system: str + + :arg tag: the name of the clustering method + Default is ``"cluster"`` + :type tag: str + + :returns: list of exported filenames + :rtype: list[str] + + Example usage: + >>> pdb = prody.parsePDB("structure.pdb") + >>> dcd = prody.Trajectory("trajectory.dcd") + >>> _, aligned_coords = prody.alignTrajectory(pdb, dcd, select='resname IOA') + >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) + >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) + >>> prody. writeClusters(pdb, dcd, distance_matrix, cluster_ids, + system="type1_RUN23", tag="hier") + """ + + from prody import writeDCD + + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + cluster_ids = np.asarray(cluster_ids) + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") + if cluster_ids.size == 0: + raise ValueError("cluster_ids is empty.") + if len(cluster_ids) != trajectory.numFrames(): + raise ValueError("cluster_ids must have one entry per trajectory frame.") + + if not isinstance(write_dcd, bool): + raise TypeError(f"write_dcd must be a bool, but got {type(write_dcd).__name__}") + + clusters = np.unique(cluster_ids) + clusters = clusters[clusters > 0] + num_clusters = clusters.size + if num_clusters == 0: + raise ValueError("No clusters were found. All frames are labeled as noise") + + exported_files = [] + + # NOTE: Loads the entire aligned trajectory into memory + _, aligned_coords = alignTrajectory(atoms, trajectory, align=align, select="all") + + for cluster in clusters: + cluster_indices = getCluster(cluster_ids, cluster) + medoid = getClusterMedoid(distance_matrix, cluster_indices = cluster_indices) + + cluster_coords = aligned_coords[cluster_indices] + + cluster_atoms = atoms.copy() + cluster_atoms.setCoords(cluster_coords[0]) + + if len(cluster_coords) > 1: + cluster_atoms.addCoordset(cluster_coords[1:]) + + if write_dcd: + dcd_filename = f"{system}_{tag}_n{num_clusters}_cluster{cluster}.dcd" + writeDCD(dcd_filename, cluster_atoms) + exported_files.append(dcd_filename) + + medoid_atoms = atoms.copy() + medoid_atoms.setCoords(aligned_coords[medoid["global"]]) + + pdb_filename = f"{system}_{tag}_n{num_clusters}_cluster{cluster}_medoid.pdb" + writePDB(pdb_filename, medoid_atoms) + exported_files.append(pdb_filename) + + return exported_files + + + class Interactions(object): From eaf54d437862b1243e04af7539cbd40b7d732f2e Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Thu, 6 Aug 2026 15:42:59 +0200 Subject: [PATCH 02/10] New clustering function - DBSCAN, OPTICS for Reachability Plot --- prody/proteins/interactions.py | 397 ++++++++++++++++++++++++++++++++- 1 file changed, 395 insertions(+), 2 deletions(-) diff --git a/prody/proteins/interactions.py b/prody/proteins/interactions.py index 01aace7c0..1c7b4ebcc 100644 --- a/prody/proteins/interactions.py +++ b/prody/proteins/interactions.py @@ -59,7 +59,7 @@ 'calcClusterPopulations', 'getCluster', 'getClusterMedoid', 'calcClusterStatistics', 'calcAllClusterStatistics', 'showClusterStatisticsTable', 'showClusterRMSDComparison', 'clusterHierarchical', 'showDendrogram', - 'clusterKMedoids', 'writeClusters'] + 'clusterKMedoids', 'clusterDBSCAN', 'showReachabilityPlot', 'writeClusters'] def cleanNumbers(listContacts): @@ -4397,7 +4397,7 @@ def showRMSDfromReference(rmsd_array, title = "Frame-to-Reference RMSD", :type ax: :class:`matplotlib.axes.Axes` :returns: the Matplotlib axes containing the plot. - :rtype: :class:`matplotlib.axes.Axes`, None + :rtype: :class:`matplotlib.axes.Axes` Example usage: >>> import matplotlib.pyplot as plt @@ -4904,6 +4904,10 @@ def calcAllClusterStatistics(distance_matrix, cluster_ids): all_stats = [] for c in clusters: + # Ignore noise + if c <= 0: + continue + cluster_stats = calcClusterStatistics(distance_matrix, cluster_ids = cluster_ids, cluster_number = int(c)) all_stats.append(cluster_stats) @@ -5591,6 +5595,395 @@ def _clusterKMedoidsSklearn(distance_matrix, k, method_sklearn, seed, max_iter, return best_labels, best_medoids, best_cost +def clusterDBSCAN(distance_matrix, eps = None, minPts = None, method='custom'): + """ + Performs DBSCAN clustering using various algorithms. + + This function acts as a facade, routing the clustering task to the specified backend + ('custom' or 'sklearn'). 'custom' uses a built-in implementation without additional + dependencies beyond NumPy. 'sklearn' needs the sklearn.cluster module to be installed. + + Note that the resulting cluster IDs are 1-based (starting at 1, not 0). + Noise points are labeled -1. + + + :arg distance_matrix: square, symmetric matrix of pairwise distances between all objects + :type distance_matrix: :class:`numpy.ndarray` + + :arg eps: the "radius" of the neighborhood within which we count neighbors + Default is ``None`` and automatically the median of the distances + in the distance matrix is used. + Ideally use :func:`showReachabilityPlot` to determine + manually the most suitable eps. + :type eps: float + + :arg minPts: the minimum number of neighbors required for a point to be considered + a core point + Default is ``None`` and automatically the 5% of the total objects, + or for less than 20 2 is used. + Ideally choose manually the most suitable minPts + :type minPts: int + + :arg method: the clustering algorithm to use + Options are 'custom' and 'sklearn' + Default is 'custom' because it needs no module installation + :type method: str + + :returns: a tuple of: + * a one-dimensional array containing the cluster ID for each object. + IDs are 1-indexed and noise corresponds to -1. + * a one-dimensional array of the frame indices corresponding to noise + :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) + + Example usage: + >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) + >>> cluster_ids, _ = prody.clusterDBSCAN(distance_matrix, eps = 1.8, minPts = 30) + """ + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + if eps is None: + eps = _calcAutoEps(distance_matrix) + _validateEps(eps, distance_matrix) + + if minPts is None: + minPts = _calcAutoMinPts(distance_matrix) + _validateMinPts(minPts, distance_matrix) + + if method == 'custom': + return _clusterDBSCANCustom(distance_matrix, eps, minPts) + elif method == 'sklearn': + return _clusterDBSCANSklearn(distance_matrix, eps, minPts) + else: + raise ValueError(f"method can be either 'custom' or 'sklearn', but got {method}.") + + +def _clusterDBSCANCustom(distance_matrix, eps, minPts): + """DBSCAN clustering with custom algorithm""" + total_points = distance_matrix.shape[0] + + labels = np.zeros(total_points, dtype = int) + cluster_id = 0 + + for p in range(total_points): + + if labels[p] != 0: + continue + + neighbors, = np.where(distance_matrix[p] <= eps) + + if len(neighbors) < minPts: + labels[p] = -1 + else: + cluster_id += 1 + labels[p] = cluster_id + + candidate_set = [n for n in neighbors if n != p] + while candidate_set: + q = candidate_set.pop() + + if labels[q] == -1: + labels[q] = cluster_id + + if labels[q] != 0: + continue + + labels[q] = cluster_id + + q_neighbors, = np.where(distance_matrix[q] <= eps) + if len(q_neighbors) >= minPts: + for n in q_neighbors: + if labels[n] == 0: + candidate_set.append(n) + elif labels[n] == -1: + labels[n] = cluster_id + + noise_frames, = np.where(labels == -1) + + return labels, noise_frames + + +def _clusterDBSCANSklearn(distance_matrix, eps, minPts): + """DBSCAN clustering with sklearn""" + + try: + from sklearn.cluster import DBSCAN + except ImportError: + raise ImportError("The 'sklearn' package is required for this DBSCAN approach. " + "Please install it using 'pip install scikit-learn'.") + + dbscan = DBSCAN(eps=eps, min_samples = minPts, metric = "precomputed") + labels = dbscan.fit_predict(distance_matrix) + + cluster_ids = np.copy(labels) + cluster_ids[cluster_ids >= 0] += 1 + + noise_frames, = np.where(labels == -1) + + return cluster_ids, noise_frames + + +def showReachabilityPlot(distance_matrix, minPts = None, method = 'custom', eps = None, + title = "OPTICS Reachability Plot", xlabel = "Frames (Sorted by OPTICS)", + ylabel = "Reachability Distance [Å]", label = None, + color = "#36454F", lw = 1.5, fill = True, ax = None): + """ + Plots the reachability plot using a simplified OPTICS algorithm. + The reachability plot should be used to determine the most suitable eps + parameter for DBSCAN. + + This function acts as a facade, routing the ordering task to the specified backend + ('custom' or 'sklearn'). 'custom' uses a built-in implementation without additional + dependencies beyond NumPy. 'sklearn' needs the sklearn.cluster module to be installed. + + + :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg minPts: the minimum number of neighbors required for a point to be considered + a core point + Default is ``None`` and automatically the 5% of the total objects, + or for less than 20 2 is used. + Ideally choose manually the most suitable minPts + :type minPts: int + + :arg method: the OPTICS algorithm to use + Options are 'custom' and 'sklearn'. + Default is 'custom' because it needs no module installation + :type method: str + + :arg eps: the "radius" of the neighborhood within which we count neighbors. + It is drawn as a horizontal line to help visualize the DBSCAN + ``eps`` threshold. + If `'auto'` the median of the pairwise distances is used. + Default is ``None``. + :type eps: float, str + + :arg title: title of the generated plot. + Default is ``"OPTICS Reachability Plot"`` + :type title: str + + :arg xlabel: Frame count label for the x-axis + Default is ``"Frames (Sorted by OPTICS)"`` + :type xlabel: str + + :arg ylabel: Reachability distance label for the y-axis + Default is ``"Reachability Distance [Å]"`` + :type ylabel: str + + :arg label: the label for the plot. + Default is ``None`` + :type label: str + + :arg color: the color of the plot. + Default is ``"#36454F"`` + :type color: str + + :arg lw: the width of the line. + Default is ``1.5`` + :type lw: float + + :arg fill: whether to fill the area beneath the curve and if applicable + the valleys below the eps threshold. + :type fill: bool + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :returns: the Matplotlib axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> plt.figure(figsize = (8, 6)) + >>> showReachabilityPlot(distance_matrix, minPts = 20, eps = 1.8) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + if minPts is None: + minPts = _calcAutoMinPts(distance_matrix) + _validateMinPts(minPts, distance_matrix) + + if not isinstance(fill, bool): + raise TypeError(f"fill must be a bool, but got {type(fill).__name__}") + + if method == 'custom': + reachability, ordering = _orderOPTICSCustom(distance_matrix, minPts) + elif method == 'sklearn': + reachability, ordering = _orderOPTICSSklearn(distance_matrix, minPts) + else: + raise ValueError(f"method must be either 'custom' or 'sklearn', but got {method}") + + if ax is None: + ax = plt.gca() + + ax.set_title(title) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + + y = reachability[ordering] + x = np.arange(len(y)) + + ax.plot(x, y, color = color, lw = lw, label = label) + + if fill: + ax.fill_between(x, 0, y, color = 'black', alpha = 0.4) + + if eps == 'auto': + eps = _calcAutoEps(distance_matrix) + + if eps is not None: + _validateEps(eps, distance_matrix) + ax.axhline(y = eps, color = 'black', linestyle = '--', linewidth = lw) + + if fill: + below_eps_mask = (y <= eps) + segments = [] + start = None + + for i, below in enumerate(below_eps_mask): + if below and start is None: + start = i + elif not below and start is not None: + segments.append((start, i)) + start = None + + if start is not None: + segments.append((start, len(below_eps_mask))) + + colors = plt.cm.tab10(np.linspace(0, 1, len(segments))) + + for (start, end), valley_color in zip(segments, colors): + ax.fill_between(x[start:end], y[start:end], eps, color=valley_color) + + if label is not None: + ax.legend() + + ax.grid(axis = 'y', linestyle = '--', linewidth = 0.8, alpha = 0.3) + + return ax + + +def _orderOPTICSCustom(distance_matrix, minPts): + """OPTICS algorithm with built-in modules and NumPy""" + + import heapq + + total_points = distance_matrix.shape[0] + + sorted_distances = np.sort(distance_matrix, axis = 1) + core_distances = sorted_distances[:, minPts - 1] + + reachability = np.full(total_points, np.inf) + processed = np.zeros(total_points, dtype = bool) + ordering = [] + + def _updateSeeds(idx): + new_reaches = np.maximum(core_distances[idx], distance_matrix[idx, :]) + update_mask = (~processed) & (new_reaches < reachability) + points_to_update, = np.where(update_mask) + reachability[update_mask] = new_reaches[update_mask] + + for j in points_to_update: + heapq.heappush(seeds, (reachability[j], j)) + + for i in range(total_points): + if processed[i]: + continue + + processed[i] = True + ordering.append(i) + seeds = [] + + _updateSeeds(i) + + while seeds: + current_reach, q = heapq.heappop(seeds) + + if processed[q]: + continue + + processed[q] = True + ordering.append(q) + + _updateSeeds(q) + + return reachability, np.array(ordering) + + +def _orderOPTICSSklearn(distance_matrix, minPts): + """OPTICS algorithm with sklearn""" + + try: + from sklearn.cluster import OPTICS + except ImportError: + raise ImportError("The 'sklearn' package is required for this OPTICS approach. " + "Please install it using 'pip install scikit-learn'.") + + optics = OPTICS(min_samples = minPts, metric = 'precomputed') + optics.fit(distance_matrix) + reachability = optics.reachability_ + ordering = optics.ordering_ + + return reachability, ordering + + +def _calcAutoEps(distance_matrix): + """Automatically determine the DBSCAN eps parameter""" + + import warnings + + pairwise_distances = distance_matrix[np.triu_indices_from(distance_matrix, k = 1)] + eps = float(np.median(pairwise_distances)) + warnings.warn(f"No eps provided. Automatically chosen at {eps:.3f}. Ideally provide your own value.") + return eps + + +def _validateEps(eps, distance_matrix): + """Validate the DBSCAN eps parameter""" + + import warnings + + maxDistance = np.max(distance_matrix) + + if not isinstance(eps, (float, int, np.floating, np.integer)): + raise TypeError(f"eps must be a numeric value, but got {type(eps).__name__}") + + if eps <= 0: + raise ValueError("eps must be positive") + elif eps > maxDistance: + warnings.warn(f"eps ({eps}) is greater than the maximum pairwise distance ({maxDistance:.3f}).\n" + "All frames will be clustered together with no noise.") + + +def _calcAutoMinPts(distance_matrix): + """Automatically determine the minPts parameter""" + + import warnings + + total_points = distance_matrix.shape[0] + minPts = max(2, int(total_points // 20)) + warnings.warn(f"No minPts provided. Automatically chosen at {minPts}. Ideally provide your own value.") + return minPts + + +def _validateMinPts(minPts, distance_matrix): + """Validate minPts parameter""" + total_points = distance_matrix.shape[0] + + if not isinstance(minPts, (int, np.integer)): + raise TypeError(f"minPts must be a positive integer, but got {type(minPts).__name__}") + + if minPts <= 0 or minPts > total_points: + raise ValueError(f"minPts must be between 1 and {total_points}.") + + def writeClusters(atoms, trajectory, distance_matrix, cluster_ids, write_dcd = True, align ="protein and backbone", system = "system", tag ="cluster"): """ From b91f0e7eae4385c92b00bad8bcacc0fd493ab9da Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Mon, 10 Aug 2026 15:13:01 +0200 Subject: [PATCH 03/10] Created rmsd_clustering module - Relocated functions - Cleaned up Redundant Code --- prody/dynamics/__init__.py | 4 + prody/dynamics/plotting.py | 306 ++++- prody/dynamics/rmsd_clustering.py | 1491 +++++++++++++++++++++++ prody/proteins/interactions.py | 1854 +---------------------------- prody/utilities/catchall.py | 108 +- 5 files changed, 1908 insertions(+), 1855 deletions(-) create mode 100644 prody/dynamics/rmsd_clustering.py diff --git a/prody/dynamics/__init__.py b/prody/dynamics/__init__.py index 9b07d994b..e0da0e780 100644 --- a/prody/dynamics/__init__.py +++ b/prody/dynamics/__init__.py @@ -392,3 +392,7 @@ from . import anmd from .anmd import * __all__.extend(anmd.__all__) + +from . import rmsd_clustering +from .rmsd_clustering import * +__all__.extend(rmsd_clustering.__all__) \ No newline at end of file diff --git a/prody/dynamics/plotting.py b/prody/dynamics/plotting.py index 51687e0a0..b52535410 100644 --- a/prody/dynamics/plotting.py +++ b/prody/dynamics/plotting.py @@ -13,6 +13,7 @@ from prody import LOGGER, SETTINGS, PY3K from prody.utilities import showFigure, addEnds, showMatrix from prody.atomic import AtomGroup, Selection, Atomic, sliceAtoms, sliceAtomicData +from prody.measure import calcRMSD from .nma import NMA from .gnm import GNMBase, GNM @@ -25,6 +26,8 @@ from .compare import calcOverlap from .lda import LDA from .logistic import LRA +from .rmsd_clustering import _validateDistanceMatrix, calcPairwiseRMSD + __all__ = ['showContactMap', 'showCrossCorr', 'showCovarianceMatrix', 'showCumulOverlap', 'showFractVars', @@ -37,7 +40,8 @@ 'showPairDeformationDist','showMeanMechStiff', 'showPerturbResponse', 'showTree', 'showTree_networkx', 'showAtomicMatrix', 'pimshow', 'showAtomicLines', 'pplot', - 'showDomainBar', 'showAtomicBars', 'showSelectionMatrix'] + 'showDomainBar', 'showAtomicBars', 'showSelectionMatrix', + 'showRMSDEvolution', 'showPairwiseRMSDHeatmap', 'showClusterRMSDComparison'] def showEllipsoid(modes, onto=None, n_std=2, scale=1., *args, **kwargs): @@ -2703,3 +2707,303 @@ def showSelectionMatrix(matrix, atoms, selstr_x=None, selstr_y=None, **kwargs): _, atoms_y = sliceAtoms(atoms, selstr_y) return showAtomicMatrix(matrix, atoms=[atoms_x, atoms_y], **kwargs) + + +def showRMSDEvolution(rmsd_array=None, ref_coords=None, aligned_coords=None, *args, **kwargs): + """ + Plots RMSD of an aligned trajectory from a reference structure. + It accepts either a pre-calculated one-dimensional array of the RMSD over the frames, or + the reference coordinates and the aligned coordinates of the trajectory and performs the + RMSD calculation itself. + + + :arg rmsd_array: one-dimensional array containing the RMSD values [in Å] + for each trajectory frame relative to the reference structure. + :type rmsd_array: :class:`numpy.ndarray` or atomic object + + :arg ref_coords: coordinates of the reference structure + either as an array with shape ``(n_atoms, 3)`` or ``(1, n_atoms, 3)``, + or as an atomic object with a ``getCoords`` method. + :type ref_coords: :class:`numpy.ndarray` + + :arg aligned_coords: coordinates of the aligned trajectory. + :type aligned_coords: :class:`numpy.ndarray` + + :arg *args: positional arguments passed to Matplotlib's ``plot``. + :type *args: tuple + + :arg title: title of the plot. + Default is ``"Frame-to-Reference RMSD"``. + :type title: str + + :arg xlabel: label for the x-axis. + Default is ``"# of Frame"``. + :type xlabel: str + + :arg ylabel: label for the y-axis. + Default is ``"RMSD from Reference [Å]"``. + :type ylabel: str + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :arg **kwargs: keyword arguments passed to Matplotlib's ``plot`` function, + excluding ``title`` and ``ax``. + :type **kwargs: dict + + :return: axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> ref_coords, aligned_coords = alignTrajectory(pdb, dcd, select='resname IOA') + >>> plt.figure() + >>> showRMSDEvolution(ref_coords=ref_coords, aligned_coords=aligned_coords) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + + + has_array = rmsd_array is not None + has_coords = ref_coords is not None or aligned_coords is not None + + if has_array == has_coords: + raise ValueError("Provide either 'rmsd_array' or both 'ref_coords' and 'aligned_coords'.") + + if not has_array: + if hasattr(ref_coords, 'getCoords'): + ref_coords = ref_coords.getCoords() + ref_coords = np.asarray(ref_coords) + + if ref_coords.ndim == 3: + if ref_coords.shape[0] != 1: + raise ValueError("reference coordinates must have shape (n_atoms, 3) or (1, n_atoms, 3).") + ref_coords = ref_coords[0] + + if ref_coords.ndim != 2 or ref_coords.shape[1] != 3: + raise ValueError("reference coordinates must have shape (n_atoms, 3) or (1, n_atoms, 3).") + + aligned_coords = np.asarray(aligned_coords) + + if aligned_coords.ndim != 3 or aligned_coords.shape[2] != 3: + raise ValueError(f"aligned_coords must have shape (n_frames, n_atoms, 3), but got {aligned_coords.shape}.") + + if ref_coords.shape != aligned_coords.shape[1:]: + raise ValueError(f"Incompatible shapes: reference is {ref_coords.shape}, but aligned frames have {aligned_coords.shape[1:]}.") + + rmsd_array = calcRMSD(ref_coords, target=aligned_coords) + + rmsd_array = np.asarray(rmsd_array) + + if rmsd_array.ndim != 1: + raise ValueError(f"rmsd_array must be one-dimensional, but got shape {rmsd_array.shape}") + if rmsd_array.size == 0: + raise ValueError("rmsd_array cannot be empty.") + + ax = kwargs.pop('ax', None) + if ax is None: + ax = plt.gca() + + xlabel = kwargs.pop('xlabel', "# Frame") + ylabel = kwargs.pop('ylabel', "RMSD from Reference [Å]") + title = kwargs.pop('title', "Frame-to-Reference RMSD") + label = kwargs.get('label', None) + + if 'color' not in kwargs and 'c' not in kwargs and not args: + kwargs['color'] = "#36454F" + + if 'lw' not in kwargs and 'linewidth' not in kwargs: + kwargs['lw'] = 0.75 + + ax.plot(rmsd_array, *args, **kwargs) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.set_title(title) + ax.grid(alpha=0.3) + + if label is not None: + ax.legend() + + return ax + + +def showPairwiseRMSDHeatmap(distance_matrix=None, aligned_coords=None, *args, **kwargs): + """ + Plots the heatmap of the frame-to-frame pairwise RMSDs using the distance matrix. + It accepts either a pre-computed pairwise distance matrix, or + the aligned coordinates of the frames and calculates the matrix itself. + + + :arg distance_matrix: two-dimensional array constituting the pairwise distance matrix. + Typically generated using :func:`calcPairwiseRMSD`. + :type distance_matrix: :class:`numpy.ndarray` + + :arg aligned_coords: coordinates of the aligned trajectory. + :type aligned_coords: :class:`numpy.ndarray` + + :arg *args: positional arguments passed to ProDy's ``showMatrix`` function. + :type *args: tuple + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :arg **kwargs: keyword arguments passed directly to ProDy's ``showMatrix`` function + :type **kwargs: dict + + :returns: the Matplotlib axes containing the heatmap. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> _, aligned_coords = alignTrajectory(pdb, dcd, select='resname IOA') + >>> distance_matrix = calcPairwiseRMSD(aligned_coords) + >>> plt.figure() + >>> showPairwiseRMSDHeatmap(distance_matrix) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + + + has_matrix = distance_matrix is not None + has_coords = aligned_coords is not None + + if has_matrix == has_coords: + raise ValueError("Provide exactly one of 'distance_matrix' or 'aligned_coords'.") + + if has_coords: + distance_matrix = calcPairwiseRMSD(aligned_coords) + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + ax = kwargs.pop('ax', None) + if ax is not None: + plt.sca(ax) + else: + ax = plt.gca() + + kwargs.setdefault('cmap', 'viridis') + kwargs.setdefault('origin', 'upper') + showMatrix(distance_matrix, *args, **kwargs) + + title = kwargs.pop('title', 'Pairwise RMSD Distance Matrix') + ax.set_title(title) + + return ax + + +def showClusterRMSDComparison(all_stats, *args, **kwargs): + """ + Overlays the internal RMSD distributions of all clusters onto a single plot + on the current axis. + + Expects the list of cluster statistics dictionaries generated by :func:`calcAllClusterStatistics`. + + + :arg all_stats: list of cluster statistics dictionaries. Each dictionary + must contain the keys ``"cluster"`` and ``"distances"``. + :type all_stats: list of dict, or dict + + :arg *args: positional arguments passed directly to Seaborn's ``histplot`` function. + :type *args: tuple + + :arg title: the title of the generated plot. + Default is ``'Intra-Cluster RMSD Distributions'`` + :type title: str + + :arg xlabel: the label for the x-axis. + Default is ``'RMSD to Medoid [Å]'`` + :type xlabel: str + + :arg ylabel: the label for the y-axis. + Default is ``'Frequency'`` + :type ylabel: str + + :arg grid: whether to display horizontal grid lines. + Default is ``True``. + :type grid: bool + + :arg label: optional label prefix applied to each cluster. + For example, ``label='Run 1'`` produces labels such as ``'Run 1 - Cluster 1'``. + :type label: str + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :arg **kwargs: keyword arguments passed directly to Seaborn's ``histplot`` function + (e.g., ``bins``, ``kde``, ``element``, ``stat``, ``alpha``, ``lw``). + :type **kwargs: dict + + :returns: the Matplotlib axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> all_stats = calcAllClusterStatistics(distance_matrix, cluster_ids) + >>> plt.figure() + >>> showClusterRMSDComparison(all_stats, label="Simulation A") + >>> plt.show() + """ + + import matplotlib.pyplot as plt + try: + import seaborn as sns + except ImportError: + raise ImportError("The 'seaborn' package is required to display the histogram." + "\nPlease install it using 'pip install seaborn'." + "\nAlternatively, use standard matplotlib.pyplot.hist() for basic plots.") + + + if isinstance(all_stats, dict): + all_stats = [all_stats] + + if not all_stats: + raise ValueError("all_stats cannot be empty.") + + ax = kwargs.pop('ax', None) + if ax is None: + ax = plt.gca() + + title = kwargs.pop('title', 'Intra-Cluster RMSD Distributions') + xlabel = kwargs.pop('xlabel', 'RMSD to Medoid [Å]') + ylabel = kwargs.pop('ylabel', 'Frequency') + grid = kwargs.pop('grid', True) + + user_label = kwargs.pop('label', None) + + # Seaborn Defaults + kwargs.setdefault('bins', 50) + kwargs.setdefault('element', 'step') + kwargs.setdefault('stat', 'count') + kwargs.setdefault('kde', True) + kwargs.setdefault('alpha', 0.5) + + if 'lw' in kwargs: + kwargs['linewidth'] = kwargs.pop('lw') + else: + kwargs.setdefault('linewidth', 1.5) + + for cluster_stats in all_stats: + if "cluster" not in cluster_stats or "distances" not in cluster_stats: + raise ValueError("Each cluster statistics dictionary must contain 'cluster' and 'distances'.") + + cluster_label = f"Cluster {cluster_stats['cluster']}" + if user_label: + cluster_label = f"{user_label} - {cluster_label}" + + sns.histplot(cluster_stats["distances"], *args, label=cluster_label, ax=ax, **kwargs) + + ax.set_title(title) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + + ax.legend() + + if grid: + ax.grid(axis='y', alpha=0.3) + + return ax \ No newline at end of file diff --git a/prody/dynamics/rmsd_clustering.py b/prody/dynamics/rmsd_clustering.py new file mode 100644 index 000000000..192fcc738 --- /dev/null +++ b/prody/dynamics/rmsd_clustering.py @@ -0,0 +1,1491 @@ +# -*- coding: utf-8 -*- + +""" +This module defines functions for clustering. It provides tools for statistical analysis of the results +""" + +__author__ = 'Konstantinos Theodosiadis' +__credits__ = ['James Krieger', 'Karolina Mikulska-Ruminska', 'Konstantinos Theodosiadis'] +__email__ = ['karolamik@fizyka.umk.pl', 'jamesmkrieger@gmail.com', 'konst.theodosiadis@gmail.com'] + + +import numpy as np +from prody.atomic import AtomGroup, Atom, Atomic, Selection, Select +from prody.utilities import getCoords +from prody.proteins import writePDB +from prody.trajectory import writeDCD +from prody.measure import calcTransformation + +__all__ = ['alignTrajectory', 'calcPairwiseRMSD', 'calcClusterPopulations', 'getCluster', 'getClusterMedoid', + 'calcClusterStatistics', 'calcAllClusterStatistics', 'showClusterStatisticsTable', 'writeClusters', + 'clusterHierarchical', 'showDendrogram', 'clusterKMedoids', 'clusterDBSCAN', 'showReachabilityPlot'] + + +def alignTrajectory(atoms, trajectory, align='protein and backbone', select='all'): + """ + Aligns each trajectory frame to the reference structure and returns a tuple of the + reference coordinates and the aligned coordinates of the selected atoms. + + The trajectory frames are aligned to the reference structure using the atoms specified by + ``align``. After alignment, the coordinates of the atoms specified by ``select`` are + extracted for every frame. + + + :arg atoms: reference structure used for the alignment. + :type atoms: :class:`prody.Atomic` + + :arg trajectory: trajectory containing the coordinate sets to align + :type trajectory: :class:`prody.Trajectory` + + :arg align: atom selection used to calculate the alignment transformation. + Must be a valid ProDy selection string. + Default is ``"protein and backbone"`` + :type align: str + + :arg select: atom selection whose coordinates are returned. + Must be a valid ProDy selection string. + Default is ``"all"`` + :type select: str + + :returns: a tuple containing: + * ref_coords (numpy.ndarray): coordinates of the selected atoms in the reference structure. + * aligned_coords (numpy.ndarray): aligned coordinates of the selected atoms for every + trajectory frame with shape ``(n_frames, n_atoms, 3)``. + :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) + + Example usage: + >>> pdb = parsePDB("structure.pdb") + >>> dcd = Trajectory("trajectory.dcd") + >>> ref_coords, aligned_coords = alignTrajectory(pdb, dcd, select='resname IOA') + """ + + if trajectory.numAtoms() != atoms.numAtoms(): + raise ValueError("Trajectory atoms count does not match structure atoms count.") + + # Save original coordinates to restore in the end + orig_coords = atoms.getCoords().copy() + + # Atoms we want to align + atom_align = atoms.select(align) + if atom_align is None: + raise ValueError(f"No atoms match '{align}' in the structure.") + ref_align = atom_align.copy() + + # Atoms we want to return + atom_select = atoms.select(select) + if atom_select is None: + raise ValueError(f"No atoms match '{select}' in the structure.") + ref_coords = atom_select.getCoords().copy() + + trajectory.link(atoms) # linking trajectory to update coordinates frame-by-frame + + n_frames = trajectory.numFrames() + n_atoms = atom_select.numAtoms() + trajectory.reset() + aligned_coords = np.zeros((n_frames, n_atoms, 3)) + + try: + for i, frame in enumerate(trajectory): + trans = calcTransformation(atom_align, ref_align) + trans.apply(atom_select) + aligned_coords[i] = atom_select.getCoords() + finally: + atoms.setCoords(orig_coords) + trajectory.reset() + + return ref_coords, aligned_coords + + +def calcPairwiseRMSD(aligned_coords): + """ + Calculates the frame-to-frame pairwise RMSD matrix using aligned structures + and returns a symmetric distance matrix. + Uses a vectorized approach for better efficiency. + + + :arg aligned_coords: aligned coordinates with shape ``(n_frames, n_atoms, 3)``. + Recommended to generate them using :func:`alignTrajectory`. + :type aligned_coords: :class:`numpy.ndarray` + + :returns: symmetric matrix containing the pairwise RMSD between all frames with shape + ``(n_frames, n_frames)``. + :rtype: :class:`numpy.ndarray` + + Example usage: + >>> ref_coords, aligned_coords = alignTrajectory(pdb, dcd, select='resname IOA') + >>> distance_matrix = calcPairwiseRMSD(aligned_coords) + """ + + from scipy.spatial.distance import cdist + + + aligned_coords = np.asarray(aligned_coords) + + if aligned_coords.ndim != 3 or aligned_coords.shape[2] != 3: + raise ValueError(f"aligned_coords must have shape (n_frames, n_atoms, 3), but got {aligned_coords.shape}.") + + if aligned_coords.shape[0] == 0: + raise ValueError("aligned_coords contains no frames.") + + if aligned_coords.shape[1] == 0: + raise ValueError("aligned_coords contains no atoms.") + + + n_frames, n_atoms, _ = aligned_coords.shape + + # Flatten frame coordinates (n_frames, n_atoms, 3) -> (n_frames, n_atoms * 3) + coords_flat = aligned_coords.reshape(n_frames, 3 * n_atoms) + + euclidean_dists = cdist(coords_flat, coords_flat, metric = 'euclidean') + distance_matrix = euclidean_dists / np.sqrt(n_atoms) + + return distance_matrix + + +def calcClusterPopulations(cluster_ids): + """ + Uses the cluster IDs array, assigning each frame to a cluster, to calculate the population + of each cluster and its corresponding percentage. + + + :arg cluster_ids: a one-dimensional array matching each frame to a cluster. + :type cluster_ids: :class:`numpy.ndarray` + + :returns: a dictionary mapping each cluster ID to its population statistics. + Each value is a dictionary with the keys ``"count"`` and ``"pct"``. + :rtype: dict + + Example usage: + >>> cluster_ids, _ = clusterHierarchical(distance_matrix) + >>> populations = calcClusterPopulations(cluster_ids) + >>> print(f"In cluster 1: {populations[1]['pct']:.2f}% of the objects") + """ + + cluster_ids = np.asarray(cluster_ids) + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}") + if cluster_ids.size == 0: + raise ValueError("cluster_ids is empty.") + + total_frames = len(cluster_ids) + clusters, frequencies = np.unique(cluster_ids, return_counts=True) + + population_data = {} + for c, n in zip(clusters, frequencies): + population_data[int(c)] = {'count': int(n), 'pct': (n / total_frames) * 100} + + return population_data + + +def getCluster(cluster_ids, cluster_number): + """ + Returns the indices of the frames grouped in a specific cluster. + + + :arg cluster_ids: a one-dimensional array assigning each frame to a cluster. + :type cluster_ids: :class:`numpy.ndarray` + + :arg cluster_number: the ID of the cluster whose members we want to return. + :type cluster_number: int + + :returns: an array of the trajectory indices assigned to the specified cluster. + :rtype: :class:`numpy.ndarray` + + Example usage: + >>> cluster_ids, _ = clusterHierarchical(distance_matrix) + >>> cluster1_indices = getCluster(cluster_ids, 1) + """ + + cluster_ids = np.asarray(cluster_ids) + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") + if cluster_ids.size == 0: + raise ValueError("cluster_ids is empty.") + + if isinstance(cluster_number, bool) or not isinstance(cluster_number, (int, np.integer)): + raise TypeError(f"cluster_number must be an integer, but got {type(cluster_number).__name__}.") + + cluster_indices = np.where(cluster_ids == cluster_number)[0] + + # Fail-safe + if cluster_indices.size == 0: + raise ValueError(f'No frames belong to Cluster {cluster_number}') + + return cluster_indices + + +def getClusterMedoid(distance_matrix, cluster_indices=None, cluster_ids=None, cluster_number=None): + """ + Determines the medoid (representative) of the cluster from the minimum total pairwise + distance to all other cluster members. + Returns a dictionary containing the medoid index within the cluster ("local") and in + the full trajectory ("global"). + + The cluster can be specified either by: + 1. cluster_indices (recommended and can be generated by :func:`getCluster`) or + 2. cluster_ids and cluster_number + + + :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg cluster_indices: a one-dimensional array of the indices of the cluster members + Default is ``None`` + :type cluster_indices: :class:`numpy.ndarray` + + :arg cluster_ids: a one-dimensional array of the cluster IDs per element + Default is ``None`` + :type cluster_ids: :class:`numpy.ndarray` + + :arg cluster_number: the ID (number) of the cluster we want to investigate + Default is ``None`` + :type cluster_number: int + + :returns: a dictionary of the medoid index, with keys "local" for the index within the + cluster and "global" for the index within the full trajectory + :rtype: dict + + Example usage: + >>> cluster_ids, _ = clusterHierarchical(distance_matrix) + >>> cluster1_indices = getCluster(cluster_ids, 1) + >>> cluster1_medoid = getClusterMedoid(distance_matrix, cluster1_indices) + >>> print(f"Medoid of cluster 1 is # {cluster1_medoid['local']} in the cluster.") + """ + + cluster_indices = _resolveClusterIndices(cluster_indices=cluster_indices, + cluster_ids=cluster_ids, + cluster_number=cluster_number) + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + n_frames = distance_matrix.shape[0] + if np.any((cluster_indices < 0) | (cluster_indices >= n_frames)): + raise ValueError("cluster_indices contain indices outside the bounds of the distance matrix.") + + cluster_distance_matrix = distance_matrix[np.ix_(cluster_indices, cluster_indices)] + + sum_dist = cluster_distance_matrix.sum(axis=1) + medoid_local = int(np.argmin(sum_dist)) + medoid_global = int(cluster_indices[medoid_local]) + + return {"local": medoid_local, "global": medoid_global} + + +def calcClusterStatistics(distance_matrix, cluster_indices=None, cluster_medoid=None, + cluster_ids=None, cluster_number=None): + """ + Calculates descriptive statistics for a single cluster. + Returns a dictionary with the statistics quantity as a key. + + The cluster can be specified either by: + 1. cluster_indices (recommended and can be generated by :func:`getCluster`) or + 2. cluster_ids and cluster_number + + + :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg cluster_indices: a one-dimensional array of the indices of the cluster members + Default is ``None`` + :type cluster_indices: :class:`numpy.ndarray` + + :arg cluster_medoid: the dictionary of the cluster medoid. + Expects the 'cluster_medoid' dictionary generated by the function + :func:`getClusterMedoid`. + Default is ``None`` + :type cluster_medoid: dict + + :arg cluster_ids: a one-dimensional array of the cluster IDs per element + Default is ``None`` + :type cluster_ids: :class:`numpy.ndarray` + + :arg cluster_number: the ID (number) of the cluster we want to investigate + Default is ``None`` + :type cluster_number: int + + :returns: a dictionary of statistics quantities with their names as keys + :rtype: dict + + Example usage: + >>> cluster_ids, _ = clusterHierarchical(distance_matrix) + >>> cluster1_indices = getCluster(cluster_ids, 1) + >>> cluster_stats1 = calcClusterStatistics(distance_matrix, cluster1_indices) + >>> print(f"Mean distance from cluster Medoid: {cluster_stats1['mean']:.2f} ± {cluster_stats1['std']:.2f} [Å]") + """ + + cluster_indices = _resolveClusterIndices(cluster_indices, cluster_ids, cluster_number) + population = len(cluster_indices) + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + if cluster_medoid is None: + cluster_medoid = getClusterMedoid(distance_matrix, cluster_indices=cluster_indices) + else: + if not isinstance(cluster_medoid, dict): + raise TypeError("cluster_medoid must be a dictionary with keys 'global' and 'local'.") + + required_keys = {"global", "local"} + if not required_keys.issubset(cluster_medoid.keys()): + missing = required_keys - cluster_medoid.keys() + raise ValueError(f"Expected keys are 'global' and 'local'.\ncluster_medoid is missing required key(s): {missing}.") + + if cluster_medoid["global"] not in cluster_indices: + raise ValueError("Index not found.\ncluster_medoid does not belong to the specified cluster.") + + if not (0 <= cluster_medoid["local"] < population): + raise ValueError("Local index is out of bounds.") + + if cluster_indices[cluster_medoid["local"]] != cluster_medoid["global"]: + raise ValueError("cluster_medoid['local'] and cluster_medoid['global'] are inconsistent.") + + distance_to_medoid = distance_matrix[cluster_indices, cluster_medoid["global"]] + total_frames = distance_matrix.shape[0] + + stats = { + "population" : population, + "pct" : (population / total_frames) * 100, + "medoid_global" : cluster_medoid["global"], + "medoid_local" : cluster_medoid["local"], + "distances" : distance_to_medoid, + "mean" : np.mean(distance_to_medoid), + "std" : np.std(distance_to_medoid), + "median" : np.median(distance_to_medoid), + "iqr" : np.percentile(distance_to_medoid, 75) - np.percentile(distance_to_medoid, 25), + "p95" : np.percentile(distance_to_medoid, 95), + "max" : np.max(distance_to_medoid) + } + + if cluster_number is not None: + stats["cluster"] = cluster_number + + return stats + + +def calcAllClusterStatistics(distance_matrix, cluster_ids): + """ + Calculates descriptive statistics for all clusters. + Returns a list of dictionaries, each corresponding to a cluster with the statistic quantities + as keys. + + + :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg cluster_ids: a one-dimensional array of the cluster IDs per element + :type cluster_ids: :class:`numpy.ndarray` + + :returns: a list of dictionaries, one dictionary for each cluster with the statistic + quantities as keys. + :rtype: list of dict + + Example usage: + >>> cluster_ids, _= clusterHierarchical(distance_matrix, 4) + >>> all_stats = calcAllClusterStatistics(distance_matrix, cluster_ids) + """ + + cluster_ids = np.asarray(cluster_ids) + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + if len(cluster_ids) != distance_matrix.shape[0]: + raise ValueError(f"Dimension mismatch: cluster_ids length ({len(cluster_ids)}) " + f"does not match distance_matrix frames ({distance_matrix.shape[0]}).") + + clusters = np.unique(cluster_ids) + all_stats = [] + + for c in clusters: + if c <= 0: + continue + + cluster_stats = calcClusterStatistics(distance_matrix, cluster_ids = cluster_ids, cluster_number = int(c)) + all_stats.append(cluster_stats) + + return all_stats + + +def _resolveClusterIndices(cluster_indices, cluster_ids, cluster_number): + """ + Checks whether the input is cluster_indices or cluster_ids and cluster_number, and for the latter returns + the cluster indices. + """ + + has_indices = cluster_indices is not None + has_ids_info = (cluster_ids is not None) and (cluster_number is not None) + + if has_indices == has_ids_info: + raise ValueError("Provide exactly one of 'cluster_indices' or both 'cluster_ids' and 'cluster_number'.") + + if has_indices: + cluster_indices = np.asarray(cluster_indices) + + if cluster_indices.ndim != 1: + raise ValueError(f"cluster_indices must be a 1D array, but shape {cluster_indices.shape} was given.") + if cluster_indices.size == 0: + raise ValueError("cluster_indices must not be empty.") + else: + cluster_ids = np.asarray(cluster_ids) + cluster_indices = getCluster(cluster_ids, cluster_number) + + return cluster_indices + + +def _validateDistanceMatrix(distance_matrix): + """ + Validates that a distance matrix is a non-empty square 2D NumPy array. + """ + + distance_matrix = np.asarray(distance_matrix) + if distance_matrix.ndim != 2: + raise ValueError(f"The distance matrix must be a 2D array, but got shape {distance_matrix.shape}.") + if distance_matrix.shape[0] != distance_matrix.shape[1]: + raise ValueError(f"distance_matrix must be square, but got shape {distance_matrix.shape}.") + if distance_matrix.size == 0: + raise ValueError("distance_matrix cannot be empty.") + + return distance_matrix + + +def showClusterStatisticsTable(all_stats, dissimilarity="RMSD", units="Å", show=True, **kwargs): + """ + Prints the cluster statistics in a table, where each column corresponds to a cluster. + Expects either the 'stats' dictionary generated by :func:`calcClusterStatistics` or + the 'all_stats' list of dictionaries generated by :func:`calcAllClusterStatistics`. + + + :arg all_stats: a dictionary, or list of dictionaries, with the cluster descriptive statistics + Recommended to generate the dictionaries from the functions :func:`calcClusterStatistics` + of a single cluster or :func:`calcAllClusterStatistics`. + :type all_stats: dict, or list of dict + + :arg dissimilarity: the dissimilarity measure + Default is `"RMSD"` + :type dissimilarity: str + + :arg units: the units of the dissimilarity measure + Default is `"Å"` + :type units: str + + :arg show: whether to print the table + Default is `True` + :type show: bool + + :arg **kwargs: keyword arguments passed directly to the ``tabulate`` function + :type **kwargs: dict + + :returns: a dictionary containing the table, the row labels and the headers, with their names as keys + :rtype: dict + + Example usage: + >>> cluster_ids, _ = clusterHierarchical(distance_matrix) + >>> all_stats = calcAllClusterStatistics(distance_matrix, cluster_ids) + >>> showClusterStatisticsTable(all_stats) + """ + + try: + from tabulate import tabulate + except ImportError: + raise ImportError("The 'tabulate' package is required to display the table. " + "Please install it using 'pip install tabulate'.") + + + if isinstance(all_stats, dict): + all_stats = [all_stats] + elif not isinstance(all_stats, list): + raise TypeError(f"all_stats must be a dict or list of dicts, but got {type(all_stats).__name__}.") + + if not all_stats: + raise ValueError("all_stats cannot be empty.") + + required = {"population", "pct", "medoid_global", "medoid_local", "mean", "std", + "median", "iqr", "p95", "max"} + + headers = [] + for i, stat in enumerate(all_stats): + if not isinstance(stat, dict): + raise TypeError(f"Expected a dictionary in all_stats, but got {type(stat).__name__} at index {i}.") + + cluster_id = stat.get('cluster', i + 1) + headers.append(f"Cluster {int(cluster_id)}") + + missing = required - stat.keys() + if missing: + raise ValueError(f"Statistics dictionary is missing required keys: {sorted(missing)}.") + + metrics = [ + ("population", "Total Frames"), + ("pct", "Population Percentage (%)"), + ("medoid_global", "Medoid Frame (Global)"), + ("medoid_local", "Medoid Frame (Within Cluster)"), + ("mean", f"Mean {dissimilarity} [{units}]"), + ("std", f"Std {dissimilarity} [{units}]"), + ("median", f"Median {dissimilarity} [{units}]"), + ("iqr", f"IQR [{units}]"), + ("p95", f"95th Percentile [{units}]"), + ("max", f"Max {dissimilarity} [{units}]") + ] + + table = [] + row_labels = [] + + for key, label in metrics: + row_labels.append(label) + table.append([stat[key] for stat in all_stats]) + + if not isinstance(show, bool): + raise TypeError("show must be a bool.") + + if show: + kwargs.setdefault('tablefmt', 'fancy_grid') + kwargs.setdefault('floatfmt', '.4f') + kwargs.setdefault('stralign', 'center') + + print(tabulate(table, headers=headers, showindex=row_labels, **kwargs)) + + return {"table": table, "row_labels": row_labels, "headers": headers} + + +def writeClusters(atoms, trajectory, distance_matrix, cluster_ids, write_dcd=True, + align ="protein and backbone", system="system", tag ="cluster"): + """ + Aligns the trajectory once, then exports representative medoid structures as PDB files and + cluster-specific DCD trajectories. + Returns the name of all exported .pdb and .dcd files. + + Note: This function loads all aligned coordinates into memory. It is highly optimized + for speed, provided the trajectory fits within available system RAM. + + :arg atoms: reference structure used for the alignment. + :type atoms: :class:`prody.Atomic` + + :arg trajectory: trajectory containing the coordinate sets to align + :type trajectory: :class:`prody.Trajectory` + + :arg distance_matrix: square, symmetric matrix of pairwise distances between all objects + :type distance_matrix: :class:`numpy.ndarray` + + :arg cluster_ids: a one-dimensional array of the cluster IDs per element + :type cluster_ids: :class:`numpy.ndarray` + + :arg write_dcd: determines whether to save a .dcd file of the cluster frames + Default is ``True`` + :type write_dcd: bool + + :arg align: atom selection used to calculate the alignment transformation. + Must be a valid ProDy selection string. + Default is ``"protein and backbone"`` + :type align: str + + :arg system: the name of the system under investigation + Default is ``"system"`` + :type system: str + + :arg tag: the name of the clustering method + Default is ``"cluster"`` + :type tag: str + + :returns: list of exported filenames + :rtype: list[str] + + Example usage: + >>> pdb = parsePDB("structure.pdb") + >>> dcd = Trajectory("trajectory.dcd") + >>> _, aligned_coords = alignTrajectory(pdb, dcd, select='resname IOA') + >>> distance_matrix = calcPairwiseRMSD(aligned_coords) + >>> cluster_ids, _ = clusterHierarchical(distance_matrix) + >>> writeClusters(pdb, dcd, distance_matrix, cluster_ids, + system="type1_RUN23", tag="hier") + """ + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + cluster_ids = np.asarray(cluster_ids) + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") + if cluster_ids.size == 0: + raise ValueError("cluster_ids is empty.") + if len(cluster_ids) != trajectory.numFrames(): + raise ValueError("cluster_ids must have one entry per trajectory frame.") + + if not isinstance(write_dcd, bool): + raise TypeError(f"write_dcd must be a bool, but got {type(write_dcd).__name__}") + + clusters = np.unique(cluster_ids) + clusters = clusters[clusters > 0] + num_clusters = clusters.size + if num_clusters == 0: + raise ValueError("No clusters were found. All frames are labeled as noise") + + exported_files = [] + + # NOTE: Loads the entire aligned trajectory into memory + _, aligned_coords = alignTrajectory(atoms, trajectory, align=align, select="all") + + for cluster in clusters: + cluster_indices = getCluster(cluster_ids, cluster) + medoid = getClusterMedoid(distance_matrix, cluster_indices = cluster_indices) + + cluster_coords = aligned_coords[cluster_indices] + + cluster_atoms = atoms.copy() + cluster_atoms.setCoords(cluster_coords[0]) + + if len(cluster_coords) > 1: + cluster_atoms.addCoordset(cluster_coords[1:]) + + if write_dcd: + dcd_filename = f"{system}_{tag}_n{num_clusters}_cluster{cluster}.dcd" + writeDCD(dcd_filename, cluster_atoms) + exported_files.append(dcd_filename) + + medoid_atoms = atoms.copy() + medoid_atoms.setCoords(aligned_coords[medoid["global"]]) + + pdb_filename = f"{system}_{tag}_n{num_clusters}_cluster{cluster}_medoid.pdb" + writePDB(pdb_filename, medoid_atoms) + exported_files.append(pdb_filename) + + return exported_files + + +def clusterHierarchical(distance_matrix, method='average', cutoff=None): + """ + Performs bottom-up hierarchical clustering from a pairwise distance matrix. + + Note that the resulting cluster IDs are 1-based (starting at 1, not 0). + + If cutoff is ``None`` or ``'auto'``, it is automatically chosen as the midpoint + of the largest gap between consecutive linkage distances. + + Recommendation: Plot the dendrogram using showDendrogram() and choose the cutoff + manually whenever possible. + + + :arg distance_matrix: either a one-dimensional condensed distance matrix or + a two-dimensional pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg method: linkage criterion used when constructing the linkage matrix. + Default is ``'average'`` + :type method: str + + :arg cutoff: the cutoff distance that determines the number of clusters. + If ``None`` or ``'auto'``, the midpoint of the largest linkage gap + is calculated automatically. + Default is ``None`` + :type cutoff: float, str, None + + :returns: + * a one-dimensional array containing the cluster ID for each object. + IDs are 1-indexed (1 to number of clusters for the cutoff) + * the hierarchical clustering linkage matrix + :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) + + Example usage: + >>> distance_matrix = calcPairwiseRMSD(aligned_coords) + >>> cluster_ids, linkage_matrix = clusterHierarchical(distance_matrix, cutoff='auto') + """ + + from scipy.cluster.hierarchy import linkage, fcluster + + + condensed_distance_matrix = _condenseDistanceMatrix(distance_matrix) + linkage_matrix = linkage(condensed_distance_matrix, method=method) + + if cutoff is None or cutoff == "auto": + cutoff = _calcAutoCutoff(linkage_matrix) + elif isinstance(cutoff, str): + raise ValueError(f"Invalid string for cutoff: '{cutoff}'. Use 'auto', None, or a numeric value.") + + cluster_ids = fcluster(linkage_matrix, t=cutoff, criterion='distance') + + return cluster_ids, linkage_matrix + + +def showDendrogram(distance_matrix=None, linkage_matrix=None, *args, **kwargs): + """ + Plots a hierarchical clustering dendrogram using either a pairwise distance matrix + or a pre-computed linkage matrix. + + By default, the complete dendrogram is shown. To display a truncated dendrogram, + pass ``truncate_mode`` and ``p`` via kwargs. + + Exactly one of ``distance_matrix`` or ``linkage_matrix`` must be provided. + + + :arg distance_matrix: one-dimensional condensed or two-dimensional square pairwise + distance matrix. + Default is ``None`` + :type distance_matrix: :class:`numpy.ndarray` + + :arg linkage_matrix: pre-computed hierarchical clustering linkage matrix. + Default is ``None`` + :type linkage_matrix: :class:`numpy.ndarray` + + :arg *args: positional arguments passed directly to SciPy's ``dendrogram`` function. + :type *args: tuple + + :arg method: linkage criterion used when constructing the linkage matrix from + ``distance_matrix``. (Passed via kwargs). + Default is ``"average"`` + :type method: str + + :arg cutoff: dendrogram color threshold. If ``"auto"``, the cutoff is placed at the midpoint + of the largest linkage gap. (Passed via kwargs). + Default is ``None`` + :type cutoff: float, str + + :arg title: title of the generated plot. (Passed via kwargs). + Default is ``"Hierarchical Clustering Dendrogram"`` + :type title: str + + :arg xlabel: label for the x-axis. (Passed via kwargs). + Default is ``"Index"`` + :type xlabel: str + + :arg ylabel: label for the y-axis. (Passed via kwargs). + Default is ``"Distance"`` + :type ylabel: str + + :arg ax: axes on which to draw the plot. (Passed via kwargs). + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :arg **kwargs: additional keyword arguments passed to SciPy's ``dendrogram`` + (e.g., ``truncate_mode``, ``p``, ``color_threshold``, ``no_labels``). + :type **kwargs: dict + + :returns: the matplotlib axes and the linkage matrix. + :rtype: tuple(:class:`matplotlib.axes.Axes`, :class:`numpy.ndarray`) + + Example usage: + >>> import matplotlib.pyplot as plt + >>> distance_matrix = calcPairwiseRMSD(aligned_coords) + >>> plt.figure(figsize=(12, 8)) + >>> ax, linkage = showDendrogram(distance_matrix, cutoff='auto', truncate_mode='lastp', p=30) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + from scipy.cluster.hierarchy import dendrogram, linkage + + + if (distance_matrix is None) == (linkage_matrix is None): + raise ValueError("Provide exactly one of 'distance_matrix' or 'linkage_matrix'.") + + method = kwargs.pop('method', 'average') + cutoff = kwargs.pop('cutoff', None) + + if linkage_matrix is None: + condensed_distance_matrix = _condenseDistanceMatrix(distance_matrix) + linkage_matrix = linkage(condensed_distance_matrix, method=method) + else: + linkage_matrix = np.asarray(linkage_matrix) + if linkage_matrix.ndim != 2 or linkage_matrix.shape[1] != 4 or linkage_matrix.shape[0] == 0: + raise ValueError("linkage_matrix must be a non-empty array with shape (n_samples-1, 4).") + + if cutoff == "auto": + cutoff = _calcAutoCutoff(linkage_matrix) + elif isinstance(cutoff, str): + raise ValueError(f"Invalid string for cutoff: '{cutoff}'. Use 'auto' or a numeric value.") + + ax = kwargs.pop('ax', None) + if ax is None: + ax = plt.gca() + + title = kwargs.pop('title', "Hierarchical Clustering Dendrogram") + xlabel = kwargs.pop('xlabel', "Index") + ylabel = kwargs.pop('ylabel', "Distance") + + # SciPy Dendogram Defaults + kwargs.setdefault('truncate_mode', None) + kwargs.setdefault('p', 30) + kwargs.setdefault('no_labels', True) + + if cutoff is not None: + kwargs.setdefault('color_threshold', cutoff) + + with plt.rc_context({"lines.linewidth": 0.6}): + dendrogram(linkage_matrix, *args, ax=ax, **kwargs) + + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.set_title(title) + + if cutoff is not None: + ax.axhline(cutoff, color="black", linestyle="--", linewidth=1.2, label=f"Cutoff ({cutoff:.2f})") + ax.legend() + + return ax, linkage_matrix + + +def _calcAutoCutoff(linkage_matrix): + """ + Calculates the optimal distance cutoff based on the midpoint of the largest gap + in the linkage matrix. + """ + + linkage_matrix = np.asarray(linkage_matrix) + merge_dist = linkage_matrix[:, 2] # locations of the merges + + if len(merge_dist) == 0: + raise ValueError("At least two frames are required for clustering") + elif len(merge_dist) == 1: + return merge_dist[0] + 1.0 + + gaps = np.diff(merge_dist) + largest_gap_idx = np.argmax(gaps) + + return (merge_dist[largest_gap_idx] + merge_dist[largest_gap_idx + 1]) / 2.0 + + +def _condenseDistanceMatrix(distance_matrix): + """ + Converts a pairwise distance matrix into condensed form. + + Accepts either: + * a condensed one-dimensional distance matrix, or + * a square two-dimensional distance matrix. + + Returns the condensed distance matrix suitable for scipy.cluster.hierarchy.linkage. + """ + from scipy.spatial.distance import squareform + + + distance_matrix = np.asarray(distance_matrix) + + if distance_matrix.size == 0: + raise ValueError("distance_matrix cannot be empty.") + + if distance_matrix.ndim == 2: + if distance_matrix.shape[0] != distance_matrix.shape[1]: + raise ValueError(f"distance_matrix must be square, but got shape {distance_matrix.shape}.") + condensed_distance_matrix = squareform(distance_matrix) + elif distance_matrix.ndim == 1: + condensed_distance_matrix = distance_matrix + else: + raise ValueError("distance_matrix must be either a condensed 1D array or a square 2D array.") + + return condensed_distance_matrix + + +def clusterKMedoids(distance_matrix, k, method='alternate', method_sklearn='pam', + initial_medoids=None, seed=None, max_iter=100, n_init=10): + """ + Performs K-Medoids clustering using various algorithms. + + This function acts as a facade, routing the clustering task to the specified backend + ('alternate', 'pam', or 'sklearn'). 'alternate' is generally faster, while 'pam' typically + finds solutions with lower cost. 'sklearn' needs the sklearn_extra.cluster module to be installed. + + Note that the resulting cluster IDs are 1-based (starting at 1, not 0). + + + :arg distance_matrix: square, symmetric matrix of pairwise distances between all objects + :type distance_matrix: :class:`numpy.ndarray` + + :arg k: prespecified number of clusters to form + :type k: int + + :arg method: the clustering algorithm to use. + Options are 'alternate' (custom Alternating Medoids), 'pam' (custom PAM), 'sklearn' (sklearn_extra) + Default is 'alternate' + :type method: str + + :arg method_sklearn: the specific method to pass to sklearn_extra if method = 'sklearn' is chosen. + Options are 'pam', 'alternate' + Default is 'pam' + :type method_sklearn: str + + :arg initial_medoids: one-dimensional array of indices to use as starting medoids. + Default is ``None`` and the starting medoids are picked randomly. + It is not supported for the 'sklearn' method. + :type initial_medoids: :class:`numpy.ndarray` + + :arg seed: random seed for reproducibility + Default is ``None`` + :type seed: int + + :arg max_iter: maximum number of iterations per a single run + Default is ``100`` + :type max_iter: int + + :arg n_init: number of times the algorithm will be run with different initial medoids + and the best result with the lowest cost is returned + + :returns: a tuple of: + * an one-dimensional array containing the cluster ID for each object. + IDs are 1-indexed (1 to k) + * a one-dimensional array of shape (k,) containing the indices of the final cluster + medoids + * the final sum of distances from each point to its nearest medoid + :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`, float) + + Example usage: + >>> distance_matrix = calcPairwiseRMSD(aligned_coords) + >>> cluster_ids, medoids, _ = clusterKMedoids(distance_matrix, 4, method = 'alternate', + seed = 42, n_init = 30) + """ + + distance_matrix = _validateDistanceMatrix(distance_matrix) + n = distance_matrix.shape[0] + + if not isinstance(k, (int, np.integer)): + raise TypeError(f"k must be an integer, got {type(k).__name__}") + if k <= 0 or k > n: + raise ValueError(f"k must be between 1 and {n} (the number of points).") + + if not isinstance(max_iter, (int, np.integer)) or max_iter <= 0: + raise ValueError("max_iter must be a positive integer.") + + if not isinstance(n_init, (int, np.integer)) or n_init <= 0: + raise ValueError("n_init must be a positive integer.") + + if seed is not None and not isinstance(seed, (int, np.integer)): + raise TypeError("seed must be an integer or None.") + + if initial_medoids is not None: + if method == 'sklearn': + raise ValueError("The 'sklearn' method does not support custom 'initial_medoids'. Use 'alternate' or 'pam' instead.") + + initial_medoids = np.asarray(initial_medoids, dtype=int) + if initial_medoids.ndim != 1: + raise ValueError("initial_medoids must be a 1D array.") + if len(initial_medoids) != k: + raise ValueError(f"Number of initial medoids ({len(initial_medoids)}) must equal k ({k}).") + if len(np.unique(initial_medoids)) != k: + raise ValueError("initial_medoids must be unique.") + if np.any((initial_medoids < 0) | (initial_medoids >= n)): + raise ValueError(f"initial_medoids contain indices out of bounds (must be 0 to {n-1}).") + + n_init = 1 + + if method == 'alternate': + return _clusterKMedoidsAlternating(distance_matrix, k, initial_medoids, seed, max_iter, n_init) + elif method == 'pam': + return _clusterKMedoidsPAM(distance_matrix, k, initial_medoids, seed, max_iter, n_init) + elif method == 'sklearn': + if method_sklearn not in ['alternate', 'pam']: + raise ValueError(f"method_sklearn '{method_sklearn}' is not valid. Options: 'alternate', 'pam'") + return _clusterKMedoidsSklearn(distance_matrix, k, method_sklearn, seed, max_iter, n_init) + else: + raise ValueError(f"Method '{method}' is not valid. Options: 'alternate', 'pam', 'sklearn'") + + +def _cost(distance_matrix, medoids): + """Calculates the sum of distances from all points to their nearest medoid.""" + distances = distance_matrix[:, medoids] + closest_distance = np.min(distances, axis=1) + return np.sum(closest_distance) + + +def _clusterKMedoidsAlternating(distance_matrix, k, initial_medoids, seed, max_iter, n_init): + """K-Medoids clustering with an Alternating Medoids algorithm""" + + n = distance_matrix.shape[0] + rng = np.random.default_rng(seed) + + best_cost = np.inf + best_medoids = None + + for run in range(n_init): + if initial_medoids is None: + medoids = np.asarray(rng.choice(n, size = k, replace = False)) + else: + medoids = initial_medoids.copy() + + for iteration in range(max_iter): + distances_to_medoids = distance_matrix[:, medoids] + cluster_ids = np.argmin(distances_to_medoids, axis = 1) + + new_medoids = np.zeros_like(medoids) + + for i in range(len(medoids)): + cluster_members = np.where(cluster_ids == i)[0] + + # Handles empty clusters by retaining the old medoid + if len(cluster_members) == 0: + new_medoids[i] = medoids[i] + continue + + intra_cluster_distances = distance_matrix[np.ix_(cluster_members, cluster_members)] + sum_distances = intra_cluster_distances.sum(axis=1) + best_medoid_idx_in_cluster = np.argmin(sum_distances) + new_medoids[i] = cluster_members[best_medoid_idx_in_cluster] + + current_cost = _cost(distance_matrix, new_medoids) + + if np.array_equal(medoids, new_medoids): + break + + medoids = new_medoids + + if current_cost < best_cost: + best_cost = current_cost + best_medoids = medoids.copy() + + distances = distance_matrix[:, best_medoids] + cluster_ids = np.argmin(distances, axis=1) + 1 + + return cluster_ids, best_medoids, best_cost + + +def _clusterKMedoidsPAM(distance_matrix, k, initial_medoids, seed, max_iter, n_init): + """K-Medoids clustering with PAM (Partitioning Around Medoids)""" + + n = distance_matrix.shape[0] + rng = np.random.default_rng(seed) + + best_cost = np.inf + best_medoids = None + + for run in range(n_init): + if initial_medoids is None: + medoids = np.asarray(rng.choice(n, size=k, replace=False)) + else: + medoids = initial_medoids.copy() + + # Using a set for removal/addition operations + non_medoids = set(i for i in range(n) if i not in medoids) + current_cost = _cost(distance_matrix, medoids) + + for iteration in range(max_iter): + best_cost_swap = current_cost + best_swap = None + + for medoid_idx, old_medoid in enumerate(medoids): + for new_medoid in non_medoids: + candidate = medoids.copy() + candidate[medoid_idx] = new_medoid + candidate_cost = _cost(distance_matrix, candidate) + + if candidate_cost < best_cost_swap: + best_cost_swap = candidate_cost + best_swap = (medoid_idx, old_medoid, new_medoid) + + if best_swap is None: + break + + idx, old_medoid, new_medoid = best_swap + medoids[idx] = new_medoid + non_medoids.remove(new_medoid) + non_medoids.add(old_medoid) + current_cost = best_cost_swap + + if current_cost < best_cost: + best_cost = current_cost + best_medoids = medoids.copy() + + distances = distance_matrix[:, best_medoids] + cluster_ids = np.argmin(distances, axis=1) + 1 + + return cluster_ids, best_medoids, best_cost + + +def _clusterKMedoidsSklearn(distance_matrix, k, method_sklearn, seed, max_iter, n_init): + """K-Medoids clustering with sklearn_extra""" + try: + from sklearn_extra.cluster import KMedoids + except ImportError: + raise ImportError("The 'sklearn_extra' package is required for this K-Medoids approach. " + "Please install it using 'pip install scikit-learn-extra'.") + + + best_cost = np.inf + best_labels = None + best_medoids = None + + rng = np.random.default_rng(seed) + + for run in range(n_init): + kmedoids = KMedoids(n_clusters=k, metric='precomputed', method=method_sklearn, init='random', + max_iter=max_iter, random_state=int(rng.integers(0, 1000000))) + + kmedoids.fit(distance_matrix) + medoids = kmedoids.medoid_indices_ + + cost = np.sum(np.min(distance_matrix[:, medoids], axis=1)) + + if cost < best_cost: + best_cost = cost + best_labels = kmedoids.labels_ + 1 + best_medoids = medoids.copy() + + return best_labels, best_medoids, best_cost + + +def clusterDBSCAN(distance_matrix, eps=None, minPts=None, method='custom'): + """ + Performs DBSCAN clustering using various algorithms. + + This function acts as a facade, routing the clustering task to the specified backend + ('custom' or 'sklearn'). 'custom' uses a built-in implementation without additional + dependencies beyond NumPy. 'sklearn' needs the sklearn.cluster module to be installed. + + Note that the resulting cluster IDs are 1-based (starting at 1, not 0). + Noise points are labeled -1. + + + :arg distance_matrix: square, symmetric matrix of pairwise distances between all objects + :type distance_matrix: :class:`numpy.ndarray` + + :arg eps: the "radius" of the neighborhood within which we count neighbors + Default is ``None`` and automatically the median of the distances + in the distance matrix is used. + Ideally use :func:`showReachabilityPlot` to determine + manually the most suitable eps. + :type eps: float + + :arg minPts: the minimum number of neighbors required for a point to be considered + a core point + Default is ``None`` and automatically the 5% of the total objects, + or for less than 20 2 is used. + Ideally choose manually the most suitable minPts + :type minPts: int + + :arg method: the clustering algorithm to use + Options are 'custom' and 'sklearn' + Default is 'custom' because it needs no module installation + :type method: str + + :returns: a tuple of: + * a one-dimensional array containing the cluster ID for each object. + IDs are 1-indexed and noise corresponds to -1. + * a one-dimensional array of the frame indices corresponding to noise + :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) + + Example usage: + >>> distance_matrix = calcPairwiseRMSD(aligned_coords) + >>> cluster_ids, _ = clusterDBSCAN(distance_matrix, eps = 1.8, minPts = 30) + """ + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + if eps is None: + eps = _calcAutoEps(distance_matrix) + _validateEps(eps, distance_matrix) + + if minPts is None: + minPts = _calcAutoMinPts(distance_matrix) + _validateMinPts(minPts, distance_matrix) + + if method == 'custom': + return _clusterDBSCANCustom(distance_matrix, eps, minPts) + elif method == 'sklearn': + return _clusterDBSCANSklearn(distance_matrix, eps, minPts) + else: + raise ValueError(f"method can be either 'custom' or 'sklearn', but got {method}.") + + +def _clusterDBSCANCustom(distance_matrix, eps, minPts): + """DBSCAN clustering with custom algorithm""" + total_points = distance_matrix.shape[0] + + labels = np.zeros(total_points, dtype=int) + cluster_id = 0 + + for p in range(total_points): + + if labels[p] != 0: + continue + + neighbors, = np.where(distance_matrix[p] <= eps) + + if len(neighbors) < minPts: + labels[p] = -1 + else: + cluster_id += 1 + labels[p] = cluster_id + + candidate_set = [n for n in neighbors if n != p] + while candidate_set: + q = candidate_set.pop() + + if labels[q] == -1: + labels[q] = cluster_id + + if labels[q] != 0: + continue + + labels[q] = cluster_id + + q_neighbors, = np.where(distance_matrix[q] <= eps) + if len(q_neighbors) >= minPts: + for n in q_neighbors: + if labels[n] == 0: + candidate_set.append(n) + elif labels[n] == -1: + labels[n] = cluster_id + + noise_frames, = np.where(labels == -1) + + return labels, noise_frames + + +def _clusterDBSCANSklearn(distance_matrix, eps, minPts): + """DBSCAN clustering with sklearn""" + + try: + from sklearn.cluster import DBSCAN + except ImportError: + raise ImportError("The 'sklearn' package is required for this DBSCAN approach. " + "Please install it using 'pip install scikit-learn'.") + + + dbscan = DBSCAN(eps = eps, min_samples = minPts, metric = "precomputed") + labels = dbscan.fit_predict(distance_matrix) + + cluster_ids = np.copy(labels) + cluster_ids[cluster_ids >= 0] += 1 + + noise_frames, = np.where(labels == -1) + + return cluster_ids, noise_frames + + +def _calcAutoEps(distance_matrix): + """Automatically determine the DBSCAN eps parameter""" + + import warnings + + + pairwise_distances = distance_matrix[np.triu_indices_from(distance_matrix, k = 1)] + eps = float(np.median(pairwise_distances)) + warnings.warn(f"No eps provided. Automatically chosen at {eps:.3f}. Ideally provide your own value.") + return eps + + +def _validateEps(eps, distance_matrix): + """Validate the DBSCAN eps parameter""" + + import warnings + + + maxDistance = np.max(distance_matrix) + + if not isinstance(eps, (float, int, np.floating, np.integer)): + raise TypeError(f"eps must be a numeric value, but got {type(eps).__name__}") + + if eps <= 0: + raise ValueError("eps must be positive") + elif eps > maxDistance: + warnings.warn(f"eps ({eps}) is greater than the maximum pairwise distance ({maxDistance:.3f}).\n" + "All frames will be clustered together with no noise.") + + +def _calcAutoMinPts(distance_matrix): + """ Automatically determine the minPts parameter""" + + import warnings + + + total_points = distance_matrix.shape[0] + minPts = max(2, int(total_points // 20)) + warnings.warn(f"No minPts provided. Automatically chosen at {minPts}. Ideally provide your own value.") + return minPts + + +def _validateMinPts(minPts, distance_matrix): + """Validate minPts parameter""" + total_points = distance_matrix.shape[0] + + if not isinstance(minPts, (int, np.integer)): + raise TypeError(f"minPts must be a positive integer, but got {type(minPts).__name__}") + + if minPts <= 0 or minPts > total_points: + raise ValueError(f"minPts must be between 1 and {total_points}.") + + +def showReachabilityPlot(distance_matrix, *args, minPts=None, method='custom', eps=None, fill=True, **kwargs): + """ + Plots the reachability plot using a simplified OPTICS algorithm. + The reachability plot should be used to determine the most suitable eps + parameter for DBSCAN. + + This function acts as a facade, routing the ordering task to the specified backend + ('custom' or 'sklearn'). 'custom' uses a built-in implementation without additional + dependencies beyond NumPy. 'sklearn' needs the sklearn.cluster module to be installed. + + + :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. + :type distance_matrix: :class:`numpy.ndarray` + + :arg *args: positional arguments passed directly to Matplotlib's ``plot`` function. + :type *args: tuple + + :arg minPts: the minimum number of neighbors required for a point to be considered + a core point. + Default is ``None`` and automatically the 5% of the total objects, + or for less than 20, 2 is used. + :type minPts: int + + :arg method: the OPTICS algorithm to use ('custom' or 'sklearn'). + Default is 'custom'. + :type method: str + + :arg eps: the "radius" of the neighborhood within which we count neighbors. + If `'auto'`, the median of the pairwise distances is used. + Default is ``None``. + :type eps: float, str + + :arg fill: whether to fill the area beneath the curve and the valleys below eps. + Default is ``True``. + :type fill: bool + + :arg **kwargs: additional keyword arguments passed to Matplotlib's ``plot`` function + :type **kwargs: dict + + :returns: the Matplotlib axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> plt.figure() + >>> showReachabilityPlot(distance_matrix, minPts=20, eps=1.8) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + + + distance_matrix = _validateDistanceMatrix(distance_matrix) + + if minPts is None: + minPts = _calcAutoMinPts(distance_matrix) + _validateMinPts(minPts, distance_matrix) + + if not isinstance(fill, bool): + raise TypeError(f"fill must be a bool, but got {type(fill).__name__}") + + if method == 'custom': + reachability, ordering = _orderOPTICSCustom(distance_matrix, minPts) + elif method == 'sklearn': + reachability, ordering = _orderOPTICSSklearn(distance_matrix, minPts) + else: + raise ValueError(f"method must be either 'custom' or 'sklearn', but got {method}") + + ax = kwargs.pop('ax', None) + if ax is None: + ax = plt.gca() + + title = kwargs.pop('title', "OPTICS Reachability Plot") + xlabel = kwargs.pop('xlabel', "Frames (Sorted by OPTICS)") + ylabel = kwargs.pop('ylabel', "Reachability Distance") + + ax.set_title(title) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + + y = reachability[ordering] + x = np.arange(len(y)) + + kwargs.setdefault('color', "#36454F") + kwargs.setdefault('lw', 1.5) + + ax.plot(x, y, *args, **kwargs) + + if fill: + ax.fill_between(x, 0, y, color='black', alpha=0.4) + + if eps == 'auto': + eps = _calcAutoEps(distance_matrix) + + if eps is not None: + _validateEps(eps, distance_matrix) + + line_width = kwargs.get('lw', kwargs.get('linewidth', 1.5)) + ax.axhline(y=eps, color='black', linestyle='--', linewidth=line_width) + + if fill: + below_eps_mask = (y <= eps) + segments = [] + start = None + + for i, below in enumerate(below_eps_mask): + if below and start is None: + start = i + elif not below and start is not None: + segments.append((start, i)) + start = None + + if start is not None: + segments.append((start, len(below_eps_mask))) + + colors = plt.cm.tab10(np.linspace(0, 1, len(segments))) + + for (start, end), valley_color in zip(segments, colors): + ax.fill_between(x[start:end], y[start:end], eps, color=valley_color) + + if 'label' in kwargs: + ax.legend() + + ax.grid(axis='y', linestyle='--', linewidth=0.8, alpha=0.3) + + return ax + + +def _orderOPTICSCustom(distance_matrix, minPts): + """OPTICS algorithm with built-in modules and NumPy""" + + import heapq + + + total_points = distance_matrix.shape[0] + + sorted_distances = np.sort(distance_matrix, axis=1) + core_distances = sorted_distances[:, minPts - 1] + + reachability = np.full(total_points, np.inf) + processed = np.zeros(total_points, dtype=bool) + ordering = [] + + def _updateSeeds(idx): + new_reaches = np.maximum(core_distances[idx], distance_matrix[idx, :]) # Reachability distance definition + update_mask = (~processed) & (new_reaches < reachability) + points_to_update, = np.where(update_mask) + reachability[update_mask] = new_reaches[update_mask] + + for j in points_to_update: + heapq.heappush(seeds, (reachability[j], j)) + + for i in range(total_points): + if processed[i]: + continue + + processed[i] = True + ordering.append(i) + seeds = [] + + _updateSeeds(i) + + while seeds: + current_reach, q = heapq.heappop(seeds) + + if processed[q]: + continue + + processed[q] = True + ordering.append(q) + + _updateSeeds(q) + + return reachability, np.array(ordering) + + +def _orderOPTICSSklearn(distance_matrix, minPts): + """OPTICS algorithm with sklearn""" + + try: + from sklearn.cluster import OPTICS + except ImportError: + raise ImportError("The 'sklearn' package is required for this OPTICS approach. " + "Please install it using 'pip install scikit-learn'.") + + optics = OPTICS(min_samples = minPts, metric = 'precomputed') + optics.fit(distance_matrix) + reachability = optics.reachability_ + ordering = optics.ordering_ + + return reachability, ordering \ No newline at end of file diff --git a/prody/proteins/interactions.py b/prody/proteins/interactions.py index 1c7b4ebcc..a4208f62f 100644 --- a/prody/proteins/interactions.py +++ b/prody/proteins/interactions.py @@ -53,13 +53,7 @@ 'calcSminaBindingAffinity', 'calcSminaPerAtomInteractions', 'calcSminaTermValues', 'showSminaTermValues', 'showPairEnergy', 'checkNonstandardResidues', 'saveInteractionsAsDummyAtoms', 'createFoldseekAlignment', 'runFoldseek', 'runDali', - 'runBLAST', 'extractMultiModelPDB', 'calcSignatureInteractions', - 'alignTrajectory', 'calcRMSDfromReference', 'showRMSDfromReference', - 'calcPairwiseRMSD', 'showPairwiseRMSDHeatmap', 'showRMSDHistogram', - 'calcClusterPopulations', 'getCluster', 'getClusterMedoid', - 'calcClusterStatistics', 'calcAllClusterStatistics', 'showClusterStatisticsTable', - 'showClusterRMSDComparison', 'clusterHierarchical', 'showDendrogram', - 'clusterKMedoids', 'clusterDBSCAN', 'showReachabilityPlot', 'writeClusters'] + 'runBLAST', 'extractMultiModelPDB', 'calcSignatureInteractions'] def cleanNumbers(listContacts): @@ -4247,1852 +4241,6 @@ def calcSignatureInteractions(PDB_folder, **kwargs): - -def alignTrajectory(atoms, trajectory, align = 'protein and backbone', select = 'all'): - """ - Aligns each trajectory frame to the reference structure and returns a tuple of the - reference coordinates and the aligned coordinates of the selected atoms. - - The trajectory frames are aligned to the reference structure using the atoms specified by - ``align``. After alignment, the coordinates of the atoms specified by ``select`` are - extracted for every frame. - - - :arg atoms: reference structure used for the alignment. - :type atoms: :class:`prody.Atomic` - - :arg trajectory: trajectory containing the coordinate sets to align - :type trajectory: :class:`prody.Trajectory` - - :arg align: atom selection used to calculate the alignment transformation. - Must be a valid ProDy selection string. - Default is ``"protein and backbone"`` - :type align: str - - :arg select: atom selection whose coordinates are returned. - Must be a valid ProDy selection string. - Default is ``"all"`` - :type select: str - - :returns: a tuple containing: - * ref_coords (numpy.ndarray): coordinates of the selected atoms in the reference structure. - * aligned_coords (numpy.ndarray): aligned coordinates of the selected atoms for every - trajectory frame with shape ``(n_frames, n_atoms, 3)``. - :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) - - Example usage: - >>> pdb = prody.parsePDB("structure.pdb") - >>> dcd = prody.Trajectory("trajectory.dcd") - >>> ref_coords, aligned_coords = prody.alignTrajectory(pdb, dcd, select='resname IOA') - """ - - if trajectory.numAtoms() != atoms.numAtoms(): - raise ValueError("Trajectory atoms count does not match structure atoms count.") - - # Save original coordinates to restore in the end - orig_coords = atoms.getCoords().copy() - - atom_align = atoms.select(align) - if atom_align is None: - raise ValueError(f"No atoms match '{align}' in the structure.") - ref_align = atom_align.copy() - - atom_select = atoms.select(select) - if atom_select is None: - raise ValueError(f"No atoms match '{select}' in the structure.") - ref_coords = atom_select.getCoords().copy() - - trajectory.link(atoms) # linking trajectory to update coordinates frame-by-frame - - n_frames = trajectory.numFrames() - n_atoms = atom_select.numAtoms() - - trajectory.reset() - aligned_coords = np.zeros((n_frames, n_atoms, 3)) - - try: - for i, frame in enumerate(trajectory): - trans = calcTransformation(atom_align, ref_align) - trans.apply(atom_select) - aligned_coords[i] = atom_select.getCoords() - finally: - atoms.setCoords(orig_coords) - trajectory.reset() - - return ref_coords, aligned_coords - - -def calcRMSDfromReference(reference_coords, aligned_coords): - """ - Calculates the RMSD of each aligned trajectory frame from the reference coordinates. - Uses a vectorized approach for better efficiency. - - - :arg reference_coords: reference coordinates with shape ``(n_atoms, 3)``. - :type reference_coords: :class:`numpy.ndarray` - - :arg aligned_coords: aligned coordinates with shape ``(n_frames, n_atoms, 3)``. - Recommended to generate them using :func:`alignTrajectory`. - :type aligned_coords: :class:`numpy.ndarray` - - :returns: an array containing the RMSD of each frame from the reference - :rtype: :class:`numpy.ndarray` - - Example usage: - >>> pdb = prody.parsePDB("structure.pdb") - >>> dcd = prody.Trajectory("trajectory.dcd") - >>> ref_coords, aligned_coords = prody.alignTrajectory(pdb, dcd, select='resname IOA') - >>> rmsd = prody.calcRMSDfromReference(ref_coords, aligned_coords) - """ - - aligned_coords = np.asarray(aligned_coords) - reference_coords = np.asarray(reference_coords) - - if reference_coords.ndim != 2: - raise ValueError(f"reference_coords must be a 2D array of shape (n_atoms, 3), but got {reference_coords.shape}") - - if aligned_coords.ndim != 3: - raise ValueError(f"aligned_coords must have shape (n_frames, n_atoms, 3), but got {aligned_coords.shape}.") - - if reference_coords.shape != aligned_coords.shape[1:]: - raise ValueError(f"Incompatible shapes: reference is {reference_coords.shape}, but aligned frames have {aligned_coords.shape[1:]}.") - - # Vectorized RMSD calculation - sd = (aligned_coords - reference_coords)**2 - msd = np.mean(sd.sum(axis=2), axis=1) - rmsd = np.sqrt(msd) - - return rmsd - - -def showRMSDfromReference(rmsd_array, title = "Frame-to-Reference RMSD", - color = "#36454F", lw = 0.75, label = None, ax = None): - """ - Plots the RMSD of an aligned trajectory from the reference on the current axis. - - - :arg rmsd_array: one-dimensional array containing the RMSD values (in Å) of each - trajectory frame relative to the reference structure. - This array is typically generated using :func:`calcRMSDfromReference`. - :type rmsd_array: :class:`numpy.ndarray` - - :arg title: the title of the generated plot. - Default is ``"Frame-to-Reference RMSD"`` - :type title: str - - :arg color: the color of the plot. - Default is ``"#36454F"`` - :type color: str - - :arg lw: the width of the line. - Default is ``0.75`` - :type lw: float - - :arg label: the label for the plot. - Default is ``None`` - :type label: str - - :arg ax: axes on which to draw the plot. - Default is ``None`` and the current axes are used. - :type ax: :class:`matplotlib.axes.Axes` - - :returns: the Matplotlib axes containing the plot. - :rtype: :class:`matplotlib.axes.Axes` - - Example usage: - >>> import matplotlib.pyplot as plt - >>> rmsd_run1 = prody.calcRMSDfromReference(ref_coords, aligned_coords1) - >>> rmsd_run2 = prody.calcRMSDfromReference(ref_coords, aligned_coords2) - >>> plt.figure(figsize=(8, 6)) - >>> prody.showRMSDfromReference(rmsd_run1, label="Run 1") - >>> prody.showRMSDfromReference(rmsd_run2, color="red", label="Run 2") - >>> plt.show() - """ - - import matplotlib.pyplot as plt - - - if ax is None: - ax = plt.gca() - - ax.plot(rmsd_array, color = color, lw = lw, label = label) - ax.set_xlabel("# Frame") - ax.set_ylabel("RMSD from Reference [Å]") - ax.set_title(title) - ax.grid(alpha=0.3) - - if label is not None: - ax.legend() - - return ax - - -def calcPairwiseRMSD(aligned_coords): - """ - Calculates the frame-to-frame pairwise RMSD matrix using aligned structures - and returns a symmetric distance matrix. - Uses a vectorized approach for better efficiency. - - - :arg aligned_coords: aligned coordinates with shape ``(n_frames, n_atoms, 3)``. - Recommended to generate them using :func:`alignTrajectory`. - :type aligned_coords: :class:`numpy.ndarray` - - :returns: symmetric matrix containing the pairwise RMSD between all frames with shape - ``(n_frames, n_frames)``. - :rtype: :class:`numpy.ndarray` - - Example usage: - >>> ref_coords, aligned_coords = prody.alignTrajectory(pdb, dcd, select='resname IOA') - >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) - """ - - from scipy.spatial.distance import cdist - - - aligned_coords = np.asarray(aligned_coords) - - if aligned_coords.ndim != 3 or aligned_coords.shape[2] != 3: - raise ValueError(f"aligned_coords must have shape (n_frames, n_atoms, 3), but got {aligned_coords.shape}.") - - if aligned_coords.shape[0] == 0: - raise ValueError("aligned_coords contains no frames.") - - if aligned_coords.shape[1] == 0: - raise ValueError("aligned_coords contains no atoms.") - - - n_frames, n_atoms, _ = aligned_coords.shape - - # Flatten frame coordinates (n_frames, n_atoms, 3) -> (n_frames, n_atoms * 3) - coords_flat = aligned_coords.reshape(n_frames, 3 * n_atoms) - - euclidean_dists = cdist(coords_flat, coords_flat, metric = 'euclidean') - distance_matrix = euclidean_dists / np.sqrt(n_atoms) - - return distance_matrix - - -def showPairwiseRMSDHeatmap(distance_matrix, title = "Pairwise RMSD Distance Matrix", - label = "RMSD [Å]", cmap = "viridis", ax = None): - """ - Plots the heatmap of the frame-to-frame pairwise RMSDs using the distance matrix - on the current axis. - - - :arg distance_matrix: two-dimensional array constituting the pairwise distance matrix. - Typically generated using :func:`calcPairwiseRMSD`. - :type distance_matrix: :class:`numpy.ndarray` - - :arg title: the title of the generated plot. - Default is ``"Pairwise RMSD Distance Matrix"`` - :type title: str - - :arg label: the label for the heatmap legend. - Default is ``"RMSD [Å]"`` - :type label: str - - :arg cmap: the colormap for the heatmap. - Default is ``"viridis"`` - :type cmap: str - - :arg ax: axes on which to draw the plot. - Default is ``None`` and the current axes are used. - :type ax: :class:`matplotlib.axes.Axes` - - :returns: the Matplotlib axes containing the heatmap. - :rtype: :class:`matplotlib.axes.Axes` - - Example usage: - >>> import matplotlib.pyplot as plt - >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) - >>> plt.figure(figsize=(8, 6)) - >>> prody.showPairwiseRMSDHeatmap(distance_matrix) - >>> plt.show() - """ - - import matplotlib.pyplot as plt - - - if distance_matrix.ndim != 2: - raise ValueError(f"distance_matrix must be a 2D array, but got shape {distance_matrix.shape}.") - if distance_matrix.shape[0] != distance_matrix.shape[1]: - raise ValueError("distance_matrix must be square.") - - if ax is None: - ax = plt.gca() - - im = ax.imshow(distance_matrix, cmap = cmap) - ax.figure.colorbar(im, ax = ax, label = label) - - ax.set_xlabel("# Frame") - ax.set_ylabel("# Frame") - ax.set_title(title) - - return ax - - -def showRMSDHistogram(rmsd_data, bins = 50, title = 'Distribution of RMSDs', - xlabel = 'RMSD [Å]', ylabel = 'Frequency', label = None, - kde = False, element = 'bars', stat = 'count', - color = 'teal', edgecolor = 'black', lw = 0.8, alpha = 0.5, ax = None): - """ - Plots the distribution of RMSD values on the current axis. - The input may be either a one-dimensional RMSD array (e.g., frame-to-reference or intra-cluster RMSDs) - or a 2D pairwise RMSD distance matrix. - - - :arg rmsd_data: 1D array of RMSDs or 2D pairwise distance matrix. - If 2D, the upper triangle (excluding diagonal) is automatically extracted. - Recommended to generate them using :func:`calcRMSDfromReference` or - :func:`calcPairwiseRMSD` respectively. - :type rmsd_data: :class:`numpy.ndarray` - - :arg bins: number of histogram bins. - Default is ``50`` - :type bins: int - - :arg title: the title of the generated plot. - Default is ``'Distribution of RMSDs'`` - :type title: str - - :arg xlabel: the label for the X axis. - Default is ``'RMSD [Å]'`` - :type xlabel: str - - :arg ylabel: the label for the Y axis. - Default is ``'Frequency'`` - :type ylabel: str - - :arg label: the label for the plot legend. - Default is ``None`` - :type label: str - - :arg kde: whether to plot a kernel density estimate. - Default is ``False`` - :type kde: bool - - :arg element: visual representation of the histogram bins ('bars', 'step', or 'poly'). - Default is ``'bars'`` - :type element: str - - :arg stat: aggregate statistic to compute in each bin ('count', 'frequency', 'probability', 'percent', 'density'). - Default is ``'count'`` - :type stat: str - - :arg color: the color of the plot. - Default is ``'teal'`` - :type color: str - - :arg edgecolor: the color of the bin edges. - Default is ``'black'`` - :type edgecolor: str - - :arg lw: the line width of the bin edges. - Default is ``0.8`` - :type lw: float - - :arg alpha: the transparency of the bins. - Default is ``0.5`` - :type alpha: float - - :arg ax: axes on which to draw the plot. - Default is ``None`` and the current axes are used. - :type ax: :class:`matplotlib.axes.Axes` - - :returns: the Matplotlib axes containing the plot. - :rtype: :class:`matplotlib.axes.Axes` - - Example usage: - >>> import matplotlib.pyplot as plt - >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) - >>> plt.figure(figsize = (8, 6)) - >>> prody.showRMSDHistogram(distance_matrix, bins=70, kde=True) - >>> plt.show() - """ - - import matplotlib.pyplot as plt - try: - import seaborn as sns - except ImportError: - raise ImportError("The 'seaborn' package is required to display the Histogram." - "\nPlease install it using 'pip install seaborn'." - "\nAlternatively, use standard matplotlib.pyplot.hist() for basic plots.") - - - data_array = np.asarray(rmsd_data) - - if data_array.ndim == 2: - if data_array.shape[0] != data_array.shape[1]: - raise ValueError(f"Pairwise RMSD matrix must be square, but got shape {data_array.shape}") - rmsd_values = data_array[np.triu_indices_from(data_array, k=1)] - elif data_array.ndim == 1: - rmsd_values = data_array - else: - raise ValueError(f"Expected 1D or 2D array, but got shape {data_array.shape}") - - if ax is None: - ax = plt.gca() - - sns.histplot(rmsd_values, bins = bins, element = element, stat = stat, kde = kde, - alpha = alpha, color = color, edgecolor = edgecolor, linewidth = lw, - label = label, ax = ax) - - ax.set_title(title) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - - if label is not None: - ax.legend() - - ax.grid(axis = 'y', alpha = 0.3) - - return ax - - -def calcClusterPopulations(cluster_ids): - """ - Uses the cluster IDs array, assigning each frame to a cluster, to calculate the population - of each cluster and its corresponding percentage. - - - :arg cluster_ids: a one-dimensional array matching each frame to a cluster. - :type cluster_ids: :class:`numpy.ndarray` - - :returns: a dictionary mapping each cluster ID to its population statistics. - Each value is a dictionary with the keys ``"count"`` and ``"pct"``. - :rtype: dict - - Example usage: - >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) - >>> populations = prody.calcClusterPopulations(cluster_ids) - >>> print(f"In cluster 1: {populations[1]['pct']:.2f}% of the objects") - """ - - cluster_ids = np.asarray(cluster_ids) - if cluster_ids.ndim != 1: - raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}") - if cluster_ids.size == 0: - raise ValueError("cluster_ids is empty.") - - total_frames = len(cluster_ids) - clusters, frequencies = np.unique(cluster_ids, return_counts=True) - - population_data = {} - for c, n in zip(clusters, frequencies): - population_data[int(c)] = {'count': int(n), 'pct': (n / total_frames) * 100} - - return population_data - - -def getCluster(cluster_ids, cluster_number): - """ - Returns the indices of the frames grouped in a specific cluster. - - - :arg cluster_ids: a one-dimensional array assigning each frame to a cluster. - :type cluster_ids: :class:`numpy.ndarray` - - :arg cluster_number: the ID of the cluster whose members we want to return. - :type cluster_number: int - - :returns: an array of the trajectory indices assigned to the specified cluster. - :rtype: :class:`numpy.ndarray` - - Example usage: - >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) - >>> cluster1_indices = prody.getCluster(cluster_ids, 1) - """ - - cluster_ids = np.asarray(cluster_ids) - if cluster_ids.ndim != 1: - raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") - if cluster_ids.size == 0: - raise ValueError("cluster_ids is empty.") - - if isinstance(cluster_number, bool) or not isinstance(cluster_number, (int, np.integer)): - raise TypeError(f"cluster_number must be an integer, but got {type(cluster_number).__name__}.") - - cluster_indices = np.where(cluster_ids == cluster_number)[0] - - if cluster_indices.size == 0: - raise ValueError(f'No frames belong to Cluster {cluster_number}') - - return cluster_indices - - -def getClusterMedoid(distance_matrix, cluster_indices = None, cluster_ids = None, cluster_number = None): - """ - Determines the medoid (representative) of the cluster from the minimum total pairwise - distance to all other cluster members. - Returns a dictionary containing the medoid index within the cluster ("local") and in - the full trajectory ("global"). - - The cluster can be specified either by: - 1. cluster_indices (recommended and can be generated by :func:`getCluster`) or - 2. cluster_ids and cluster_number - - - :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. - :type distance_matrix: :class:`numpy.ndarray` - - :arg cluster_indices: a one-dimensional array of the indices of the cluster members - Default is ``None`` - :type cluster_indices: :class:`numpy.ndarray` - - :arg cluster_ids: a one-dimensional array of the cluster IDs per element - Default is ``None`` - :type cluster_ids: :class:`numpy.ndarray` - - :arg cluster_number: the ID (number) of the cluster we want to investigate - Default is ``None`` - :type cluster_number: int - - :returns: a dictionary of the medoid index, with keys "local" for the index within the - cluster and "global" for the index within the full trajectory - :rtype: dict - - Example usage: - >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) - >>> cluster1_indices = prody.getCluster(cluster_ids, 1) - >>> cluster1_medoid = prody.getClusterMedoid(distance_matrix, cluster1_indices) - >>> print(f"Medoid of cluster 1 is # {cluster1_medoid['local']} in the cluster.") - """ - - cluster_indices = _resolveClusterIndices(cluster_indices=cluster_indices, - cluster_ids=cluster_ids, - cluster_number=cluster_number) - - distance_matrix = _validateDistanceMatrix(distance_matrix) - - n_frames = distance_matrix.shape[0] - if np.any((cluster_indices < 0) | (cluster_indices >= n_frames)): - raise ValueError("cluster_indices contain indices outside the bounds of the distance matrix.") - - cluster_distance_matrix = distance_matrix[np.ix_(cluster_indices, cluster_indices)] - - sum_dist = cluster_distance_matrix.sum(axis=1) - medoid_local = int(np.argmin(sum_dist)) - medoid_global = int(cluster_indices[medoid_local]) - - return {"local": medoid_local, "global": medoid_global} - - -def calcClusterStatistics(distance_matrix, cluster_indices = None, cluster_medoid = None, - cluster_ids = None, cluster_number = None): - """ - Calculates descriptive statistics for a single cluster. - Returns a dictionary with the statistics quantity as a key. - - The cluster can be specified either by: - 1. cluster_indices (recommended and can be generated by :func:`getCluster`) or - 2. cluster_ids and cluster_number - - - :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. - :type distance_matrix: :class:`numpy.ndarray` - - :arg cluster_indices: a one-dimensional array of the indices of the cluster members - Default is ``None`` - :type cluster_indices: :class:`numpy.ndarray` - - :arg cluster_medoid: the dictionary of the cluster medoid. - Expects the 'cluster_medoid' dictionary generated by the function - :func:`getClusterMedoid`. - Default is ``None`` - :type cluster_medoid: dict - - :arg cluster_ids: a one-dimensional array of the cluster IDs per element - Default is ``None`` - :type cluster_ids: :class:`numpy.ndarray` - - :arg cluster_number: the ID (number) of the cluster we want to investigate - Default is ``None`` - :type cluster_number: int - - :returns: a dictionary of statistics quantities with their names as keys - :rtype: dict - - Example usage: - >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) - >>> cluster1_indices = prody.getCluster(cluster_ids, 1) - >>> cluster_stats1 = prody.calcClusterStatistics(distance_matrix, cluster1_indices) - >>> print(f"Mean distance from cluster Medoid: {cluster_stats1['mean']:.2f} ± {cluster_stats1['std']:.2f} [Å]") - """ - - cluster_indices = _resolveClusterIndices(cluster_indices, cluster_ids, cluster_number) - population = len(cluster_indices) - - distance_matrix = _validateDistanceMatrix(distance_matrix) - - if cluster_medoid is None: - cluster_medoid = getClusterMedoid(distance_matrix, cluster_indices=cluster_indices) - else: - if not isinstance(cluster_medoid, dict): - raise TypeError("cluster_medoid must be a dictionary with keys 'global' and 'local'.") - - required_keys = {"global", "local"} - if not required_keys.issubset(cluster_medoid.keys()): - missing = required_keys - cluster_medoid.keys() - raise ValueError(f"Expected keys are 'global' and 'local'.\ncluster_medoid is missing required key(s): {missing}.") - - if cluster_medoid["global"] not in cluster_indices: - raise ValueError("Index not found.\ncluster_medoid does not belong to the specified cluster.") - - if not (0 <= cluster_medoid["local"] < population): - raise ValueError("Local index is out of bounds.") - - if cluster_indices[cluster_medoid["local"]] != cluster_medoid["global"]: - raise ValueError("cluster_medoid['local'] and cluster_medoid['global'] are inconsistent.") - - distance_to_medoid = distance_matrix[cluster_indices, cluster_medoid["global"]] - total_frames = distance_matrix.shape[0] - - stats = { - "population" : population, - "pct" : (population / total_frames) * 100, - "medoid_global" : cluster_medoid["global"], - "medoid_local" : cluster_medoid["local"], - "distances" : distance_to_medoid, - "mean" : np.mean(distance_to_medoid), - "std" : np.std(distance_to_medoid), - "median" : np.median(distance_to_medoid), - "iqr" : np.percentile(distance_to_medoid, 75) - np.percentile(distance_to_medoid, 25), - "p95" : np.percentile(distance_to_medoid, 95), - "max" : np.max(distance_to_medoid) - } - - if cluster_number is not None: - stats["cluster"] = cluster_number - - return stats - - -def calcAllClusterStatistics(distance_matrix, cluster_ids): - """ - Calculates descriptive statistics for all clusters. - Returns a list of dictionaries, each corresponding to a cluster with the statistic quantities - as keys. - - - :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. - :type distance_matrix: :class:`numpy.ndarray` - - :arg cluster_ids: a one-dimensional array of the cluster IDs per element - :type cluster_ids: :class:`numpy.ndarray` - - :returns: a list of dictionaries, one dictionary for each cluster with the statistic - quantities as keys. - :rtype: list of dict - - Example usage: - >>> cluster_ids, _= prody.clusterHierarchical(distance_matrix, 4) - >>> all_stats = prody.calcAllClusterStatistics(distance_matrix, cluster_ids) - """ - - cluster_ids = np.asarray(cluster_ids) - if cluster_ids.ndim != 1: - raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") - - distance_matrix = _validateDistanceMatrix(distance_matrix) - - if len(cluster_ids) != distance_matrix.shape[0]: - raise ValueError(f"Dimension mismatch: cluster_ids length ({len(cluster_ids)}) " - f"does not match distance_matrix frames ({distance_matrix.shape[0]}).") - - clusters = np.unique(cluster_ids) - all_stats = [] - - for c in clusters: - # Ignore noise - if c <= 0: - continue - - cluster_stats = calcClusterStatistics(distance_matrix, cluster_ids = cluster_ids, cluster_number = int(c)) - all_stats.append(cluster_stats) - - return all_stats - - -def _resolveClusterIndices(cluster_indices, cluster_ids, cluster_number): - """ - Checks whether the input is cluster_indices or cluster_ids and cluster_number, and for the latter returns - the cluster indices. - """ - - has_indices = cluster_indices is not None - has_ids_info = (cluster_ids is not None) and (cluster_number is not None) - - if has_indices == has_ids_info: - raise ValueError("Provide exactly one of 'cluster_indices' or both 'cluster_ids' and 'cluster_number'.") - - if has_indices: - cluster_indices = np.asarray(cluster_indices) - - if cluster_indices.ndim != 1: - raise ValueError(f"cluster_indices must be a 1D array, but shape {cluster_indices.shape} was given.") - if cluster_indices.size == 0: - raise ValueError("cluster_indices must not be empty.") - else: - cluster_ids = np.asarray(cluster_ids) - cluster_indices = getCluster(cluster_ids, cluster_number) - - return cluster_indices - - -def _validateDistanceMatrix(distance_matrix): - """ - Validates that a distance matrix is a non-empty square 2D NumPy array. - """ - - distance_matrix = np.asarray(distance_matrix) - if distance_matrix.ndim != 2: - raise ValueError(f"The distance matrix must be a 2D array, but got shape {distance_matrix.shape}.") - if distance_matrix.shape[0] != distance_matrix.shape[1]: - raise ValueError(f"distance_matrix must be square, but got shape {distance_matrix.shape}.") - if distance_matrix.size == 0: - raise ValueError("distance_matrix cannot be empty.") - - return distance_matrix - - -def showClusterStatisticsTable(all_stats, dissimilarity = "RMSD", units = "Å", show = True): - """ - Prints the cluster statistics in a table, where each column corresponds to a cluster. - Expects either the 'stats' dictionary generated by :func:`calcClusterStatistics` or - the 'all_stats' list of dictionaries generated by :func:`calcAllClusterStatistics`. - - - :arg all_stats: a dictionary, or list of dictionaries, with the cluster descriptive statistics - Recommended to generate the dictionaries from the functions :func:`calcClusterStatistics` - of a single cluster or :func:`calcAllClusterStatistics`. - :type all_stats: dict, or list of dict - - :arg dissimilarity: the dissimilarity measure - Default is `"RMSD"` - :type dissimilarity: str - - :arg units: the units of the dissimilarity measure - Default is `"Å"` - :type units: str - - :arg show: whether to print the table - Default is `True` - :type show: bool - - :returns: a dictionary containing the table, the row labels and the headers, with their names as keys - :rtype: dict - - Example usage: - >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) - >>> all_stats = prody.calcAllClusterStatistics(distance_matrix, cluster_ids) - >>> prody.showClusterStatisticsTable(all_stats); - """ - - try: - from tabulate import tabulate - except ImportError: - raise ImportError("The 'tabulate' package is required to display the table. " - "Please install it using 'pip install tabulate'.") - - - if isinstance(all_stats, dict): - all_stats = [all_stats] - elif not isinstance(all_stats, list): - raise TypeError(f"all_stats must be a dict or list of dicts, but got {type(all_stats).__name__}.") - - if not all_stats: - raise ValueError("all_stats cannot be empty.") - - required = {"population", "pct", "medoid_global", "medoid_local", "mean", "std", - "median", "iqr", "p95", "max"} - - headers = [] - for i, stat in enumerate(all_stats): - if not isinstance(stat, dict): - raise TypeError(f"Expected a dictionary in all_stats, but got {type(stat).__name__} at index {i}.") - - cluster_id = stat.get('cluster', i + 1) - headers.append(f"Cluster {int(cluster_id)}") - - missing = required - stat.keys() - if missing: - raise ValueError(f"Statistics dictionary is missing required keys: {sorted(missing)}.") - - metrics = [ - ("population", "Total Frames"), - ("pct", "Population Percentage (%)"), - ("medoid_global", "Medoid Frame (Global)"), - ("medoid_local", "Medoid Frame (Within Cluster)"), - ("mean", f"Mean {dissimilarity} [{units}]"), - ("std", f"Std {dissimilarity} [{units}]"), - ("median", f"Median {dissimilarity} [{units}]"), - ("iqr", f"IQR [{units}]"), - ("p95", f"95th Percentile [{units}]"), - ("max", f"Max {dissimilarity} [{units}]") - ] - - table = [] - row_labels = [] - - for key, label in metrics: - row_labels.append(label) - table.append([stat[key] for stat in all_stats]) - - if not isinstance(show, bool): - raise TypeError("show must be a bool.") - - if show: - print(tabulate(table, headers=headers, showindex=row_labels, - tablefmt="fancy_grid", floatfmt=".4f", stralign='center')) - - return {"table": table, "row_labels": row_labels, "headers": headers} - - -def showClusterRMSDComparison(all_stats, bins = 50, title = 'RMSD Distributions', - xlabel = 'RMSD to Medoid [Å]', ylabel = 'Frequency', - kde = True, element = 'step', stat = 'count', - alpha = 0.5, lw = 1.5, ax = None): - """ - Overlays the internal RMSD distributions of all clusters onto a single plot - on the current axis. - - Expects the list of cluster statistics dictionaries generated by :func:`calcAllClusterStatistics`. - - - :arg all_stats: list of cluster statistics dictionaries. Each dictionary - must contain the keys ``"cluster"`` and ``"distances"``. - :type all_stats: list of dict - - :arg bins: number of histogram bins. - Default is ``50`` - :type bins: int - - :arg title: the title of the generated plot. - Default is ``'Intra-Cluster RMSD Distributions'`` - :type title: str - - :arg xlabel: the label for the X axis. - Default is ``'RMSD to Medoid [Å]'`` - :type xlabel: str - - :arg ylabel: the label for the Y axis. - Default is ``'Frequency'`` - :type ylabel: str - - :arg kde: whether to plot kernel density estimates. - Default is ``True`` - :type kde: bool - - :arg element: visual representation of the histogram bins ('step', 'bars', or 'poly'). - Default is ``'step'`` - :type element: str - - :arg stat: aggregate statistic to compute in each bin. - Default is ``'count'`` - :type stat: str - - :arg alpha: the transparency of the bins. - Default is ``0.5`` - :type alpha: float - - :arg lw: line width of the distribution lines. - Default is ``1.5`` - :type lw: float - - :arg ax: axes on which to draw the plot. - Default is ``None`` and the current axes are used. - :type ax: :class:`matplotlib.axes.Axes` - - :returns: the Matplotlib axes containing the plot. - :rtype: :class:`matplotlib.axes.Axes` - - Example usage: - >>> import matplotlib.pyplot as plt - >>> all_stats = prody.calcAllClusterStatistics(distance_matrix, clusterIDs_array) - >>> plt.figure(figsize = (8, 6)) - >>> prody.showClusterRMSDComparison(all_stats) - >>> plt.show() - """ - - import matplotlib.pyplot as plt - try: - import seaborn as sns - except ImportError: - raise ImportError("The 'seaborn' package is required to display the Histogram." - "\nPlease install it using 'pip install seaborn'." - "\nAlternatively, use standard matplotlib.pyplot.hist() for basic plots.") - - - if len(all_stats) == 0: - raise ValueError("all_stats cannot be empty.") - - if ax is None: - ax = plt.gca() - - for cluster_stats in all_stats: - sns.histplot(cluster_stats["distances"], bins = bins, element = element, stat = stat, kde = kde, - alpha = alpha, linewidth = lw, label = f"Cluster {cluster_stats['cluster']}", - ax = ax) - - ax.set_title(title) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - - ax.legend(title = "Clusters") - ax.grid(axis = 'y', alpha = 0.3) - - return ax - - -def clusterHierarchical(distance_matrix, method='average', cutoff=None): - """ - Performs bottom-up hierarchical clustering from a pairwise distance matrix. - - Note that the resulting cluster IDs are 1-based (starting at 1, not 0). - - If cutoff is ``None``, it is automatically chosen as the midpoint of the largest - gap between consecutive linkage distances. - - Recommendation: Plot the dendrogram using showDendrogram() and choose the cutoff - manually whenever possible. - - - :arg distance_matrix: either a one-dimensional condensed distance matrix or - a two-dimensional pairwise distance matrix. - :type distance_matrix: :class:`numpy.ndarray` - - :arg method: linkage criterion used when constructing the linkage matrix. - Default is ``'average'`` - :type method: str - - :arg cutoff: the cutoff distance that determines the number of clusters. - Default is ``None`` and the midpoint of the largest linkage gap - is used automatically. - :type cutoff: float - - :returns: - * a one-dimensional array containing the cluster ID for each object. - IDs are 1-indexed (1 to number of clusters for the cutoff) - * the hierarchical clustering linkage matrix - :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) - - Example usage: - >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) - >>> cluster_ids, linkage_matrix = prody.clusterHierarchical(distance_matrix) - """ - - from scipy.cluster.hierarchy import linkage, fcluster - - - condensed_distance_matrix = _condenseDistanceMatrix(distance_matrix) - linkage_matrix = linkage(condensed_distance_matrix, method=method) - - if cutoff is None: - cutoff = _calcAutoCutoff(linkage_matrix) - elif isinstance(cutoff, str): - raise ValueError(f"Invalid input for cutoff: '{cutoff}'. Use numeric value or None.") - - cluster_ids = fcluster(linkage_matrix, t=cutoff, criterion='distance') - - return cluster_ids, linkage_matrix - - -def showDendrogram(distance_matrix = None, linkage_matrix = None, method = 'average', cutoff = None, - truncate_mode = None, p = 30, title = "Hierarchical Clustering Dendrogram", - ylabel = "RMSD [Å]", ax = None): - """ - Plots a hierarchical clustering dendrogram using either a pairwise RMSD distance matrix - or a pre-computed linkage matrix. - - By default, the complete dendrogram is shown. To display a truncated dendrogram, - use ``truncate_mode`` together with ``p``. - - Exactly one of ``distance_matrix`` or ``linkage_matrix`` must be provided. - - - :arg distance_matrix: one-dimensional condensed or two-dimensional square pairwise - distance matrix. Typically generated using :func:`calcPairwiseRMSD`. - Default is ``None`` - :type distance_matrix: :class:`numpy.ndarray` - - :arg linkage_matrix: pre-computed hierarchical clustering linkage matrix. - Default is ``None`` - :type linkage_matrix: :class:`numpy.ndarray` - - :arg method: linkage criterion used when constructing the linkage matrix from - ``distance_matrix``. - Default is ``"average"`` - :type method: str - - :arg cutoff: dendrogram color threshold. If ``"auto"``, the cutoff is placed at the midpoint - of the largest linkage gap. - Default is ``None`` - :type cutoff: float, str - - :arg truncate_mode: dendrogram truncation mode passed to :func:`scipy.cluster.hierarchy.dendrogram`. - Default is ``None`` - :type truncate_mode: str, None - - :arg p: truncation parameter used together with ``truncate_mode``. - Default is ``30`` - :type p: int - - :arg title: title of the generated plot. - Default is ``"Hierarchical Clustering Dendrogram"`` - :type title: str - - :arg ylabel: Distance matric label for the y-axis - Default is ``"RMSD [Å]"`` - :type ylabel: str - - :arg ax: axes on which to draw the plot. - Default is ``None`` and the current axes are used. - :type ax: :class:`matplotlib.axes.Axes` - - :returns: the matplotlib axes and the linkage matrix. - :rtype: tuple(:class:`matplotlib.axes.Axes`, :class:`numpy.ndarray`) - - Example usage: - >>> import matplotlib.pyplot as plt - >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) - >>> plt.figure(figsize=(12, 8)) - >>> ax, linkage = prody.showDendrogram(distance_matrix, cutoff='auto', truncate_mode='lastp', p=30) - >>> plt.show() - """ - - import matplotlib.pyplot as plt - from scipy.cluster.hierarchy import dendrogram, linkage - - - if (distance_matrix is None) == (linkage_matrix is None): - raise ValueError("Provide exactly one of 'distance_matrix' or 'linkage_matrix'.") - - if linkage_matrix is None: - condensed_distance_matrix = _condenseDistanceMatrix(distance_matrix) - linkage_matrix = linkage(condensed_distance_matrix, method=method) - else: - linkage_matrix = np.asarray(linkage_matrix) - - if linkage_matrix.ndim != 2 or linkage_matrix.shape[1] != 4 or linkage_matrix.shape[0] == 0: - raise ValueError("linkage_matrix must be a non-empty array with shape (n_samples-1, 4).") - - if cutoff == "auto": - cutoff = _calcAutoCutoff(linkage_matrix) - elif isinstance(cutoff, str): - raise ValueError(f"Invalid string for cutoff: '{cutoff}'. Use 'auto' or a numeric value.") - - if ax is None: - ax = plt.gca() - - with plt.rc_context({"lines.linewidth": 0.6}): - dendrogram(linkage_matrix, truncate_mode = truncate_mode, p = p, color_threshold = cutoff, - no_labels = True, ax = ax) - - ax.set_xlabel("# Frame") - ax.set_ylabel(ylabel) - ax.set_title(title) - - if cutoff is not None: - ax.axhline(cutoff, color = "black", linestyle = "--", linewidth = 1.2, label = f"Cutoff ({cutoff:.2f} Å)") - ax.legend() - - return ax, linkage_matrix - - -def _calcAutoCutoff(linkage_matrix): - """ - Calculates the optimal distance cutoff based on the midpoint of the largest gap - in the linkage matrix. - """ - - linkage_matrix = np.asarray(linkage_matrix) - merge_dist = linkage_matrix[:, 2] # locations of the merges - - if len(merge_dist) == 0: - raise ValueError("At least two frames are required for clustering") - elif len(merge_dist) == 1: - return merge_dist[0] + 1.0 - - gaps = np.diff(merge_dist) - largest_gap_idx = np.argmax(gaps) - - return (merge_dist[largest_gap_idx] + merge_dist[largest_gap_idx + 1]) / 2.0 - - -def _condenseDistanceMatrix(distance_matrix): - """ - Converts a pairwise distance matrix into condensed form. - - Accepts either: - * a condensed one-dimensional distance matrix, or - * a square two-dimensional distance matrix. - - Returns the condensed distance matrix suitable for scipy.cluster.hierarchy.linkage. - """ - from scipy.spatial.distance import squareform - - - distance_matrix = np.asarray(distance_matrix) - - if distance_matrix.size == 0: - raise ValueError("distance_matrix cannot be empty.") - - if distance_matrix.ndim == 2: - if distance_matrix.shape[0] != distance_matrix.shape[1]: - raise ValueError(f"distance_matrix must be square, but got shape {distance_matrix.shape}.") - condensed_distance_matrix = squareform(distance_matrix) - elif distance_matrix.ndim == 1: - condensed_distance_matrix = distance_matrix - else: - raise ValueError("distance_matrix must be either a condensed 1D array or a square 2D array.") - - return condensed_distance_matrix - - -def clusterKMedoids(distance_matrix, k, method = 'alternate', method_sklearn = 'pam', - initial_medoids = None, seed = None, max_iter = 100, n_init = 10): - """ - Performs K-Medoids clustering using various algorithms. - - This function acts as a facade, routing the clustering task to the specified backend - ('alternate', 'pam', or 'sklearn'). 'alternate' is generally faster, while 'pam' typically - finds solutions with lower cost. 'sklearn' needs the sklearn_extra.cluster module to be installed. - - Note that the resulting cluster IDs are 1-based (starting at 1, not 0). - - - :arg distance_matrix: square, symmetric matrix of pairwise distances between all objects - :type distance_matrix: :class:`numpy.ndarray` - - :arg k: prespecified number of clusters to form - :type k: int - - :arg method: the clustering algorithm to use. - Options are 'alternate' (custom Alternating Medoids), 'pam' (custom PAM), 'sklearn' (sklearn_extra) - Default is 'alternate' - :type method: str - - :arg method_sklearn: the specific method to pass to sklearn_extra if method = 'sklearn' is chosen. - Options are 'pam', 'alternate' - Default is 'pam' - :type method_sklearn: str - - :arg initial_medoids: one-dimensional array of indices to use as starting medoids. - Default is ``None`` and the starting medoids are picked randomly. - It is not supported for the 'sklearn' method. - :type initial_medoids: :class:`numpy.ndarray` - - :arg seed: random seed for reproducibility - Default is ``None`` - :type seed: int - - :arg max_iter: maximum number of iterations per a single run - Default is ``100`` - :type max_iter: int - - :arg n_init: number of times the algorithm will be run with different initial medoids - and the best result with the lowest cost is returned - - :returns: a tuple of: - * an one-dimensional array containing the cluster ID for each object. - IDs are 1-indexed (1 to k) - * a one-dimensional array of shape (k,) containing the indices of the final cluster - medoids - * the final sum of distances from each point to its nearest medoid - :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`, float) - - Example usage: - >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) - >>> cluster_ids, medoids, _ = prody.clusterKMedoids(distance_matrix, 4, method = 'alternate', - seed = 42, n_init = 30) - """ - - distance_matrix = _validateDistanceMatrix(distance_matrix) - n = distance_matrix.shape[0] - - if not isinstance(k, (int, np.integer)): - raise TypeError(f"k must be an integer, got {type(k).__name__}") - if k <= 0 or k > n: - raise ValueError(f"k must be between 1 and {n} (the number of points).") - - if not isinstance(max_iter, (int, np.integer)) or max_iter <= 0: - raise ValueError("max_iter must be a positive integer.") - - if not isinstance(n_init, (int, np.integer)) or n_init <= 0: - raise ValueError("n_init must be a positive integer.") - - if seed is not None and not isinstance(seed, (int, np.integer)): - raise TypeError("seed must be an integer or None.") - - if initial_medoids is not None: - if method == 'sklearn': - raise ValueError("The 'sklearn' method does not support custom 'initial_medoids'. Use 'alternate' or 'pam' instead.") - - initial_medoids = np.asarray(initial_medoids, dtype=int) - if initial_medoids.ndim != 1: - raise ValueError("initial_medoids must be a 1D array.") - if len(initial_medoids) != k: - raise ValueError(f"Number of initial medoids ({len(initial_medoids)}) must equal k ({k}).") - if len(np.unique(initial_medoids)) != k: - raise ValueError("initial_medoids must be unique.") - if np.any((initial_medoids < 0) | (initial_medoids >= n)): - raise ValueError(f"initial_medoids contain indices out of bounds (must be 0 to {n-1}).") - - n_init = 1 - - if method == 'alternate': - return _clusterKMedoidsAlternating(distance_matrix, k, initial_medoids, seed, max_iter, n_init) - elif method == 'pam': - return _clusterKMedoidsPAM(distance_matrix, k, initial_medoids, seed, max_iter, n_init) - elif method == 'sklearn': - if method_sklearn not in ['alternate', 'pam']: - raise ValueError(f"method_sklearn '{method_sklearn}' is not valid. Options: 'alternate', 'pam'") - return _clusterKMedoidsSklearn(distance_matrix, k, method_sklearn, seed, max_iter, n_init) - else: - raise ValueError(f"Method '{method}' is not valid. Options: 'alternate', 'pam', 'sklearn'") - - -def _cost(distance_matrix, medoids): - """Calculates the sum of distances from all points to their nearest medoid.""" - distances = distance_matrix[:, medoids] - closest_distance = np.min(distances, axis=1) - return np.sum(closest_distance) - - -def _clusterKMedoidsAlternating(distance_matrix, k, initial_medoids, seed, max_iter, n_init): - """K-Medoids clustering with an Alternating Medoids algorithm""" - - n = distance_matrix.shape[0] - rng = np.random.default_rng(seed) - - best_cost = np.inf - best_medoids = None - - for run in range(n_init): - if initial_medoids is None: - medoids = np.asarray(rng.choice(n, size = k, replace = False)) - else: - medoids = initial_medoids.copy() - - for iteration in range(max_iter): - distances_to_medoids = distance_matrix[:, medoids] - cluster_ids = np.argmin(distances_to_medoids, axis = 1) - - new_medoids = np.zeros_like(medoids) - - for i in range(len(medoids)): - cluster_members = np.where(cluster_ids == i)[0] - - # Handles empty clusters by retaining the old medoid - if len(cluster_members) == 0: - new_medoids[i] = medoids[i] - continue - - intra_cluster_distances = distance_matrix[np.ix_(cluster_members, cluster_members)] - sum_distances = intra_cluster_distances.sum(axis=1) - best_medoid_idx_in_cluster = np.argmin(sum_distances) - new_medoids[i] = cluster_members[best_medoid_idx_in_cluster] - - current_cost = _cost(distance_matrix, new_medoids) - - if np.array_equal(medoids, new_medoids): - break - - medoids = new_medoids - - if current_cost < best_cost: - best_cost = current_cost - best_medoids = medoids.copy() - - distances = distance_matrix[:, best_medoids] - cluster_ids = np.argmin(distances, axis=1) + 1 - - return cluster_ids, best_medoids, best_cost - - -def _clusterKMedoidsPAM(distance_matrix, k, initial_medoids, seed, max_iter, n_init): - """K-Medoids clustering with PAM (Partitioning Around Medoids)""" - - n = distance_matrix.shape[0] - rng = np.random.default_rng(seed) - - best_cost = np.inf - best_medoids = None - - for run in range(n_init): - if initial_medoids is None: - medoids = np.asarray(rng.choice(n, size=k, replace=False)) - else: - medoids = initial_medoids.copy() - - # Using a set for removal/addition operations - non_medoids = set(i for i in range(n) if i not in medoids) - current_cost = _cost(distance_matrix, medoids) - - for iteration in range(max_iter): - best_cost_swap = current_cost - best_swap = None - - for medoid_idx, old_medoid in enumerate(medoids): - for new_medoid in non_medoids: - candidate = medoids.copy() - candidate[medoid_idx] = new_medoid - candidate_cost = _cost(distance_matrix, candidate) - - if candidate_cost < best_cost_swap: - best_cost_swap = candidate_cost - best_swap = (medoid_idx, old_medoid, new_medoid) - - if best_swap is None: - break - - idx, old_medoid, new_medoid = best_swap - medoids[idx] = new_medoid - non_medoids.remove(new_medoid) - non_medoids.add(old_medoid) - current_cost = best_cost_swap - - if current_cost < best_cost: - best_cost = current_cost - best_medoids = medoids.copy() - - distances = distance_matrix[:, best_medoids] - cluster_ids = np.argmin(distances, axis=1) + 1 - - return cluster_ids, best_medoids, best_cost - - -def _clusterKMedoidsSklearn(distance_matrix, k, method_sklearn, seed, max_iter, n_init): - """K-Medoids clustering with sklearn_extra""" - try: - from sklearn_extra.cluster import KMedoids - except ImportError: - raise ImportError("The 'sklearn_extra' package is required for this K-Medoids approach. " - "Please install it using 'pip install scikit-learn-extra'.") - - - best_cost = np.inf - best_labels = None - best_medoids = None - - rng = np.random.default_rng(seed) - - for run in range(n_init): - kmedoids = KMedoids(n_clusters=k, metric='precomputed', method=method_sklearn, init='random', - max_iter=max_iter, random_state=int(rng.integers(0, 1000000))) - - kmedoids.fit(distance_matrix) - medoids = kmedoids.medoid_indices_ - - cost = np.sum(np.min(distance_matrix[:, medoids], axis=1)) - - if cost < best_cost: - best_cost = cost - best_labels = kmedoids.labels_ + 1 - best_medoids = medoids.copy() - - return best_labels, best_medoids, best_cost - - -def clusterDBSCAN(distance_matrix, eps = None, minPts = None, method='custom'): - """ - Performs DBSCAN clustering using various algorithms. - - This function acts as a facade, routing the clustering task to the specified backend - ('custom' or 'sklearn'). 'custom' uses a built-in implementation without additional - dependencies beyond NumPy. 'sklearn' needs the sklearn.cluster module to be installed. - - Note that the resulting cluster IDs are 1-based (starting at 1, not 0). - Noise points are labeled -1. - - - :arg distance_matrix: square, symmetric matrix of pairwise distances between all objects - :type distance_matrix: :class:`numpy.ndarray` - - :arg eps: the "radius" of the neighborhood within which we count neighbors - Default is ``None`` and automatically the median of the distances - in the distance matrix is used. - Ideally use :func:`showReachabilityPlot` to determine - manually the most suitable eps. - :type eps: float - - :arg minPts: the minimum number of neighbors required for a point to be considered - a core point - Default is ``None`` and automatically the 5% of the total objects, - or for less than 20 2 is used. - Ideally choose manually the most suitable minPts - :type minPts: int - - :arg method: the clustering algorithm to use - Options are 'custom' and 'sklearn' - Default is 'custom' because it needs no module installation - :type method: str - - :returns: a tuple of: - * a one-dimensional array containing the cluster ID for each object. - IDs are 1-indexed and noise corresponds to -1. - * a one-dimensional array of the frame indices corresponding to noise - :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) - - Example usage: - >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) - >>> cluster_ids, _ = prody.clusterDBSCAN(distance_matrix, eps = 1.8, minPts = 30) - """ - - distance_matrix = _validateDistanceMatrix(distance_matrix) - - if eps is None: - eps = _calcAutoEps(distance_matrix) - _validateEps(eps, distance_matrix) - - if minPts is None: - minPts = _calcAutoMinPts(distance_matrix) - _validateMinPts(minPts, distance_matrix) - - if method == 'custom': - return _clusterDBSCANCustom(distance_matrix, eps, minPts) - elif method == 'sklearn': - return _clusterDBSCANSklearn(distance_matrix, eps, minPts) - else: - raise ValueError(f"method can be either 'custom' or 'sklearn', but got {method}.") - - -def _clusterDBSCANCustom(distance_matrix, eps, minPts): - """DBSCAN clustering with custom algorithm""" - total_points = distance_matrix.shape[0] - - labels = np.zeros(total_points, dtype = int) - cluster_id = 0 - - for p in range(total_points): - - if labels[p] != 0: - continue - - neighbors, = np.where(distance_matrix[p] <= eps) - - if len(neighbors) < minPts: - labels[p] = -1 - else: - cluster_id += 1 - labels[p] = cluster_id - - candidate_set = [n for n in neighbors if n != p] - while candidate_set: - q = candidate_set.pop() - - if labels[q] == -1: - labels[q] = cluster_id - - if labels[q] != 0: - continue - - labels[q] = cluster_id - - q_neighbors, = np.where(distance_matrix[q] <= eps) - if len(q_neighbors) >= minPts: - for n in q_neighbors: - if labels[n] == 0: - candidate_set.append(n) - elif labels[n] == -1: - labels[n] = cluster_id - - noise_frames, = np.where(labels == -1) - - return labels, noise_frames - - -def _clusterDBSCANSklearn(distance_matrix, eps, minPts): - """DBSCAN clustering with sklearn""" - - try: - from sklearn.cluster import DBSCAN - except ImportError: - raise ImportError("The 'sklearn' package is required for this DBSCAN approach. " - "Please install it using 'pip install scikit-learn'.") - - dbscan = DBSCAN(eps=eps, min_samples = minPts, metric = "precomputed") - labels = dbscan.fit_predict(distance_matrix) - - cluster_ids = np.copy(labels) - cluster_ids[cluster_ids >= 0] += 1 - - noise_frames, = np.where(labels == -1) - - return cluster_ids, noise_frames - - -def showReachabilityPlot(distance_matrix, minPts = None, method = 'custom', eps = None, - title = "OPTICS Reachability Plot", xlabel = "Frames (Sorted by OPTICS)", - ylabel = "Reachability Distance [Å]", label = None, - color = "#36454F", lw = 1.5, fill = True, ax = None): - """ - Plots the reachability plot using a simplified OPTICS algorithm. - The reachability plot should be used to determine the most suitable eps - parameter for DBSCAN. - - This function acts as a facade, routing the ordering task to the specified backend - ('custom' or 'sklearn'). 'custom' uses a built-in implementation without additional - dependencies beyond NumPy. 'sklearn' needs the sklearn.cluster module to be installed. - - - :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. - :type distance_matrix: :class:`numpy.ndarray` - - :arg minPts: the minimum number of neighbors required for a point to be considered - a core point - Default is ``None`` and automatically the 5% of the total objects, - or for less than 20 2 is used. - Ideally choose manually the most suitable minPts - :type minPts: int - - :arg method: the OPTICS algorithm to use - Options are 'custom' and 'sklearn'. - Default is 'custom' because it needs no module installation - :type method: str - - :arg eps: the "radius" of the neighborhood within which we count neighbors. - It is drawn as a horizontal line to help visualize the DBSCAN - ``eps`` threshold. - If `'auto'` the median of the pairwise distances is used. - Default is ``None``. - :type eps: float, str - - :arg title: title of the generated plot. - Default is ``"OPTICS Reachability Plot"`` - :type title: str - - :arg xlabel: Frame count label for the x-axis - Default is ``"Frames (Sorted by OPTICS)"`` - :type xlabel: str - - :arg ylabel: Reachability distance label for the y-axis - Default is ``"Reachability Distance [Å]"`` - :type ylabel: str - - :arg label: the label for the plot. - Default is ``None`` - :type label: str - - :arg color: the color of the plot. - Default is ``"#36454F"`` - :type color: str - - :arg lw: the width of the line. - Default is ``1.5`` - :type lw: float - - :arg fill: whether to fill the area beneath the curve and if applicable - the valleys below the eps threshold. - :type fill: bool - - :arg ax: axes on which to draw the plot. - Default is ``None`` and the current axes are used. - :type ax: :class:`matplotlib.axes.Axes` - - :returns: the Matplotlib axes containing the plot. - :rtype: :class:`matplotlib.axes.Axes` - - Example usage: - >>> import matplotlib.pyplot as plt - >>> plt.figure(figsize = (8, 6)) - >>> showReachabilityPlot(distance_matrix, minPts = 20, eps = 1.8) - >>> plt.show() - """ - - import matplotlib.pyplot as plt - - - distance_matrix = _validateDistanceMatrix(distance_matrix) - - if minPts is None: - minPts = _calcAutoMinPts(distance_matrix) - _validateMinPts(minPts, distance_matrix) - - if not isinstance(fill, bool): - raise TypeError(f"fill must be a bool, but got {type(fill).__name__}") - - if method == 'custom': - reachability, ordering = _orderOPTICSCustom(distance_matrix, minPts) - elif method == 'sklearn': - reachability, ordering = _orderOPTICSSklearn(distance_matrix, minPts) - else: - raise ValueError(f"method must be either 'custom' or 'sklearn', but got {method}") - - if ax is None: - ax = plt.gca() - - ax.set_title(title) - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - - y = reachability[ordering] - x = np.arange(len(y)) - - ax.plot(x, y, color = color, lw = lw, label = label) - - if fill: - ax.fill_between(x, 0, y, color = 'black', alpha = 0.4) - - if eps == 'auto': - eps = _calcAutoEps(distance_matrix) - - if eps is not None: - _validateEps(eps, distance_matrix) - ax.axhline(y = eps, color = 'black', linestyle = '--', linewidth = lw) - - if fill: - below_eps_mask = (y <= eps) - segments = [] - start = None - - for i, below in enumerate(below_eps_mask): - if below and start is None: - start = i - elif not below and start is not None: - segments.append((start, i)) - start = None - - if start is not None: - segments.append((start, len(below_eps_mask))) - - colors = plt.cm.tab10(np.linspace(0, 1, len(segments))) - - for (start, end), valley_color in zip(segments, colors): - ax.fill_between(x[start:end], y[start:end], eps, color=valley_color) - - if label is not None: - ax.legend() - - ax.grid(axis = 'y', linestyle = '--', linewidth = 0.8, alpha = 0.3) - - return ax - - -def _orderOPTICSCustom(distance_matrix, minPts): - """OPTICS algorithm with built-in modules and NumPy""" - - import heapq - - total_points = distance_matrix.shape[0] - - sorted_distances = np.sort(distance_matrix, axis = 1) - core_distances = sorted_distances[:, minPts - 1] - - reachability = np.full(total_points, np.inf) - processed = np.zeros(total_points, dtype = bool) - ordering = [] - - def _updateSeeds(idx): - new_reaches = np.maximum(core_distances[idx], distance_matrix[idx, :]) - update_mask = (~processed) & (new_reaches < reachability) - points_to_update, = np.where(update_mask) - reachability[update_mask] = new_reaches[update_mask] - - for j in points_to_update: - heapq.heappush(seeds, (reachability[j], j)) - - for i in range(total_points): - if processed[i]: - continue - - processed[i] = True - ordering.append(i) - seeds = [] - - _updateSeeds(i) - - while seeds: - current_reach, q = heapq.heappop(seeds) - - if processed[q]: - continue - - processed[q] = True - ordering.append(q) - - _updateSeeds(q) - - return reachability, np.array(ordering) - - -def _orderOPTICSSklearn(distance_matrix, minPts): - """OPTICS algorithm with sklearn""" - - try: - from sklearn.cluster import OPTICS - except ImportError: - raise ImportError("The 'sklearn' package is required for this OPTICS approach. " - "Please install it using 'pip install scikit-learn'.") - - optics = OPTICS(min_samples = minPts, metric = 'precomputed') - optics.fit(distance_matrix) - reachability = optics.reachability_ - ordering = optics.ordering_ - - return reachability, ordering - - -def _calcAutoEps(distance_matrix): - """Automatically determine the DBSCAN eps parameter""" - - import warnings - - pairwise_distances = distance_matrix[np.triu_indices_from(distance_matrix, k = 1)] - eps = float(np.median(pairwise_distances)) - warnings.warn(f"No eps provided. Automatically chosen at {eps:.3f}. Ideally provide your own value.") - return eps - - -def _validateEps(eps, distance_matrix): - """Validate the DBSCAN eps parameter""" - - import warnings - - maxDistance = np.max(distance_matrix) - - if not isinstance(eps, (float, int, np.floating, np.integer)): - raise TypeError(f"eps must be a numeric value, but got {type(eps).__name__}") - - if eps <= 0: - raise ValueError("eps must be positive") - elif eps > maxDistance: - warnings.warn(f"eps ({eps}) is greater than the maximum pairwise distance ({maxDistance:.3f}).\n" - "All frames will be clustered together with no noise.") - - -def _calcAutoMinPts(distance_matrix): - """Automatically determine the minPts parameter""" - - import warnings - - total_points = distance_matrix.shape[0] - minPts = max(2, int(total_points // 20)) - warnings.warn(f"No minPts provided. Automatically chosen at {minPts}. Ideally provide your own value.") - return minPts - - -def _validateMinPts(minPts, distance_matrix): - """Validate minPts parameter""" - total_points = distance_matrix.shape[0] - - if not isinstance(minPts, (int, np.integer)): - raise TypeError(f"minPts must be a positive integer, but got {type(minPts).__name__}") - - if minPts <= 0 or minPts > total_points: - raise ValueError(f"minPts must be between 1 and {total_points}.") - - -def writeClusters(atoms, trajectory, distance_matrix, cluster_ids, write_dcd = True, - align ="protein and backbone", system = "system", tag ="cluster"): - """ - Aligns the trajectory once, then exports representative medoid structures as PDB files and - cluster-specific DCD trajectories. - Returns the name of all exported .pdb and .dcd files. - - Note: This function loads all aligned coordinates into memory. It is highly optimized - for speed, provided the trajectory fits within available system RAM. - - :arg atoms: reference structure used for the alignment. - :type atoms: :class:`prody.Atomic` - - :arg trajectory: trajectory containing the coordinate sets to align - :type trajectory: :class:`prody.Trajectory` - - :arg distance_matrix: square, symmetric matrix of pairwise distances between all objects - :type distance_matrix: :class:`numpy.ndarray` - - :arg cluster_ids: a one-dimensional array of the cluster IDs per element - :type cluster_ids: :class:`numpy.ndarray` - - :arg write_dcd: determines whether to save a .dcd file of the cluster frames - Default is ``True`` - :type write_dcd: bool - - :arg align: atom selection used to calculate the alignment transformation. - Must be a valid ProDy selection string. - Default is ``"protein and backbone"`` - :type align: str - - :arg system: the name of the system under investigation - Default is ``"system"`` - :type system: str - - :arg tag: the name of the clustering method - Default is ``"cluster"`` - :type tag: str - - :returns: list of exported filenames - :rtype: list[str] - - Example usage: - >>> pdb = prody.parsePDB("structure.pdb") - >>> dcd = prody.Trajectory("trajectory.dcd") - >>> _, aligned_coords = prody.alignTrajectory(pdb, dcd, select='resname IOA') - >>> distance_matrix = prody.calcPairwiseRMSD(aligned_coords) - >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix) - >>> prody. writeClusters(pdb, dcd, distance_matrix, cluster_ids, - system="type1_RUN23", tag="hier") - """ - - from prody import writeDCD - - - distance_matrix = _validateDistanceMatrix(distance_matrix) - - cluster_ids = np.asarray(cluster_ids) - if cluster_ids.ndim != 1: - raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") - if cluster_ids.size == 0: - raise ValueError("cluster_ids is empty.") - if len(cluster_ids) != trajectory.numFrames(): - raise ValueError("cluster_ids must have one entry per trajectory frame.") - - if not isinstance(write_dcd, bool): - raise TypeError(f"write_dcd must be a bool, but got {type(write_dcd).__name__}") - - clusters = np.unique(cluster_ids) - clusters = clusters[clusters > 0] - num_clusters = clusters.size - if num_clusters == 0: - raise ValueError("No clusters were found. All frames are labeled as noise") - - exported_files = [] - - # NOTE: Loads the entire aligned trajectory into memory - _, aligned_coords = alignTrajectory(atoms, trajectory, align=align, select="all") - - for cluster in clusters: - cluster_indices = getCluster(cluster_ids, cluster) - medoid = getClusterMedoid(distance_matrix, cluster_indices = cluster_indices) - - cluster_coords = aligned_coords[cluster_indices] - - cluster_atoms = atoms.copy() - cluster_atoms.setCoords(cluster_coords[0]) - - if len(cluster_coords) > 1: - cluster_atoms.addCoordset(cluster_coords[1:]) - - if write_dcd: - dcd_filename = f"{system}_{tag}_n{num_clusters}_cluster{cluster}.dcd" - writeDCD(dcd_filename, cluster_atoms) - exported_files.append(dcd_filename) - - medoid_atoms = atoms.copy() - medoid_atoms.setCoords(aligned_coords[medoid["global"]]) - - pdb_filename = f"{system}_{tag}_n{num_clusters}_cluster{cluster}_medoid.pdb" - writePDB(pdb_filename, medoid_atoms) - exported_files.append(pdb_filename) - - return exported_files - - - - - class Interactions(object): """Class for Interaction analysis of proteins.""" diff --git a/prody/utilities/catchall.py b/prody/utilities/catchall.py index c07bc6d8c..1d15aa559 100644 --- a/prody/utilities/catchall.py +++ b/prody/utilities/catchall.py @@ -10,7 +10,7 @@ __all__ = ['calcTree', 'writeTree', 'parseTree', 'clusterMatrix', - 'showLines', 'showMatrix', 'showBars', + 'showLines', 'showMatrix', 'showBars', 'showHistogram', 'reorderMatrix', 'findSubgroups', 'getCoords', 'getLinkage', 'getTreeFromLinkage', 'clusterSubfamilies', 'calcRMSDclusters', 'calcGromosClusters', 'calcGromacsClusters', @@ -1283,3 +1283,109 @@ def calcKmedoidClusters(coordsets, nClusters): labels = c.labels_ _, counts = np.unique(labels, return_counts=True) return c.medoid_indices_, labels, counts + + +def showHistogram(data, *args, **kwargs): + """ + Plots the distribution of values on the current axis. + The input may be either a one-dimensional array or a symmetric square 2D matrix. + + + :arg data: 1D array or 2D square symmetric matrix. + If 2D, the upper triangle (excluding diagonal) is automatically extracted. + :type data: :class:`numpy.ndarray` + + :arg *args: positional arguments passed directly to Seaborn's ``histplot`` function. + :type *args: tuple + + :arg title: title of the plot. + Default is ``'Distribution'``. + :type title: str + + :arg xlabel: label for the x-axis. + Default is ``'Value'``. + :type xlabel: str + + :arg ylabel: label for the y-axis. + Default is ``'Frequency'``. + :type ylabel: str + + :arg grid: whether to display horizontal grid lines. + Default is ``True``. + :type grid: bool + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :arg **kwargs: keyword arguments passed directly to Seaborn's ``histplot`` function + :type **kwargs: dict + + :returns: the Matplotlib axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> distance_matrix = calcPairwiseRMSD(aligned_coords) + >>> plt.figure() + >>> showHistogram(distance_matrix, xlabel='RMSD [Å]') + >>> plt.show() + """ + + import matplotlib.pyplot as plt + try: + import seaborn as sns + except ImportError: + raise ImportError("The 'seaborn' package is required to display the histogram." + "\nPlease install it using 'pip install seaborn'." + "\nAlternatively, use standard matplotlib.pyplot.hist() for basic plots.") + + + data_array = np.asarray(data) + + if data_array.ndim == 2: + if data_array.shape[0] != data_array.shape[1]: + raise ValueError(f"2D data matrix must be square, but got shape {data_array.shape}") + values = data_array[np.triu_indices_from(data_array, k=1)] + elif data_array.ndim == 1: + values = data_array + else: + raise ValueError(f"Expected 1D or 2D array, but got shape {data_array.shape}") + + ax = kwargs.pop('ax', None) + if ax is None: + ax = plt.gca() + + title = kwargs.pop('title', 'Distribution') + xlabel = kwargs.pop('xlabel', 'Value') + ylabel = kwargs.pop('ylabel', 'Frequency') + grid = kwargs.pop('grid', True) + label = kwargs.get('label', None) + + # Seaborn Defaults + kwargs.setdefault('bins', 50) + kwargs.setdefault('element', 'bars') + kwargs.setdefault('stat', 'count') + kwargs.setdefault('kde', False) + kwargs.setdefault('alpha', 0.5) + kwargs.setdefault('color', 'teal') + kwargs.setdefault('edgecolor', 'black') + + if 'lw' in kwargs: + kwargs['linewidth'] = kwargs.pop('lw') + else: + kwargs.setdefault('linewidth', 0.8) + + sns.histplot(values, *args, ax=ax, **kwargs) + + ax.set_title(title) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + + if label is not None: + ax.legend() + + if grid: + ax.grid(axis='y', alpha=0.3) + + return ax \ No newline at end of file From e09413d610a5e9f5019cea3a5cf315b0af7925cf Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Tue, 11 Aug 2026 10:29:52 +0200 Subject: [PATCH 04/10] Fixed title bug in showPairwiseRMSDHeatmap --- prody/dynamics/plotting.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/prody/dynamics/plotting.py b/prody/dynamics/plotting.py index b52535410..6a8cfcc7e 100644 --- a/prody/dynamics/plotting.py +++ b/prody/dynamics/plotting.py @@ -2879,7 +2879,9 @@ def showPairwiseRMSDHeatmap(distance_matrix=None, aligned_coords=None, *args, ** distance_matrix = _validateDistanceMatrix(distance_matrix) + title = kwargs.pop('title', 'Pairwise RMSD Distance Matrix') ax = kwargs.pop('ax', None) + if ax is not None: plt.sca(ax) else: @@ -2889,7 +2891,6 @@ def showPairwiseRMSDHeatmap(distance_matrix=None, aligned_coords=None, *args, ** kwargs.setdefault('origin', 'upper') showMatrix(distance_matrix, *args, **kwargs) - title = kwargs.pop('title', 'Pairwise RMSD Distance Matrix') ax.set_title(title) return ax From b594caf69138539237fda0ac73c3071e1c84a84f Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Wed, 12 Aug 2026 11:34:42 +0200 Subject: [PATCH 05/10] Fixed int appearence on showClusterStatisticsTable - General minor fixes and changes --- prody/dynamics/__init__.py | 2 +- prody/dynamics/plotting.py | 4 +- prody/dynamics/rmsd_clustering.py | 87 +++++++++++++++++++++---------- prody/proteins/interactions.py | 1 - prody/utilities/catchall.py | 36 ++++++++++++- 5 files changed, 98 insertions(+), 32 deletions(-) diff --git a/prody/dynamics/__init__.py b/prody/dynamics/__init__.py index e0da0e780..c2009f308 100644 --- a/prody/dynamics/__init__.py +++ b/prody/dynamics/__init__.py @@ -395,4 +395,4 @@ from . import rmsd_clustering from .rmsd_clustering import * -__all__.extend(rmsd_clustering.__all__) \ No newline at end of file +__all__.extend(rmsd_clustering.__all__) diff --git a/prody/dynamics/plotting.py b/prody/dynamics/plotting.py index 6a8cfcc7e..11a60106b 100644 --- a/prody/dynamics/plotting.py +++ b/prody/dynamics/plotting.py @@ -2794,8 +2794,8 @@ def showRMSDEvolution(rmsd_array=None, ref_coords=None, aligned_coords=None, *ar raise ValueError(f"Incompatible shapes: reference is {ref_coords.shape}, but aligned frames have {aligned_coords.shape[1:]}.") rmsd_array = calcRMSD(ref_coords, target=aligned_coords) - - rmsd_array = np.asarray(rmsd_array) + else: + rmsd_array = np.asarray(rmsd_array) if rmsd_array.ndim != 1: raise ValueError(f"rmsd_array must be one-dimensional, but got shape {rmsd_array.shape}") diff --git a/prody/dynamics/rmsd_clustering.py b/prody/dynamics/rmsd_clustering.py index 192fcc738..f2af4ceab 100644 --- a/prody/dynamics/rmsd_clustering.py +++ b/prody/dynamics/rmsd_clustering.py @@ -475,7 +475,8 @@ def showClusterStatisticsTable(all_stats, dissimilarity="RMSD", units="Å", show :arg **kwargs: keyword arguments passed directly to the ``tabulate`` function :type **kwargs: dict - :returns: a dictionary containing the table, the row labels and the headers, with their names as keys + :returns: a dictionary containing the raw numeric table data, the row labels, + and the headers, with their names as keys :rtype: dict Example usage: @@ -484,13 +485,6 @@ def showClusterStatisticsTable(all_stats, dissimilarity="RMSD", units="Å", show >>> showClusterStatisticsTable(all_stats) """ - try: - from tabulate import tabulate - except ImportError: - raise ImportError("The 'tabulate' package is required to display the table. " - "Please install it using 'pip install tabulate'.") - - if isinstance(all_stats, dict): all_stats = [all_stats] elif not isinstance(all_stats, list): @@ -514,37 +508,63 @@ def showClusterStatisticsTable(all_stats, dissimilarity="RMSD", units="Å", show if missing: raise ValueError(f"Statistics dictionary is missing required keys: {sorted(missing)}.") + + floatfmt = kwargs.pop('floatfmt', '.2f') + + # The third element marks whether the metric is a floating point number (True) or an exact integer (False) metrics = [ - ("population", "Total Frames"), - ("pct", "Population Percentage (%)"), - ("medoid_global", "Medoid Frame (Global)"), - ("medoid_local", "Medoid Frame (Within Cluster)"), - ("mean", f"Mean {dissimilarity} [{units}]"), - ("std", f"Std {dissimilarity} [{units}]"), - ("median", f"Median {dissimilarity} [{units}]"), - ("iqr", f"IQR [{units}]"), - ("p95", f"95th Percentile [{units}]"), - ("max", f"Max {dissimilarity} [{units}]") + ("population", "Total Frames", False), + ("pct", "Population Percentage (%)", True), + ("medoid_global", "Medoid Frame (Global)", False), + ("medoid_local", "Medoid Frame (Within Cluster)", False), + ("mean", f"Mean {dissimilarity} [{units}]", True), + ("std", f"Std {dissimilarity} [{units}]", True), + ("median", f"Median {dissimilarity} [{units}]", True), + ("iqr", f"IQR [{units}]", True), + ("p95", f"95th Percentile [{units}]", True), + ("max", f"Max {dissimilarity} [{units}]", True) ] - table = [] + table_raw = [] + table_display = [] row_labels = [] - for key, label in metrics: + for key, label, is_float in metrics: row_labels.append(label) - table.append([stat[key] for stat in all_stats]) + raw_row = [stat[key] for stat in all_stats] + table_raw.append(raw_row) + + fmt = floatfmt if is_float else ".0f" + + formatted_row = [] + for val in raw_row: + if val is None: + formatted_row.append("N/A") + else: + try: + formatted_row.append(format(val, fmt)) + except ValueError: + formatted_row.append(str(val)) + + table_display.append(formatted_row) if not isinstance(show, bool): raise TypeError("show must be a bool.") if show: + try: + from tabulate import tabulate + except ImportError: + raise ImportError("The 'tabulate' package is required to display the table. " + "Please install it using 'pip install tabulate' or set show=False.") + kwargs.setdefault('tablefmt', 'fancy_grid') - kwargs.setdefault('floatfmt', '.4f') kwargs.setdefault('stralign', 'center') + kwargs.setdefault('numalign', 'center') - print(tabulate(table, headers=headers, showindex=row_labels, **kwargs)) + print(tabulate(table_display, headers=headers, showindex=row_labels, **kwargs)) - return {"table": table, "row_labels": row_labels, "headers": headers} + return {"table": table_raw, "row_labels": row_labels, "headers": headers} def writeClusters(atoms, trajectory, distance_matrix, cluster_ids, write_dcd=True, @@ -660,8 +680,8 @@ def clusterHierarchical(distance_matrix, method='average', cutoff=None): of the largest gap between consecutive linkage distances. Recommendation: Plot the dendrogram using showDendrogram() and choose the cutoff - manually whenever possible. - + manually whenever possible. + :arg distance_matrix: either a one-dimensional condensed distance matrix or a two-dimensional pairwise distance matrix. @@ -683,6 +703,12 @@ def clusterHierarchical(distance_matrix, method='average', cutoff=None): * the hierarchical clustering linkage matrix :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`) + + .. seealso:: + :func:`clusterMatrix` + Hierarchical clustering of a distance matrix using SciPy. + + Example usage: >>> distance_matrix = calcPairwiseRMSD(aligned_coords) >>> cluster_ids, linkage_matrix = clusterHierarchical(distance_matrix, cutoff='auto') @@ -923,6 +949,12 @@ def clusterKMedoids(distance_matrix, k, method='alternate', method_sklearn='pam' * the final sum of distances from each point to its nearest medoid :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`, float) + + .. seealso:: + :func:`prody.utilities.catchall.calcKmedoidClusters` + Direct K-Medoids clustering on coordinate sets using ``scikit-learn-extra``. + + Example usage: >>> distance_matrix = calcPairwiseRMSD(aligned_coords) >>> cluster_ids, medoids, _ = clusterKMedoids(distance_matrix, 4, method = 'alternate', @@ -1488,4 +1520,5 @@ def _orderOPTICSSklearn(distance_matrix, minPts): reachability = optics.reachability_ ordering = optics.ordering_ - return reachability, ordering \ No newline at end of file + return reachability, ordering + \ No newline at end of file diff --git a/prody/proteins/interactions.py b/prody/proteins/interactions.py index a4208f62f..3edba3cc6 100644 --- a/prody/proteins/interactions.py +++ b/prody/proteins/interactions.py @@ -4240,7 +4240,6 @@ def calcSignatureInteractions(PDB_folder, **kwargs): - class Interactions(object): """Class for Interaction analysis of proteins.""" diff --git a/prody/utilities/catchall.py b/prody/utilities/catchall.py index 1d15aa559..a8dd996a0 100644 --- a/prody/utilities/catchall.py +++ b/prody/utilities/catchall.py @@ -416,6 +416,7 @@ def clusterMatrix(distance_matrix=None, similarity_matrix=None, labels=None, ret return the sorted matrix, indices used for sorting, sorted labels (if **labels** are passed), and linkage matrix (if **return_linkage** is **True**). + :arg distance_matrix: an N-by-N matrix containing some measure of distance such as 1. - seqid_matrix (Hamming distance), rmsds, or distances in PCA space :type distance_matrix: :class:`~numpy.ndarray` @@ -438,6 +439,11 @@ def clusterMatrix(distance_matrix=None, similarity_matrix=None, labels=None, ret Other arguments for :func:`~scipy.hierarchy.linkage` and :func:`~scipy.hierarchy.dendrogram` can also be provided and will be taken as **kwargs**. + + + .. seealso:: + :func:`clusterHierarchical` + Hierarchical clustering of a distance matrix using SciPy. """ import scipy.cluster.hierarchy as sch @@ -1273,6 +1279,33 @@ def calcRMSDclusters(rmsd_matrix, c, labels=None): calcGromacsClusters = calcRMSDclusters def calcKmedoidClusters(coordsets, nClusters): + """ + Performs K-Medoids clustering directly on coordinate sets using ``scikit-learn-extra``. + + Flattens the input coordinate sets into a 2D matrix and fits a K-Medoids model. + Note that this function requires the ``scikit-learn-extra`` package to be installed. + + The returned cluster labels are 0-indexed. + + + :arg coordsets: coordinate sets of shape ``(n_conformations, n_atoms, 3)`` or any array + where the first dimension represents individual frames/conformations. + :type coordsets: :class:`numpy.ndarray` + + :arg nClusters: prespecified number of clusters to form. + :type nClusters: int + + :returns: a tuple of: + * a one-dimensional array of shape ``(nClusters,)`` containing the indices of the medoid frames + * a one-dimensional array containing the 0-indexed cluster label for each frame + * a one-dimensional array containing the number of items in each cluster + :rtype: tuple(:class:`numpy.ndarray`, :class:`numpy.ndarray`, :class:`numpy.ndarray`) + + + .. seealso:: + :func:`clusterKMedoids` + K-Medoids clustering on distance matrices supporting multiple backends and 1-based indexing. + """ try: from sklearn_extra.cluster import KMedoids except ImportError: @@ -1388,4 +1421,5 @@ def showHistogram(data, *args, **kwargs): if grid: ax.grid(axis='y', alpha=0.3) - return ax \ No newline at end of file + return ax + \ No newline at end of file From cb41a01a935ce7fa797ba037e1dcc4ae3ee722fe Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Wed, 12 Aug 2026 11:45:07 +0200 Subject: [PATCH 06/10] Revert changes to interactions.py --- prody/proteins/interactions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/prody/proteins/interactions.py b/prody/proteins/interactions.py index 3edba3cc6..dfca7a1d3 100644 --- a/prody/proteins/interactions.py +++ b/prody/proteins/interactions.py @@ -4237,7 +4237,7 @@ def calcSignatureInteractions(PDB_folder, **kwargs): # Proceed with plotting plot_barh(result, bond_type, n_per_plot=n_per_plot, min_height=min_height) - + class Interactions(object): From 488be8931ebd0dfd9a62bf3a2517b24563463c0b Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Wed, 12 Aug 2026 12:39:00 +0200 Subject: [PATCH 07/10] Matplotlib fallback for Seaborn histograms --- prody/dynamics/plotting.py | 59 +++++++++++++++++----------- prody/utilities/catchall.py | 78 ++++++++++++++++++++----------------- 2 files changed, 80 insertions(+), 57 deletions(-) diff --git a/prody/dynamics/plotting.py b/prody/dynamics/plotting.py index 11a60106b..56859ea83 100644 --- a/prody/dynamics/plotting.py +++ b/prody/dynamics/plotting.py @@ -2894,8 +2894,8 @@ def showPairwiseRMSDHeatmap(distance_matrix=None, aligned_coords=None, *args, ** ax.set_title(title) return ax - - + + def showClusterRMSDComparison(all_stats, *args, **kwargs): """ Overlays the internal RMSD distributions of all clusters onto a single plot @@ -2903,23 +2903,27 @@ def showClusterRMSDComparison(all_stats, *args, **kwargs): Expects the list of cluster statistics dictionaries generated by :func:`calcAllClusterStatistics`. + Uses Seaborn's :func:`sns.histplot` if available, falling back to + Matplotlib's :meth:`matplotlib.axes.Axes.hist` if Seaborn is not installed. + :arg all_stats: list of cluster statistics dictionaries. Each dictionary must contain the keys ``"cluster"`` and ``"distances"``. :type all_stats: list of dict, or dict - :arg *args: positional arguments passed directly to Seaborn's ``histplot`` function. + :arg *args: positional arguments passed directly to Seaborn's ``histplot`` + or Matplotlib's ``hist`` function. :type *args: tuple :arg title: the title of the generated plot. Default is ``'Intra-Cluster RMSD Distributions'`` :type title: str - :arg xlabel: the label for the x-axis. + :arg xlabel: the label for the X axis. Default is ``'RMSD to Medoid [Å]'`` :type xlabel: str - :arg ylabel: the label for the y-axis. + :arg ylabel: the label for the Y axis. Default is ``'Frequency'`` :type ylabel: str @@ -2935,8 +2939,8 @@ def showClusterRMSDComparison(all_stats, *args, **kwargs): Default is ``None`` and the current axes are used. :type ax: :class:`matplotlib.axes.Axes` - :arg **kwargs: keyword arguments passed directly to Seaborn's ``histplot`` function - (e.g., ``bins``, ``kde``, ``element``, ``stat``, ``alpha``, ``lw``). + :arg **kwargs: keyword arguments passed directly to Seaborn's ``histplot`` + or Matplotlib's ``hist`` function. :type **kwargs: dict :returns: the Matplotlib axes containing the plot. @@ -2953,12 +2957,10 @@ def showClusterRMSDComparison(all_stats, *args, **kwargs): import matplotlib.pyplot as plt try: import seaborn as sns + has_seaborn = True except ImportError: - raise ImportError("The 'seaborn' package is required to display the histogram." - "\nPlease install it using 'pip install seaborn'." - "\nAlternatively, use standard matplotlib.pyplot.hist() for basic plots.") - - + has_seaborn = False + if isinstance(all_stats, dict): all_stats = [all_stats] @@ -2976,18 +2978,28 @@ def showClusterRMSDComparison(all_stats, *args, **kwargs): user_label = kwargs.pop('label', None) - # Seaborn Defaults - kwargs.setdefault('bins', 50) - kwargs.setdefault('element', 'step') - kwargs.setdefault('stat', 'count') - kwargs.setdefault('kde', True) - kwargs.setdefault('alpha', 0.5) - if 'lw' in kwargs: kwargs['linewidth'] = kwargs.pop('lw') else: kwargs.setdefault('linewidth', 1.5) + kwargs.setdefault('bins', 50) + kwargs.setdefault('alpha', 0.5) + + if has_seaborn: + kwargs.setdefault('element', 'step') + kwargs.setdefault('stat', 'count') + kwargs.setdefault('kde', True) + else: # matplotlib fallback + LOGGER.info("Package 'seaborn' not found; falling back to matplotlib.pyplot.hist().") + kde = kwargs.pop('kde', None) + kwargs.pop('element', None) + kwargs.pop('stat', None) + if kde: + LOGGER.warning("Kernel density estimation (kde=True) requires 'seaborn' and will be ignored.") + + kwargs.setdefault('histtype', 'step') + for cluster_stats in all_stats: if "cluster" not in cluster_stats or "distances" not in cluster_stats: raise ValueError("Each cluster statistics dictionary must contain 'cluster' and 'distances'.") @@ -2996,8 +3008,11 @@ def showClusterRMSDComparison(all_stats, *args, **kwargs): if user_label: cluster_label = f"{user_label} - {cluster_label}" - sns.histplot(cluster_stats["distances"], *args, label=cluster_label, ax=ax, **kwargs) - + if has_seaborn: + sns.histplot(cluster_stats["distances"], *args, label=cluster_label, ax=ax, **kwargs) + else: # matplotlib fallback + ax.hist(cluster_stats["distances"], *args, label=cluster_label, **kwargs) + ax.set_title(title) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) @@ -3007,4 +3022,4 @@ def showClusterRMSDComparison(all_stats, *args, **kwargs): if grid: ax.grid(axis='y', alpha=0.3) - return ax \ No newline at end of file + return ax diff --git a/prody/utilities/catchall.py b/prody/utilities/catchall.py index a8dd996a0..94fde93a7 100644 --- a/prody/utilities/catchall.py +++ b/prody/utilities/catchall.py @@ -1323,35 +1323,34 @@ def showHistogram(data, *args, **kwargs): Plots the distribution of values on the current axis. The input may be either a one-dimensional array or a symmetric square 2D matrix. + Uses Seaborn's :func:`sns.histplot` if available, falling back to + Matplotlib's :meth:`matplotlib.axes.Axes.hist` if Seaborn is not installed. + :arg data: 1D array or 2D square symmetric matrix. If 2D, the upper triangle (excluding diagonal) is automatically extracted. :type data: :class:`numpy.ndarray` - :arg *args: positional arguments passed directly to Seaborn's ``histplot`` function. + :arg *args: positional arguments passed directly to Seaborn's ``histplot`` + or Matplotlib's ``hist`` function. :type *args: tuple - :arg title: title of the plot. - Default is ``'Distribution'``. + :arg title: title of the plot. Default is ``'Distribution'``. :type title: str - :arg xlabel: label for the x-axis. - Default is ``'Value'``. + :arg xlabel: label for the x-axis. Default is ``'Value'``. :type xlabel: str - :arg ylabel: label for the y-axis. - Default is ``'Frequency'``. + :arg ylabel: label for the y-axis. Default is ``'Frequency'``. :type ylabel: str - :arg grid: whether to display horizontal grid lines. - Default is ``True``. + :arg grid: whether to display horizontal grid lines. Default is ``True``. :type grid: bool - :arg ax: axes on which to draw the plot. - Default is ``None`` and the current axes are used. + :arg ax: axes on which to draw the plot. Default is ``None`` (current axes used). :type ax: :class:`matplotlib.axes.Axes` - :arg **kwargs: keyword arguments passed directly to Seaborn's ``histplot`` function + :arg **kwargs: keyword arguments passed to Seaborn (if installed) or Matplotlib. :type **kwargs: dict :returns: the Matplotlib axes containing the plot. @@ -1362,20 +1361,18 @@ def showHistogram(data, *args, **kwargs): >>> distance_matrix = calcPairwiseRMSD(aligned_coords) >>> plt.figure() >>> showHistogram(distance_matrix, xlabel='RMSD [Å]') - >>> plt.show() + >>> plt.show() """ - + import matplotlib.pyplot as plt try: import seaborn as sns + has_seaborn = True except ImportError: - raise ImportError("The 'seaborn' package is required to display the histogram." - "\nPlease install it using 'pip install seaborn'." - "\nAlternatively, use standard matplotlib.pyplot.hist() for basic plots.") - - + has_seaborn = False + data_array = np.asarray(data) - + if data_array.ndim == 2: if data_array.shape[0] != data_array.shape[1]: raise ValueError(f"2D data matrix must be square, but got shape {data_array.shape}") @@ -1384,40 +1381,51 @@ def showHistogram(data, *args, **kwargs): values = data_array else: raise ValueError(f"Expected 1D or 2D array, but got shape {data_array.shape}") - + ax = kwargs.pop('ax', None) if ax is None: ax = plt.gca() - + title = kwargs.pop('title', 'Distribution') xlabel = kwargs.pop('xlabel', 'Value') ylabel = kwargs.pop('ylabel', 'Frequency') grid = kwargs.pop('grid', True) label = kwargs.get('label', None) - # Seaborn Defaults - kwargs.setdefault('bins', 50) - kwargs.setdefault('element', 'bars') - kwargs.setdefault('stat', 'count') - kwargs.setdefault('kde', False) - kwargs.setdefault('alpha', 0.5) - kwargs.setdefault('color', 'teal') - kwargs.setdefault('edgecolor', 'black') - if 'lw' in kwargs: kwargs['linewidth'] = kwargs.pop('lw') else: kwargs.setdefault('linewidth', 0.8) - sns.histplot(values, *args, ax=ax, **kwargs) - + kwargs.setdefault('bins', 50) + kwargs.setdefault('color', 'teal') + kwargs.setdefault('alpha', 0.5) + kwargs.setdefault('edgecolor', 'black') + + if has_seaborn: + kwargs.setdefault('element', 'bars') + kwargs.setdefault('stat', 'count') + kwargs.setdefault('kde', False) + + sns.histplot(values, *args, ax=ax, **kwargs) + else: # matplotlib fallback + LOGGER.info("Package 'seaborn' not found; falling back to matplotlib.pyplot.hist().") + kde = kwargs.pop('kde', None) + kwargs.pop('element', None) + kwargs.pop('stat', None) + + if kde: + LOGGER.warning("Kernel density estimation (kde=True) requires 'seaborn' and will be ignored.") + + ax.hist(values, *args, **kwargs) + ax.set_title(title) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) - + if label is not None: ax.legend() - + if grid: ax.grid(axis='y', alpha=0.3) From 474e0d1947d7c91c12f0b253f27d46ff4f092dad Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Wed, 12 Aug 2026 16:22:22 +0200 Subject: [PATCH 08/10] Fixed colormap bug on showReachabilityPlot --- prody/dynamics/rmsd_clustering.py | 36 +++++++++++++++++++++++-------- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/prody/dynamics/rmsd_clustering.py b/prody/dynamics/rmsd_clustering.py index f2af4ceab..488724bb4 100644 --- a/prody/dynamics/rmsd_clustering.py +++ b/prody/dynamics/rmsd_clustering.py @@ -1330,7 +1330,8 @@ def _validateMinPts(minPts, distance_matrix): raise ValueError(f"minPts must be between 1 and {total_points}.") -def showReachabilityPlot(distance_matrix, *args, minPts=None, method='custom', eps=None, fill=True, **kwargs): +def showReachabilityPlot(distance_matrix, *args, minPts=None, method='custom', eps=None, + fill=True, colors=None, **kwargs): """ Plots the reachability plot using a simplified OPTICS algorithm. The reachability plot should be used to determine the most suitable eps @@ -1366,6 +1367,10 @@ def showReachabilityPlot(distance_matrix, *args, minPts=None, method='custom', e Default is ``True``. :type fill: bool + :arg colors: custom color mapping for the valleys (list, dict, or colormap name). + Default is ``None`` and uses Matplotlib's tab10 palette. + :type colors: list, dict, str, or None + :arg **kwargs: additional keyword arguments passed to Matplotlib's ``plot`` function :type **kwargs: dict @@ -1375,8 +1380,8 @@ def showReachabilityPlot(distance_matrix, *args, minPts=None, method='custom', e Example usage: >>> import matplotlib.pyplot as plt >>> plt.figure() - >>> showReachabilityPlot(distance_matrix, minPts=20, eps=1.8) - >>> plt.show() + >>> showReachabilityPlot(distance_matrix, minPts = 20, eps = 2.8) + >>> plt.show() """ import matplotlib.pyplot as plt @@ -1434,7 +1439,7 @@ def showReachabilityPlot(distance_matrix, *args, minPts=None, method='custom', e below_eps_mask = (y <= eps) segments = [] start = None - + for i, below in enumerate(below_eps_mask): if below and start is None: start = i @@ -1445,9 +1450,22 @@ def showReachabilityPlot(distance_matrix, *args, minPts=None, method='custom', e if start is not None: segments.append((start, len(below_eps_mask))) - colors = plt.cm.tab10(np.linspace(0, 1, len(segments))) - - for (start, end), valley_color in zip(segments, colors): + if colors is None: + cmap = plt.get_cmap('tab10') + elif isinstance(colors, str): + cmap = plt.get_cmap(colors) + else: + cmap = colors + + for i, (start, end) in enumerate(segments): + if isinstance(cmap, dict): + valley_color = cmap.get(i, 'gray') + elif isinstance(cmap, (list, tuple)): + valley_color = cmap[i % len(cmap)] + else: + num_colors = getattr(cmap, 'N', 10) + valley_color = cmap(i % num_colors) + ax.fill_between(x[start:end], y[start:end], eps, color=valley_color) if 'label' in kwargs: @@ -1455,8 +1473,8 @@ def showReachabilityPlot(distance_matrix, *args, minPts=None, method='custom', e ax.grid(axis='y', linestyle='--', linewidth=0.8, alpha=0.3) - return ax - + return ax + def _orderOPTICSCustom(distance_matrix, minPts): """OPTICS algorithm with built-in modules and NumPy""" From 19d8ed8c6711950e741759392b4f0da4e1569bc5 Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Wed, 12 Aug 2026 16:55:01 +0200 Subject: [PATCH 09/10] Added warning for cluster ID 0 handling --- prody/dynamics/rmsd_clustering.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/prody/dynamics/rmsd_clustering.py b/prody/dynamics/rmsd_clustering.py index 488724bb4..413ecbbb2 100644 --- a/prody/dynamics/rmsd_clustering.py +++ b/prody/dynamics/rmsd_clustering.py @@ -367,11 +367,11 @@ def calcAllClusterStatistics(distance_matrix, cluster_ids): Returns a list of dictionaries, each corresponding to a cluster with the statistic quantities as keys. - :arg distance_matrix: a two-dimensional square, symmetric pairwise distance matrix. :type distance_matrix: :class:`numpy.ndarray` - :arg cluster_ids: a one-dimensional array of the cluster IDs per element + :arg cluster_ids: a one-dimensional array of 1-indexed cluster IDs per element. + Values less than 1 (e.g. 0, -1) are omitted. :type cluster_ids: :class:`numpy.ndarray` :returns: a list of dictionaries, one dictionary for each cluster with the statistic @@ -379,10 +379,13 @@ def calcAllClusterStatistics(distance_matrix, cluster_ids): :rtype: list of dict Example usage: - >>> cluster_ids, _= clusterHierarchical(distance_matrix, 4) - >>> all_stats = calcAllClusterStatistics(distance_matrix, cluster_ids) + >>> cluster_ids, _ = prody.clusterHierarchical(distance_matrix, 4) + >>> all_stats = prody.calcAllClusterStatistics(distance_matrix, cluster_ids) """ + import warnings + + cluster_ids = np.asarray(cluster_ids) if cluster_ids.ndim != 1: raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") @@ -394,13 +397,19 @@ def calcAllClusterStatistics(distance_matrix, cluster_ids): f"does not match distance_matrix frames ({distance_matrix.shape[0]}).") clusters = np.unique(cluster_ids) + + if 0 in clusters: + warnings.warn("Cluster ID '0' detected in cluster_ids." + " This library uses 1-indexed clusters; ID 0 will be omitted from statistics.", + UserWarning) + all_stats = [] for c in clusters: if c <= 0: continue - cluster_stats = calcClusterStatistics(distance_matrix, cluster_ids = cluster_ids, cluster_number = int(c)) + cluster_stats = calcClusterStatistics(distance_matrix, cluster_ids=cluster_ids, cluster_number=int(c)) all_stats.append(cluster_stats) return all_stats From 1cd37b9e2e7bb41d70ef25d54441c7e3f55e92f4 Mon Sep 17 00:00:00 2001 From: kontheodosiadis Date: Thu, 13 Aug 2026 11:28:40 +0200 Subject: [PATCH 10/10] Added two new functions for tracking the clusters in the trajectory --- prody/dynamics/plotting.py | 327 ++++++++++++++++++++++++++++++++++++- 1 file changed, 325 insertions(+), 2 deletions(-) diff --git a/prody/dynamics/plotting.py b/prody/dynamics/plotting.py index 56859ea83..780fc6af0 100644 --- a/prody/dynamics/plotting.py +++ b/prody/dynamics/plotting.py @@ -41,7 +41,8 @@ 'showPerturbResponse', 'showTree', 'showTree_networkx', 'showAtomicMatrix', 'pimshow', 'showAtomicLines', 'pplot', 'showDomainBar', 'showAtomicBars', 'showSelectionMatrix', - 'showRMSDEvolution', 'showPairwiseRMSDHeatmap', 'showClusterRMSDComparison'] + 'showRMSDEvolution', 'showPairwiseRMSDHeatmap', 'showClusterRMSDComparison', + 'showClusterRMSDEvolution', 'showClusterBarcode'] def showEllipsoid(modes, onto=None, n_std=2, scale=1., *args, **kwargs): @@ -2998,7 +2999,7 @@ def showClusterRMSDComparison(all_stats, *args, **kwargs): if kde: LOGGER.warning("Kernel density estimation (kde=True) requires 'seaborn' and will be ignored.") - kwargs.setdefault('histtype', 'step') + kwargs.setdefault('histtype', 'stepfilled') for cluster_stats in all_stats: if "cluster" not in cluster_stats or "distances" not in cluster_stats: @@ -3023,3 +3024,325 @@ def showClusterRMSDComparison(all_stats, *args, **kwargs): ax.grid(axis='y', alpha=0.3) return ax + + +def showClusterRMSDEvolution(cluster_ids, rmsd_array=None, ref_coords=None, aligned_coords=None, + colors=None, mode='background', alpha=0.3, *args, **kwargs): + """ + Plots the RMSD evolution over trajectory frames, highlighting or color-coding + the clusters for visual validation. + + Accepts either ``rmsd_array`` or ``ref_coords`` and ``aligned_coords`` for + RMSD calculation. + + + :arg cluster_ids: 1D array mapping each frame to its assigned cluster ID. + Cluster IDs are expected to be 1-based. Values less than 1 + (e.g., 0 or -1) are ignored when coloring clusters. + :type cluster_ids: :class:`numpy.ndarray` + + :arg rmsd_array: 1D array containing pre-calculated RMSD values per frame. + If ``None``, ``ref_coords`` and ``aligned_coords`` are used + to calculate the RMSD values. + :type rmsd_array: :class:`numpy.ndarray` + + :arg ref_coords: reference coordinates for RMSD calculation. + :type ref_coords: :class:`numpy.ndarray` + + :arg aligned_coords: aligned trajectory coordinates for RMSD calculation. + :type aligned_coords: :class:`numpy.ndarray` + + :arg colors: custom color mapping for clusters. Can be a list or tuple of + colors, a dictionary mapping cluster IDs to colors, or a + Matplotlib colormap name. If ``None``, the ``tab10`` colormap + is used. + :type colors: list, tuple, dict, str, or None + + :arg mode: how to visualize cluster membership. + Options are ``'background'`` and ``'dots'``. + Default is ``'background'``. + :type mode: str + + :arg alpha: transparency for background spans or scatter dots. + Default is ``0.3``. + :type alpha: float + + :arg *args: positional arguments passed directly to :func:`showRMSDEvolution`. + :type *args: tuple + + :arg **kwargs: additional keyword arguments passed directly to :func:`showRMSDEvolution`. + :type **kwargs: dict + + :returns: Matplotlib axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> ref_coords, aligned_coords = alignTrajectory(pdb, dcd, select='resname DAP') + >>> rmsd_evolution = calcRMSD(ref_coords, target=aligned_coords) + >>> rmsd_distance_matrix = calcPairwiseRMSD(aligned_coords) + >>> hier_cluster_ids, _ = clusterHierarchical(rmsd_distance_matrix) + >>> plt.figure() + >>> showClusterRMSDEvolution(cluster_ids=hier_cluster_ids, rmsd_array=rmsd_evolution, mode='background') + >>> plt.show() + """ + + import matplotlib.pyplot as plt + import warnings + + + cluster_ids = np.asarray(cluster_ids) + + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") + + if cluster_ids.size == 0: + raise ValueError("cluster_ids must not be empty.") + + if mode not in ('background', 'dots'): + raise ValueError(f"Unknown mode '{mode}'. Choose either 'background' or 'dots'.") + + n_frames = len(cluster_ids) + + if rmsd_array is None: + if ref_coords is None or aligned_coords is None: + raise ValueError("Either rmsd_array or both ref_coords and aligned_coords must be provided.") + + rmsd_array = calcRMSD(ref_coords, target=aligned_coords) + else: + rmsd_array = np.asarray(rmsd_array) + + if rmsd_array.ndim != 1: + raise ValueError(f"rmsd_array must be a 1D array, but got shape {rmsd_array.shape}.") + + if len(rmsd_array) != n_frames: + raise ValueError(f"rmsd_array must contain one value per frame, but got {len(rmsd_array)} RMSD values for {n_frames} frames.") + + ax = kwargs.pop('ax', None) + if ax is None: + ax = plt.gca() + + kwargs['ax'] = ax + + ax = showRMSDEvolution(rmsd_array=rmsd_array, *args, **kwargs) + unique_clusters = np.unique(cluster_ids) + + if 0 in unique_clusters: + warnings.warn("Cluster ID '0' detected in cluster_ids. " + "This library uses 1-indexed clusters; ID 0 will not be colored.", + UserWarning) + valid_clusters = [cluster for cluster in unique_clusters if cluster > 0] + + if colors is None: + cmap = plt.get_cmap('tab10') + elif isinstance(colors, str): + cmap = plt.get_cmap(colors) + else: + cmap = colors + + if isinstance(colors, (list, tuple)) and len(colors) == 0: + raise ValueError("colors cannot be an empty list or tuple.") + + num_colors = getattr(cmap, 'N', 10) + + cluster_color_map = {} + for i, cluster_number in enumerate(valid_clusters): + + if isinstance(colors, dict): + if cluster_number in colors: + cluster_color_map[cluster_number] = colors[cluster_number] + else: + cluster_color_map[cluster_number] = cmap(i % num_colors) + + elif isinstance(colors, (list, tuple)): + cluster_color_map[cluster_number] = colors[i % len(colors)] + + else: + cluster_color_map[cluster_number] = cmap(i % num_colors) + + if mode == 'background': + change_points = np.where(cluster_ids[:-1] != cluster_ids[1:])[0] + 1 + segment_starts = np.insert(change_points, 0, 0) + segment_ends = np.append(change_points, n_frames) + drawn_labels = set() + + for start, end in zip(segment_starts, segment_ends): + cluster_number = cluster_ids[start] + if cluster_number not in cluster_color_map: + continue + cluster_color = cluster_color_map[cluster_number] + + if cluster_number not in drawn_labels: + cluster_label = f"Cluster {cluster_number}" + drawn_labels.add(cluster_number) + else: + cluster_label = None + + ax.axvspan(start - 0.5, end - 0.5, color=cluster_color, alpha=alpha, label=cluster_label, linewidth=0) + + elif mode == 'dots': + for cluster_number, cluster_color in cluster_color_map.items(): + mask = cluster_ids == cluster_number + frame_indices = np.where(mask)[0] + ax.scatter(frame_indices, rmsd_array[mask], color=cluster_color, alpha=alpha, + s=15, label=f"Cluster {cluster_number}", zorder=3) + + if cluster_color_map: + handles, labels = ax.get_legend_handles_labels() + if handles: + def legend_sort_key(pair): + label = pair[1] + if label.startswith("Cluster "): + try: + return (1, int(label.split()[1])) + except (IndexError, ValueError): + pass + return (0, label) + + sorted_pairs = sorted(zip(handles, labels), key=legend_sort_key) + sorted_handles, sorted_labels = zip(*sorted_pairs) + + ax.legend(sorted_handles, sorted_labels, bbox_to_anchor=(1.01, 1), loc='upper left') + + return ax + + +def showClusterBarcode(cluster_ids, colors=None, *args, **kwargs): + """ + Creates a 1D barcode-style timeline of trajectory frames, where each frame + is represented by a vertical line colored according to its cluster ID. + + + :arg cluster_ids: 1D array assigning each frame to a cluster. + Cluster IDs are expected to be 1-based. Values less than 1 + (e.g., 0 or -1) are ignored when coloring clusters. + :type cluster_ids: :class:`numpy.ndarray` + + :arg colors: custom color mapping for clusters. Can be a list or tuple of + colors, a dictionary mapping cluster IDs to colors, or a + Matplotlib colormap name. If ``None``, the ``tab10`` colormap + is used. + :type colors: list, tuple, dict, str, or None + + :arg ax: axes on which to draw the plot. + Default is ``None`` and the current axes are used. + :type ax: :class:`matplotlib.axes.Axes` + + :arg title: title of the plot. + Default is ``'Trajectory Cluster Timeline'``. + :type title: str + + :arg xlabel: label for the x-axis. + Default is ``'Frame Index'``. + :type xlabel: str + + :arg *args: positional arguments passed directly to Matplotlib's ``axvspan`` function. + :type *args: tuple + + :arg **kwargs: keyword arguments passed directly to Matplotlib's ``axvspan`` function. + :type **kwargs: dict + + :returns: Matplotlib axes containing the plot. + :rtype: :class:`matplotlib.axes.Axes` + + Example usage: + >>> import matplotlib.pyplot as plt + >>> hier_cluster_ids, _ = clusterHierarchical(distance_matrix) + >>> plt.figure() + >>> showClusterBarcode(hier_cluster_ids) + >>> plt.show() + """ + + import matplotlib.pyplot as plt + import warnings + + + cluster_ids = np.asarray(cluster_ids) + + if cluster_ids.ndim != 1: + raise ValueError(f"cluster_ids must be a 1D array, but got shape {cluster_ids.shape}.") + if cluster_ids.size == 0: + raise ValueError("cluster_ids must not be empty.") + + ax = kwargs.pop('ax', None) + if ax is None: + ax = plt.gca() + + title = kwargs.pop('title', 'Trajectory Cluster Timeline') + xlabel = kwargs.pop('xlabel', 'Frame Index') + + n_frames = len(cluster_ids) + + unique_clusters = np.unique(cluster_ids) + if 0 in unique_clusters: + warnings.warn("Cluster ID '0' detected. " + "This library uses 1-indexed clusters; ID 0 will not be colored.", + UserWarning) + + valid_clusters = [cluster_number for cluster_number in unique_clusters if cluster_number > 0] + + if colors is None: + cmap = plt.get_cmap('tab10') + elif isinstance(colors, str): + cmap = plt.get_cmap(colors) + else: + cmap = colors + + if isinstance(colors, (list, tuple)) and len(colors) == 0: + raise ValueError("colors cannot be an empty list or tuple.") + + num_colors = getattr(cmap, 'N', 10) + + cluster_color_map = {} + + for i, cluster_number in enumerate(valid_clusters): + if isinstance(colors, dict): + cluster_color_map[cluster_number] = colors.get(cluster_number, cmap(i % num_colors)) + elif isinstance(colors, (list, tuple)): + cluster_color_map[cluster_number] = colors[i % len(colors)] + else: + cluster_color_map[cluster_number] = cmap(i % num_colors) + + kwargs.setdefault('linewidth', 0) + + change_points = np.where(cluster_ids[:-1] != cluster_ids[1:])[0] + 1 + segment_starts = np.insert(change_points, 0, 0) + segment_ends = np.append(change_points, n_frames) + + drawn_labels = set() + + for start, end in zip(segment_starts, segment_ends): + cluster_number = cluster_ids[start] + + if cluster_number not in cluster_color_map: + continue + + if cluster_number not in drawn_labels: + cluster_label = f"Cluster {cluster_number}" + drawn_labels.add(cluster_number) + else: + cluster_label = None + + ax.axvspan(start - 0.5, end - 0.5, *args, color=cluster_color_map[cluster_number], + label=cluster_label, **kwargs) + + ax.set_xlim(-0.5, n_frames - 0.5) + ax.set_ylim(0, 1) + + ax.set_yticks([]) + ax.spines['left'].set_visible(False) + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) + + ax.set_xlabel(xlabel) + ax.set_title(title) + + if cluster_color_map: + handles, labels = ax.get_legend_handles_labels() + if handles: + sorted_pairs = sorted(zip(handles, labels), key=lambda x: int(x[1].split()[1])) + sorted_handles, sorted_labels = zip(*sorted_pairs) + + ax.legend(sorted_handles, sorted_labels, bbox_to_anchor=(1.01, 1), loc='upper left') + + return ax