Trajectory Clustering, Statistical Analysis, and Visualization Tools by kontheodosiadis · Pull Request #2264 · prody/ProDy

kontheodosiadis · 2026-08-07T11:08:03Z

The code was developed under the supervision of Prof. Karolina Mikulska-Ruminska (NCU).

Summary

This pull request introduces a set of tools for clustering and analyzing Molecular Dynamics trajectories. It provides functions for trajectory alignment, RMSD calculations, distance matrix generation, clustering, cluster statistical analysis, visualization, and export of representative structures and clustered trajectories.

Workflow

The functions follow a logical pipline:

Frame-to-Reference Analysis

Align trajectory frames to a reference structure.
Calculate RMSD for arbitrary atom selections (e.g. ligands, proteins).

Frame-to-Frame Analysis

Generate pairwise RMSD distance matrices across trajectory frames.

Clustering

Apply clustering algorithms on the pairwise distance matrix.

Analysis and Visualization

Calculate descriptive cluster statistics.
Identify representative medoid structures.
Generate visualizations for cluster interpretation.

Export

Export representative medoid structures as .pdb files.
Export aligned cluster trajectories as .dcd files.

Key Features

Clustering Algorithms

Bottom-up Hierarchical Clustering
K-Medoids Clustering
DBSCAN Clustering and OPTICS Reachability Plot

Clustering algorithms support multiple backends, where needed. Optimized implementations from scikit-learn (or SciPy, where appropriate) can be used, while dependency-free implementations based on NumPy and the Python standard library are provided whenever practical.

Analysis Utilities

Vectorized pairwise RMSD distance matrix generation
Medoid identification
Cluster summary and statistics tables
Cluster trajectory export (PDB/DCD)

Visualization

Hierarchical clustering dendrograms with optional cutoff visualization
OPTICS reachability plots with optional DBSCAN $\epsilon$ threshold visualization
RMSD histograms with optional Kernel Density Estimation (KDE) overlays

@karolamik13 @jamesmkrieger

01_rmsd_frame_to_reference.ipynb
02_rmsd_frame_to_frame.ipynb
03_hierarchical_clustering.ipynb
04_kmedoids_clustering.ipynb
05_density_clustering.ipynb

jamesmkrieger

I still need to look at the details but in general this seems like a good addition. I have two high-level comments for now.

move this code into different files instead of interactions.py as it is based on rmsd, not interactions. It may well get split across multiple files, but you will probably have one main one like prody/dynamics/rmsd_ clustering.py

There are currently rmsd clustering and low level plotting tools in prody/utilities/catchall.py, so that could be a possible location for some or all of it, but probably it would be better to have a dedicated file prody/dynamics/rmsd_ clustering.py that could take most of the new code.

Plotting tools could go in either prody/dynamics/plotting.py or prody/utilities/catchall.py depending whether they are high-level (plotting particular dynamics-related features and including the calculations and plotting together) or low-level (just providing plotting utilities like histograms, which could also be quite useful).

prody/utilities/catchall.py has clusterSubfamilies and clusterMatrix that are data agnostic and take a matrix of similarities or distances/rmsds; calcRMSDclusters, which is aliased as calcGromosClusters and calcGromacsClusters, that is based on the method used on Gromacs; and calcKmedoidClusters.

There is some reinventing of the wheel here that needs to be consolidated with existing functions. I have commented on some examples.

jamesmkrieger · 2026-08-09T06:27:54Z

+    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


prody already has rmsd calculation. Is this version quicker? If so, it should be integrated with what we already have. If not, it should be replaced here with calcRMSD or getRMSD

James, should we still avoid f-strings (.format() was usually used to support Python 2.7 users)?
If so, we have some in the code (including the one above).

Probably doesn't matter that much at this point as there won't really be any Python 2.7 users besides the website and we're not rebuilding it at the moment

…undant Code

kontheodosiadis · 2026-08-10T13:23:51Z

I moved the files to other folders according to your suggestion. I also cleaned things up by removing unnecessary functions and utilizing existing ProDy ones where appropriate.

