RoG-adaptive conformational clustering of MD trajectories

Automatic, radius-of-gyration (RoG) adaptive clustering of molecular-dynamics conformations. The trajectory is sliced into fixed-width RoG bins, each bin is clustered independently with hierarchical clustering (Ward.D), and the collected representatives are refined with k-means. The result is a small set of representative structures (frames) that span the conformational landscape, weighted by population so that densely sampled regions get more clusters.

The number of bins and the number of clusters are not hand-picked - they adapt to your data from a single, intuitive parameter (the RoG bin width).

Origin. This is an automated, RoG-adaptive extension of the clustering scripts from the gEDES project (Vargiu lab): https://github.com/vargiulab/gEDES/tree/main/2_clustering_scripts. The gEDES version sets the window count to max(10, RoG_range / 0.2), asks the user for the total cluster number, keeps a separate script and R template per CV count, and maps structures to CVs with a linear grep loop. This version instead: exposes the RoG bin width as the single knob (default 0.1 A, N_bin = ceil(RoG_range / bin_width), no lower floor); derives the cluster count automatically and allocates it proportionally to each bin's population; auto-detects the number of CVs (one script for 4 / 6 / 8 ...); embeds its R templates; replaces the linear grep search with a hash-table (awk associative-array) lookup for speed; and handles sparse / empty bins. Please credit gEDES if you use this (see the Acknowledgements below).

---

What it does

1. Bin by RoG. The RoG range of the trajectory is divided into bins of a chosen width (default 0.1 Å). The number of bins adapts automatically:

   N_bin = ceil( (RoG_max - RoG_min) / bin_width )
   

2. Cluster each bin with hierarchical clustering (hclust, Ward.D) on the collective variables (CVs). Each bin receives a number of clusters proportional to its population, with a floor of minclwin per bin, so no populated RoG region is lost and dense regions get more detail.

3. Refine with k-means seeded from the per-bin representatives, and extract the representative frame of each final cluster as a PDB.

4. Align & report. Representative frames are aligned (VMD) and their RoG / RMSD-from-reference are tabulated and sorted.

The number of collective variables (4, 5, 6, 8, …) is detected automatically from the input - one script handles any CV count.

---

Requirements

Tool	Used for
bash	driver script
R (packages fpc, cluster)	hierarchical + k-means clustering
PLUMED	computing CVs from the trajectory (gen_colvar_meta*.sh)
VMD	aligning representative frames (trajanalyze_templ.tcl)
GROMACS (gmx) or catdcd	extracting frames as PDB (lines_templatekmeans.sh)

---

Repository layout

Generic (target-independent - use as-is):

• R_clustering.sh - the main driver. The R templates (hierarchical + k-means) are embedded inside it and written out at runtime, so no separate .R files are needed.
• clean.sh - removes intermediate/output files from a previous run.

Target-specific templates (copy into a per-target working directory and edit):

• gen_colvar_meta.sh (_6CV.sh, _8CV.sh) - runs PLUMED and builds the input files (index_RoGBS.dat, CVS_clustering.dat) for 4 / 6 / 8 CVs.
• plumed-common.dat - PLUMED input defining your CVs. Ships with the public gEDES example (8 CVs: binding-site RoG + 3 pocket + 4 contact CVs); replace the atom numbers with your system's. RoG must stay the first printed CV.
• lines_templatekmeans.sh - extracts the selected frames as PDBs (set your trajectory + reference structure here).
• trajanalyze_templ.tcl - VMD script that aligns frames and writes RoG_trajectory.dat / RMSD_fromholo.dat (set your reference selection here).

Provided by you per target (not committed - they are large / system-specific):

• the trajectory (*_concat.dcd or .xtc), holo_processed.pdb, initial.pdb, and the COLVAR_apoMD produced by PLUMED.

---

Input files (in the working directory)

Produced by gen_colvar_meta*.sh:

• index_RoGBS.dat - frame_index RoG [other CVs...]. Column 2 (RoG) is what the binning uses.
• CVS_clustering.dat - frame_index cv1,cv2,...,cvN. The comma-separated CV block; the CV count N is auto-detected from here.

---

Usage

./R_clustering.sh [minclwin] [bin_width]

Both arguments are optional. Run ./R_clustering.sh -h for the built-in help.

Argument	Meaning	Default
minclwin	Minimum clusters per RoG bin. The main knob for how many clusters you get: every bin yields at least this many representatives; denser bins get proportionally more. Larger → more clusters.	1
bin_width	Width of each RoG bin, in Å. Smaller → more, finer bins. The number of bins adapts to the data.	0.1

Examples

./R_clustering.sh            # minclwin=1, 0.1 Å bins  (default)
./R_clustering.sh 3          # minclwin=3, 0.1 Å bins  (more clusters)
./R_clustering.sh 1 0.2      # minclwin=1, 0.2 Å bins  (coarser)
./R_clustering.sh 2 0.05     # minclwin=2, 0.05 Å bins (finer)

Typical workflow

./gen_colvar_meta.sh         # (edit for your trajectory first) -> builds inputs
./clean.sh                   # clear any previous run
./R_clustering.sh 1 0.1      # run the clustering

---

Output

A folder named for the actual result, e.g.:

frames_hierarchkmeans_13clusters_4CV_bw0.1_ward/
├── <frame>.pdb                 # one representative structure per cluster
├── frames.pdb                  # all representatives concatenated
├── clusters_aligned/           # VMD-aligned representatives
├── RoG_trajectory.dat          # RoG of each representative
├── RMSD_fromholo.dat           # RMSD-from-reference of each representative
├── sorted_RoG_trajectory.dat   # the above, sorted ascending
└── sorted_RMSD_fromholo.dat

The folder name encodes: number of clusters, CV count, and bin width.

---

How the numbers are decided

• The RoG bin width sets how many bins there are (N_bin = ⌈Δ/bin_width⌉).
• Each bin's cluster count is round(N · pop_bin / pop_total) with a floor of minclwin, where N = minclwin × N_bin. Denser bins → more clusters.
• Per-bin cluster counts are summed, pooled, deduplicated, and k-means is seeded from the survivors - so the final cluster count is data-driven and typically differs from any single input number.
• Bins with only 1-2 structures are handled directly (a lone structure is its own representative); empty bins are skipped.

---

License

MIT - see LICENSE.

Acknowledgements / citation

This pipeline builds directly on the gEDES clustering scripts from the Vargiu lab: https://github.com/vargiulab/gEDES. The clustering method (binding-site RoG + contact CVs, hierarchical Ward.D per RoG window, k-means representative selection) and the example plumed-common.dat originate there; this repository adds the RoG-adaptive binning, population-proportional cluster allocation, automatic CV-count detection, and template embedding.

If you use this in published work, please cite the gEDES project and this repository, together with the papers where the approach was developed and applied:

Key references

1. Toward the Prediction of Binding Events in Very Flexible, Allosteric, Multidomain Proteins. J. Chem. Inf. Model. https://doi.org/10.1021/acs.jcim.4c01810

2. Holo-like and Druggable Protein Conformations from Enhanced Sampling of Binding Pocket Volume and Shape. J. Chem. Inf. Model. 2018. https://doi.org/10.1021/acs.jcim.8b00730 (EDES method)