Trajectory Clustering, Statistical Analysis, and Visualization Tools

[kontheodosiadis]

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:

1. Frame-to-Reference Analysis

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

2. Frame-to-Frame Analysis

• Generate pairwise RMSD distance matrices across trajectory frames.

3. Clustering

• Apply clustering algorithms on the pairwise distance matrix.

4. Analysis and Visualization

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

5. 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.

1. 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.

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

[jamesmkrieger]

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

[karolamik13]

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

[jamesmkrieger]

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

[kontheodosiadis]

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]

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]

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]

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

[jamesmkrieger]

also give an option of returning the aligned trajectory object itself

[jamesmkrieger]

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]

you could increment all values to start at 1 and warn rather than skipping them

[jamesmkrieger]

I think we're more or less in a good state now

[karolamik13]

I tested all the functions. Working great.