For now, I have kept my clustering functions separately in rmsd_clustering.py. While there is some overlap with existing functions, my implementation utilizes different libraries and supports multiple backends, so I think for the time-being it is safer to have them apart.

I have attached the notebooks with the corrected functions.
01_rmsd_frame_to_reference.ipynb
02_rmsd_frame_to_frame.ipynb
03_hierarchical_clustering.ipynb
04_kmedoids_clustering.ipynb
05_density_clustering.ipynb

jamesmkrieger · 2026-08-11T14:44:29Z

I moved the files to other folders according to your suggestion. I also cleaned things up by removing unnecessary functions and utilizing existing ProDy ones where appropriate.

Thanks. This looks good now.

For now, I have kept my clustering functions separately in rmsd_clustering.py. While there is some overlap with existing functions, my implementation utilizes different libraries and supports multiple backends, so I think for the time-being it is safer to have them apart.

Yes, I agree. I suggest making mention of the others in the docstrings to help with users getting confused otherwise.

I have attached the notebooks with the corrected functions

These notebooks use a filename for a file that I don't have so I won't be able to run them, but I will still glance over them seeing as the outputs are there and see if I have any more comments before finalising my review.

Can you put a smaller version of the file in tests/datafiles (maybe with different numbers 1, 2 and 3 frames per cluster for 3 clusters) and make these into unit tests? These would also serve as illustrative examples for users of how to run the code

It would also be good to convert some of these into tutorials for the prody website at some point although I don't know whether/when we can actually update it to include them. They also serve as useful findable tutorials on the prody website repo where they also render pretty well

jamesmkrieger · 2026-08-11T15:53:25Z

The notebooks generally look great. I have a few comments on them for improvement, but these do not generally impact the integrated code. One that does is that I see that all values in the tables in notebooks 04 and 05 including the part using showClusterStatisticsTable are represented as floats with several decimal places. I'd generally reduce actual floats to have 2 or 3 decimal places and convert those that actually aren't floats at all like Total Frames and Medoid Frame to int.

The approach comparison results in notebook 01 with results reproducing to residuals within 10^-7 to 10^-8 or 10^-15 are a perfect basis for being able to compare against a reference array from one of these calculations with numpy allclose like you have in notebook 02, which could also be a unit test.
Note the rmsd 0 at frame 1003 is just because that is the start of run 3 and matches the start of run 2 and the first residuals also reflect that.
Based on notebooks 01 and 02, it is clear that there are two clear clusters in run 2 that look like they could well be populated in trajectory 3 too, just with much shorter residence times and more interconversion. Can you map out which frames end up in which clusters and whether it matches the expectations?
Does the outlier cluster of 9 frames correspond to the starts of both runs that then don't persist in the hierachical case in notebook 03 just like these frames 0,1,2,1003 show up in the DBSCAN noise frames in notebook 05? The 1D RMSD plots in notebook 01 do show minimum RMSDs from the closer cluster in the range of 3.8 so they could well do. The cutoff is a bit arbitrary really but I'd say 3.8 removing these outlier frames is good.
It would be good to look at what's different about frames 605, 854 and 855 that make them outliers too. 605 looks like it's a brief deviation from the cluster in that time period. Perhaps 854 and 855 are at the transition between that cluster and going back into the other cluster at the end of run 2.
There's also a typo in the markdown cell underneath code cell 8 that prints the frames detected as noise and in a few other places. You should spell separated with an a after the p. There are some other typos around too
The number of clusters colour bar doesn't make sense. There should be fewer colours to match the matrix. You can truncate the tab10 cmap to just cover a smaller range of values
I think you worded it a bit negatively that DBSCAN tended to prefer one cluster but you justified it quite well. It's a bit surprising that your final result gives a much smaller cluster 2 than the other methods, so probably you want to comment on that too. I'd maybe pick epsilon 1.4 or 1.5 as that should split the two main clusters instead of just splitting off a minor one. You could take one of the cases with 3 clusters so the minor one is also captured and I'd say the high noise is a real biological noise that you don't want to optimise away.

jamesmkrieger

Here are a fair few changes, but it's looking really good now.

jamesmkrieger · 2026-08-11T14:00:10Z

+
+from . import rmsd_clustering
+from .rmsd_clustering import *
+__all__.extend(rmsd_clustering.__all__)


add back empty line

jamesmkrieger · 2026-08-11T14:02:52Z

+
+        rmsd_array = calcRMSD(ref_coords, target=aligned_coords)
+
+    rmsd_array = np.asarray(rmsd_array)


calcRMSD should already return an array from getRMSD

jamesmkrieger · 2026-08-11T14:05:59Z

+    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.")


can the code use matplotlib hist as a fall back rather than raising an error and asking the user to plot it themselves?

jamesmkrieger · 2026-08-11T14:06:55Z

+    user_label = kwargs.pop('label', None)
+
+    # Seaborn Defaults
+    kwargs.setdefault('bins', 50)


can the user provide these kwargs to override these defaults e.g. if there are much fewer structures so 50 bins is too many?

jamesmkrieger · 2026-08-11T14:07:47Z


            # Proceed with plotting
            plot_barh(result, bond_type, n_per_plot=n_per_plot, min_height=min_height)
-


revert this minimal change so this file doesn't show up as changed

jamesmkrieger · 2026-08-11T14:28:08Z

+    :arg trajectory: trajectory containing the coordinate sets to align
+    :type trajectory: :class:`prody.Trajectory`
+
+    :arg align: atom selection used to calculate the alignment transformation.


include selstr or something like that in the variable name so it's immediately obvious what it is without reading the doc string

jamesmkrieger · 2026-08-11T14:34:13Z

+        atoms.setCoords(orig_coords)
+        trajectory.reset()
+
+    return ref_coords, aligned_coords


also give an option of returning the aligned trajectory object itself

jamesmkrieger · 2026-08-11T14:34:32Z

+    return ref_coords, aligned_coords
+
+
+def calcPairwiseRMSD(aligned_coords):


It would be good to have a method like this for TrajBase like Ensemble has getRMSDs(self, pairwise=False), but I suppose that is a bigger change that we maybe don't want to make right now

jamesmkrieger · 2026-08-11T15:55:07Z

+
+    if show:
+        kwargs.setdefault('tablefmt', 'fancy_grid')
+        kwargs.setdefault('floatfmt', '.4f')


This default isn't right for all rows. Some of them are truly ints and most of the others don't deserve so much precision either.

kontheodosiadis added 2 commits August 4, 2026 17:55

New clustering functions - Hierarchical, KMedoids

b08824a

New clustering function - DBSCAN, OPTICS for Reachability Plot

eaf54d4

jamesmkrieger requested changes Aug 9, 2026

View reviewed changes

jamesmkrieger self-requested a review August 9, 2026 06:52

karolamik13 self-requested a review August 10, 2026 06:46

Created rmsd_clustering module - Relocated functions - Cleaned up Red…

b91f0e7

…undant Code

Fixed title bug in showPairwiseRMSDHeatmap

e09413d

jamesmkrieger requested changes Aug 11, 2026

View reviewed changes


		rmsd_array = calcRMSD(ref_coords, target=aligned_coords)

		rmsd_array = np.asarray(rmsd_array)


		# Proceed with plotting
		plot_barh(result, bond_type, n_per_plot=n_per_plot, min_height=min_height)

		return ref_coords, aligned_coords


		def calcPairwiseRMSD(aligned_coords):

Uh oh!

Conversation

kontheodosiadis commented Aug 7, 2026

Summary

Workflow

Key Features

Clustering Algorithms

Analysis Utilities

Visualization

Uh oh!

jamesmkrieger left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

kontheodosiadis commented Aug 10, 2026

Uh oh!

jamesmkrieger commented Aug 11, 2026

Uh oh!

jamesmkrieger commented Aug 11, 2026

Uh oh!

jamesmkrieger left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants