pubs.acs.org/jcim Article On Quality Thresholds for the Clustering of Molecular Structures Xavier Daura* and Oscar Conchillo-Solé Cite This: J. Chem. Inf. Model. 2022, 62, 5738−5745 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: It has been recently suggested that diametral (so-called quality) similarity thresholds are superior to radial ones for the clustering of molecular three-dimensional structures (González-Alemán et al., 2020). The argument has been made for two clustering algorithms available in various software packages for the analysis of molecular structures from ensembles generated by computer simulations, attributed to Daura et al. (1999) (radial threshold) and Heyer et al. (1999) (diametral threshold). Here, we compare these two algorithms using the root-mean-squared difference (rmsd) between the Cartesian coordinates of selected atoms as pairwise similarity metric. We discuss formally the relation between these two methods and illustrate their behavior with two examples, a set of points in two dimensions and the coordinates of the tau polypeptide along a trajectory extracted from a replica-exchange molecular-dynamics simulation (Shea and Levine, 2016). We show that the two methods produce equally sized clusters as long as adequate choices are made for the respective thresholds. The real issue is not whether the threshold is radial or diametral but how to choose in either case a threshold value that is physically meaningful. We will argue that, when clustering molecular structures with the rmsd as a metric, the simplest best guess for a threshold is actually radial in nature. ■ INTRODUCTION Over 2 decades ago, Heyer et al.1 developed an algorithm to When applied to molecular structures using the rmsd as a metric, both QTC and RTC scan a precalculated rmsd matrix cluster open reading frames (ORFs) based on their expression in search for the molecular configuration with the largest levels, using what they came to call jackknife correlation as number of neighbors satisfying the threshold, in an iterative pairwise similarity metric. The focus of the algorithm was to process that outputs a new cluster at each step and ends when find large clusters that had a “quality” guarantee, that is, a there are no further configurations to cluster. The difference minimum jackknife correlation between any two ORFs between the two algorithms lies in the nature of the threshold belonging to the same cluster. In other words, clusters would (diametral or radial) and the procedure to count the neighbors be guaranteed to have a maximum diameter, determined by a in. Recently, González-Alemán et al.2 pointed out that, because correlation threshold, ensuring the transitive property for the of their similarity, these two algorithms were often confused in relation ‘correlation > threshold’ between element pairs in a various software implementations commonly used in the field, cluster (i.e., if correlation(a, b) > threshold and correlation(b, misleading their users. They also evaluated the performance of c) > threshold, then correlation(a, c) > threshold, for any a, b, c the two algorithms by analyzing a trajectory of the tau belonging to the cluster). The algorithm was named quality polypeptide extracted from a replica-exchange molecular- cluster (QT_Clust) and it has since been used for the dynamics simulation previously published by Shea and clustering of several other types of data, including molecular Levine.5 In doing so, they used the same rmsd value for the structures.2 In this study, we will refer to this algorithm as thresholds applied in the two algorithms. They concluded that, QTC. due to its lack of a quality threshold, the RTC algorithm The same year, Daura et al.3,4 introduced an algorithm to (referred to as Daura’s algorithm in the paper) tends to cluster cluster molecular structures using the root-mean-squared unrelated configurations together and gave examples in which difference (rmsd) between the Cartesian coordinates of different QTC clusters were found as composing a single RTC selected atoms as pairwise similarity metric. The algorithm cluster. was meant to favor the most populated cluster and was radial in nature, that is, the rmsd threshold was applied from a configuration taken as a cluster reference, meaning transitivity Received: August 25, 2022 was not ensured for the relation 'rmsd < threshold'. In this Published: October 20, 2022 study, we will refer to this algorithm as RTC. The two algorithms are described under the Computational Details section and analyzed formally under the Results and Discussion section. © 2022 The Authors. Published by American Chemical Society https://doi.org/10.1021/acs.jcim.2c01079 5738 J. Chem. Inf. Model. 2022, 62, 5738−5745 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Clearly, the results observed by González-Alemán et al. had configurations) until no other configurations fulfilling the little to do with the quality of the thresholds and much to do diameter threshold are available (outer loop). with using the same rmsd value for a radial and a diametral The clusterings were performed with inhouse software threshold. While the relation 'rmsd < threshold' is not reproducing exactly the algorithms described here. Results transitive for the set of configurations conforming a cluster using established software implementations (available free of generated by the RTC algorithm, the alternative relation 'rmsd charge) on the tau-polypeptide example are provided as the < diameter' is. This leads to the following question: can we Supporting Information (SI) for comparison. For the RTC generate clusters with a predetermined maximum diameter case, results found in the SI were obtained using the using the RTC algorithm? In other words, can we choose the McLachlan algorithm6 as implemented in ProFit v3.3 radial threshold in the RTC algorithm in such a way that the (http://www.bioinf.org.uk/software/profit/) for the rmsd maximum diameter of a cluster will be equal to the threshold calculation and the RTC algorithm as implemented in we would use with the QTC algorithm? The answer is of HADDOCK v2.07 (cluster_struc, https://www.bonvinlab. course yes, if one would actually wish to do so. org/software/haddock2.2/) for the clustering. For the QTC In the following sections, we will present formally and case, results found in the SI were obtained using the analyze in detail the characteristics of the RTC and QTC implementation published by González-Alemán et al.2 methods. Although the two methods have been available and (https://github.com/rglez/QT). The results are exactly the heavily used since over 2 decades, a detailed analysis of their same, with small differences in the QTC case due to properties has not been published. We will show that it is implementation details that are explained in the SI document indeed possible to obtain equally sized clusters with the two and conform in both cases with the QTC algorithm. algorithms and will argue that this is in fact of little importance MD Simulation Data. The trajectory of the tau because, first, the threshold is an arbitrary quantity that may polypeptide was downloaded from https://github.com/ have different “ideal” values depending on the objective of the LQCT/BitQT/blob/master/examples/aligned_original_tau_ analysis and, second, for the purposes discussed here the 6K.dcd, together with the reference PDB file https://github. simplest physically based guide to decide on the value of the com/LQCT/BitQT/blob/master/examples/aligned_tau.pdb. threshold is in fact radial in nature rather than diametral. It corresponds to the exact same trajectory used by González- Alemán et al.2 in their comparison of the QTC and RTC ■ COMPUTATIONAL DETAILS Clustering Algorithms. The two algorithms, QTC1 and clustering algorithms. The trajectory contains 6001 config- urations of the polypeptide. We used the backbone N, H, Cα, RTC,3,4 require as input the matrix of rmsds between all pairs C, and O atoms of residues Lys2 to Asp11, that is, 50 atoms in of configurations. total, for least-squares fitting and rmsd calculation. The two In its mth iteration, the QTC algorithm scans the matrix in terminal residues, Gly1 and Leu12 and their capping groups, search for the configuration with the largest number of were left out because they are relatively free to rotate and neighbors in order to generate the mth cluster. Specifically, would introduce unnecessary noise in the clustering of the rest each configuration not clustered in the previous m − 1 of the structure. Likewise, we excluded the side chains because iterations (we shall refer to them as the available config- one generally focuses on the backbone to define a fold and side urations) will be considered both as a seed of a tentative chains would only introduce noise. This selection of atoms is cluster and as a potential neighbor of all other seeds. We use clearly different from that used by González-Alemán et al.2 (all the term tentative cluster to refer to a seed and its neighbors atoms), but this is irrelevant for the questions addressed here. before the neighborhood sizes of all seeds are compared to To generate points for the example in two dimensions (not select the actual cluster. For each seed, its neighbors will be that this is important), we simply took the x and y coordinates determined as follows: From all other available configurations, of the backbone N atom of Lys2 after least-squares fitting of all the one that upon its addition extends the diameter of the configurations to configuration number 2910. A subset of 1501 cluster the least, while fulfilling the condition that the diameter elements was then constructed by selecting 1 element every 4, must remain smaller than the threshold, is taken as the next starting with element 1. neighbor and included in the seed’s tentative cluster. This process is repeated until no remaining available configuration fulfills the threshold, at which point the tentative cluster for ■ RESULTS AND DISCUSSION Theoretical Framework and Properties. We note that that seed is complete. Note that within an iteration all available throughout this article we use the term diameter in its configurations are tested as potential neighbors of each one of generalized form, that is, as the largest distance between any the available seeds. Once the tentative clusters for all available two points on the boundary of a closed geometric figure (in seeds have been obtained, the one with the largest number of this case a cluster). Likewise, we use the term sphere as a short- elements is promoted to constitute the mth cluster, and all hand for (n − 1)-sphere, defined as the (n − 1)-dimensional elements of that cluster are removed from the pool of available boundary of an (n-dimensional) n-ball. configurations, thus finalizing the mth iteration. The algorithm n is stopped when there are no more configurations available or Let Sm = {x i (xi ,1 , ..., xi , n) : i Im} be a set of new clusters fall below a preset minimum number of elements. Euclidean vectors in a Cartesian frame (for convenience, we The RTC algorithm differs from the QTC one in the way shall also refer to xi as a point in that frame) representing the the neighbors of a seed are determined at each iteration: All Nm configurations of the molecule that are available for available configurations at a distance from the seed smaller clustering at the mth iteration of the RTC or QTC algorithm, than the threshold are taken as elements of the seed’s tentative where n is the number of coordinates that will be used for the cluster. Thus, the RTC algorithm avoids the double loop per rmsd calculation, Im = {i : 1 i N1, i Jm } is the seed that characterizes the QTC algorithm �to find the next set of indices of the elements of S1 that are available for element of the tentative cluster (inner loop over available clustering at the mth iteration and Jm, with J1 = 0̷ and 5739 https://doi.org/10.1021/acs.jcim.2c01079 J. Chem. Inf. Model. 2022, 62, 5738−5745 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Jm > 1 = {j : 1 j N1, x j Cl , 1 l m 1}, i s cluster), they are strictly not. This is because the collection the set of indices of the elements of S1 that have been already Dm(θ) defined above may contain more than one set. Indeed, clustered in previous iterations, where Cl stands for the cluster when using these algorithms in practice, at any given iteration set defined in iteration l (see below). it is relatively common to see two or more seeds tie as the ones Let xk ∈ Sm be the seed for a tentative cluster of elements of generating the tentative clusters with the largest number of Sm and Am,k(θ) = {xi ∈ Sm : rmsdki < θ} the set of elements elements (this is precisely why f is needed in the definition of within a rmsd-threshold θ ∈ R>0 from xk. Note that Cm(θ)), in which case both algorithms become configuration- order sensitive (for a given function f). An illustrative example rmsdki = x i x k (na) 1/2 , where na is the number of atoms of many candidate seeds forming tentative clusters of the same involved in the rmsd calculation (in principle, n = 3 × na). size in a given iteration is given in the Supporting Information Then, we define Am(θ) = {Am,k(θ) : k ∈ Im} as the collection of (see the comment on cluster 2 in the QTC case). It is true, such sets for all available seeds and Bm(θ) = {|Am,k(θ)| : k ∈ however, that ties tend to occur between seeds that are very Im}, where |Am,k(θ)| stands for the cardinality of Am,k(θ), as the close in space, thus having a relatively small impact on the collection of corresponding set sizes. clustering. In the RTC algorithm, Am,k(θ) is the tentative cluster We shall use the term cluster shape to refer to the convex “proposed” by seed xk, which shall be then compared to the hull of the set of points that belong to the cluster. Thus, the tentative clusters “proposed” by all other seeds. Thus, we cluster shape with maximum volume is for either algorithm a define Dm(θ) = {Am,k(θ) ∈ Am(θ) : |Am,k(θ)| = max(Bm(θ))} as sphere, of radius θ for RTC clusters and θ/2 for QTC ones. the collection of sets with the largest number of elements and Note that we can define three types of geometric centers for Em(θ) = {k ∈ Im : Am,k(θ) ∈ Dm(θ)} as the collection of both RTC and QTC clusters. The first type is the center of the corresponding indices. The mth cluster (output of the mth neighbor-search volume (a sphere) and corresponds to the iteration of the algorithm) is then defined as Cm(θ) = {Am,k(θ) position of the seed (as we have seen above, both algorithms ∈ Dm(θ) : k = f(Em(θ))}, where f is a function that returns one search within Am,k(θ)). For this reason, it is customary to refer element from a set, typically the function min(), in which case to the seed, particularly in the RTC case, as the central element Cm(θ) would be the set with the lowest index from those with of the cluster, but this is, as we shall see, misleading if we give it the largest number of elements. a spatial significance. The second type is the geometric center To impose the condition that the diameter of Cm(θ) is of the cluster shape, which will only coincide with the first one smaller than θ, as done in the QTC algorithm, we need to if the cluster is spherical (RTC case) or spherical and centered define a new set Fm,k(θ) = p(Am,k(θ)), where p is an element- on the seed (QTC case). And the third type is the geometric selection procedure, that is, Fm,k(θ) ⊂ Am,k(θ), such that rmsdij center of the cluster elements, which will only coincide with < θ, ∀ xi, xj ∈ Fm,k(θ). Fm,k(θ) is, in this algorithmic context, the the second one if the spatial distribution of points in the cluster tentative cluster “proposed” by seed xk, which shall be is homogeneous (which is rare). Thus, it is in fact common for compared to the tentative clusters “proposed” by all other the seeds of even RTC clusters to be relatively distant from the seeds. Thus, as done for the tentative clusters in the RTC case, geometric center of the cluster shape and/or the geometric we define Fm(θ) = {Fm,k(θ) : k ∈ Im} as the collection of such center of the cluster elements. sets for all available seeds and redefine Bm(θ) = {|Fm,k(θ)| : k ∈ The seed of a QTC cluster tends in fact to be close to the Im} as the collection of corresponding set sizes. Accordingly, cluster’s boundary. This is due to the specific procedure p that we redefine Dm(θ) = {Fm,k(θ) ∈ Fm(θ) : |Fm,k(θ)| = selects the elements of Fm,k(θ) from Am,k(θ). As described max(Bm(θ))} and Em(θ) = {k ∈ Im : Fm,k(θ) ∈ Dm(θ)}. The under the Computational Details section, at each step in the mth cluster is then defined as Cm(θ) = {Fm,k(θ) ∈ Dm(θ) : k = process of recruiting new elements for the tentative cluster f(Em(θ))}. Fm,k(θ), the point that minimizes the increase in diameter while Note that the latter cluster definition is not specific for the fulfilling the diametral threshold θ is selected. Thus, the first QTC algorithm but general for a group of diameter-based steps of the procedure are highly determined by the local algorithms. This is because, as it is defined, the procedure p is distribution of points around the seed, which is never equal in not unique. That is, the condition rmsdij < θ, ∀ xi, xj ∈ Fm,k(θ), all directions. This will introduce an early bias or dominant Fm,k(θ) ⊂ Am,k(θ) can be satisfied by different selection direction and sense for the cluster’s growth, which may be procedures p, leading to different subsets Fm,k(θ) of Am,k(θ). more or less prominent depending on the exact spatial For example, the tentative cluster Fm,k(θ) can be grown starting distribution of points. In cases in which the distribution from its seed following the procedure by Heyer et al.1 and induces a marked directionality (which could be relatively described under the Computational Details section (the frequent as we shall see in the 2D example below) and procedure p implemented in the QTC algorithm), or could assuming a spherical cluster, the cluster will tend to grow in be grown following a specific sequence of configurations (e.g., eccentric spherical layers away from the seed, leaving the seed time sequence), that is, testing at each step if the next at or very close to the cluster’s boundary. Note that it is the configuration in sequence satisfies the diametral threshold local distribution of points around the seed that primarily instead of searching for the configuration that minimizes the determines the direction and sense of growth of Fm,k(θ), rather increase in the cluster diameter (note that this would produce than the global spatial distribution of points in Am,k(θ). different Fm,k(θ) subsets for different configuration sequences). Therefore, Fm,k(θ) is not necessarily the subset of Am,k(θ) with Another variant of p could be searching for the subset Fm,k(θ) the highest cardinality. Another consequence of this is that, with the largest number of elements, and so forth. unlike for RTC clusters, for QTC clusters the relation |Cm(θ)| While the RTC and QTC algorithms are generally presented > |Cm+1(θ)| does not need to hold: the direction of growth of as invariant to permutation (referring to the order in which the the tentative cluster around a given seed xk may in some cases algorithm evaluates the seeds or the order in which the be less optimal for set Sm than for set Sm+1, so that we may have configurations are tested for inclusion in a seed’s tentative |Fm,k(θ)| < |Fm+1,k(θ)|, which can eventually lead to a situation 5740 https://doi.org/10.1021/acs.jcim.2c01079 J. Chem. Inf. Model. 2022, 62, 5738−5745 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 1. Representation of the first 20 clusters of 1501 points in 2D. (A) RTC algorithm, θ = 1.1 Å. (B) QTC algorithm, θ = 2.2 Å. The color sequence is the same in both panels, starting with dark green for the first cluster. The seed elements of the clusters are indicated with a black circle (with the area matching the cluster color). Points belonging to clusters with an index higher than 20 are indicated as black squares. (C) Detail of the process of generation of cluster 1 from panel B (dark green): to illustrate this process, the points are shown in four different colors following their sequence of inclusion in the cluster. Thus, initial growth from the seed (dark-green circle) is indicated as a dark-green trajectory. Elements in red, light green, and blue correspond to successive phases in the growth of the cluster, in this order. in which |Fm,k(θ)| < |Fm,j(θ)|, where Fm,j(θ) = Cm(θ), and underlying cluster topology (if such should physically exist). |Fm+1,k(θ)| > |Fm,j(θ)|, where Fm+1,k(θ) = Cm+1(θ). In practice, Therefore, in order to have equally sized clusters, the relation however, the fact that Cm(θ) is selected from a collection between the rmsd thresholds θr (RTC) and θd (QTC) will be Fm(θ) of overlapping tentative clusters Fm,k(θ) makes this in practice θd/2 ≤ θr < θd. potential inversion of cluster-size order very infrequent and we A more important question, however, is how to choose the have only observed it in two cases in the examples below, for threshold. Typically, the option of choice in the literature is very small (irrelevant) clusters. trial and error: both algorithms, particularly RTC, are To adequately compare the RTC and QTC algorithms in sufficiently fast that one can actually run them several times practical applications, one needs to keep in mind that θ is a with different θ values, until the outcome satisfies any chosen radial threshold (i.e., θr) in the RTC algorithm and a diametral criteria. The criterion is often visual, that is, the structures in a threshold (i.e., θd) in the QTC algorithm. If the spatial cluster “look the same”.2 While this may fit the purpose in distribution of points is such that the diameters of the RTC some types of studies, it is actually a weak criterion from a clusters are approximately equal to two times the threshold θr physical standpoint. As a general rule, before clustering data in at least one direction, the extreme case being spherical points in an n-dimensional space, one may want to query this clusters, one should use θr = θ/2 as the RTC threshold and θd space to gather information on the distribution of data points = θ as the QTC threshold for comparable results. However, as in it. This is also what makes physical sense in this case because we shall see in the tau-polypeptide example, the distribution of the space in question is the molecule’s configuration space points has rarely these characteristics when working with n- (reduced to the number of coordinates used for the rmsd dimensional data representing molecular configurations from calculation and the given ensemble sample). As will be shown computer simulation. First, the region of this n-dimensional for the tau-polypeptide example below, the simplest effective space corresponding to a conformer of the molecule (where way to query this space is by calculating the distribution of the term conformer may be taken as a very well-defined rmsd values. This can be done for the full (half) rmsd matrix, structure or a broader structural state, depending on the but the mixing of underlying distributions generally reduces its purpose of the clustering) will generally have an irregular shape informative value. Instead, we propose that a first tentative and its immediate surrounding may be void in many directions, clustering may be performed to calculate the rmsd distribution so that even if the algorithm tries to mix in points for each of the seed elements of the most populated clusters corresponding to neighbor conformers (when the threshold (as representatives of the high-density regions of interest), that is too big) in many spatial directions there might be simply is, taking the respective full columns (or rows) from the rmsd nothing to mix in. Second, the density of points in this n- matrix. In many cases, these distributions will show a first dimensional space is typically reduced in clustering exercises region of relatively high probability density, followed by a deep due to the selection of only one structure every so many and then a second, larger increase in probability density. This simulation steps, which has the effect of blurring the deep in the probability density corresponds to the end of the 5741 https://doi.org/10.1021/acs.jcim.2c01079 J. Chem. Inf. Model. 2022, 62, 5738−5745 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 2. rmsd distributions for the seeds of the six most populated clusters, after RTC clustering with six different thresholds. Each curve corresponds to the distribution of the rmsds between the given seed and all other 6000 configurations. The distributions are cut at 0.20 nm for clarity (the inset in panel D shows the full distributions from panel A as example). Panels A to F correspond to the distributions for the seeds of clusters 1 to 6, respectively. Each panel contains the distributions for six seeds, resulting from RTC clusterings with θr values of 0.12 (black), 0.13 (red), 0.14 (green), 0.15 (blue), 0.16 (orange), and 0.17 nm (brown). The dashed vertical lines show the positions of the thresholds (with colors matching the corresponding distributions). When two distributions overlap (i.e., clusterings with different thresholds produce the same seeds for the given cluster number), the overlapping curves are replaced by circles of the corresponding color. first layer of configurations around the seed and may therefore blue, in this order. The directionality of the growth, away from be interpreted as a (spherical) region of conformational the seed in eccentric circular layers, can be clearly observed. transition. After comparing the distributions for the seed Clustering of Structures from a MD Trajectory. elements of the main clusters, this information can be used to Following the strategy described above to infer a physically assess the threshold for a second, final clustering. Note that this meaningful threshold, we first used the RTC algorithm to information is radial in nature, and it therefore leads to a value perform a clustering of the 6001 configurations of the tau- for θr rather than θd. polypeptide, in order to identify points in the 150-dimensional In the next subsections, we illustrate some of the properties space (50 atoms × 3 coordinates) that are in regions of high discussed above of the RTC and QTC algorithms with two density. In this case, we chose the seed elements of the six examples. most populated clusters to then examine the distribution of Example Case in 2D. Figure 1 shows the clustering of the rmsds between each of these elements and the other 6000 set of 1501 points in two dimensions (see the Computational configurations. In order to see how the distributions may vary Details section) performed with the RTC (panel A) and QTC depending on the threshold used for this initial clustering, we (panels B and C) algorithms. This is a good example of a major chose six evenly spaced thresholds between 0.12 and 0.17 nm. limitation of fixed-size clustering algorithms: if the data has no The results are shown in Figure 2. In all distributions, an initial region corresponding to a first particular underlying structure or the threshold is completely layer of configurations around the seed can be distinguished, inadequate, these algorithms will still partition the data after which the probability density goes down to zero or close according to the chosen threshold. The question of how to to it. The distributions for the seeds of the first, most choose a threshold is, therefore, of particular significance, even populated clusters tend to be less sensitive to the threshold if we will ignore it in this first example. used for the clustering, as expected. This is not so much Note that, as mentioned under the Theoretical Framework because the clusters are more populated but because they are and Properties section, for densely populated spaces with few generated first, that is, the following clusters are affected by void regions, a selection of thresholds such that θr = θd/2, which configurations have or have not been already taken by where θr is the threshold used with the RTC algorithm and θd the previous ones. Although the rmsd value at which the is the threshold used with the QTC algorithm, produces probability density reaches zero differs for the different equally sized clusters in the two cases. As also mentioned, distributions, 0.17 nm stands out as a possible consensus while the seeds of the RTC clusters tend to be closer to the threshold: it is a point at which the probability density either centroids of their cluster shapes, the seeds of the QTC clusters reaches zero (notably for the seeds of cluster 1) or has not yet tend to be closer to the clusters’ boundaries. recovered significantly from zero. Panel C illustrates the generation process of QTC clusters, Based on these observations, we focused on the clusterings using the first cluster from panel B (dark green) as an example. performed with θr values of 0.12 nm, that is, the lowest value The initial steps are shown as a trajectory starting from the that seems adequate for some of the distributions shown in seed (dark-green circle). Consecutive phases in the growing of Figure 2, and 0.17 nm. To define the corresponding θd the cluster are illustrated with elements in red, light green, and thresholds for the QTC algorithm, we looked at the diameter 5742 https://doi.org/10.1021/acs.jcim.2c01079 J. Chem. Inf. Model. 2022, 62, 5738−5745 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 3. Relation between the distance to the seed and corresponding angle for every pair of elements of a cluster (excluding the seed). Specifically, rmsdij + rmsdik, for all xj, xk ∈ C1(θr), j, k ≠ i, where xi is the seed of C1(θr) and θr has the values 0.12 nm (A) and 0.17 nm (B), against the angle φ between the vectors xij and xik, where xij = xj − xi. Figure 4. rmsd distributions for the first 10 clusters, using the RTC algorithm with θr = 0.12 nm (black) and the QTC algorithm with θd = 0.20 nm (red). The number of elements in the cluster is for each case indicated. Figure 5. rmsd distributions for the first 10 clusters, using the RTC algorithm with θr = 0.17 nm (black) and the QTC algorithm with θd = 0.25 nm (red). The number of elements in the cluster is for each case indicated. of cluster 1 in each of these two RTC clusterings. Note that 0.278 nm), we decided to choose θd values of 0.20 and 0.25 cluster 1 is our best guide because its diameter is not nm for the QTC algorithm. conditioned by previous clusters. The diameter was 0.199 nm Note that, as already discussed, the diameter of an RTC for CRTC RTC 1 (0.12) and 0.278 nm for C1 (0.17). With these cluster should be generally expected to be smaller than 2θr for reference values and taking into account that the probability actual data sets from simulation, as confirmed by the values density for a rmsd distance above 0.25 nm in CRTC1 (0.17) is indicated in the previous paragraph. The reason for this is very small (only 17 elements at distances between 0.25 and illustrated in Figure 3. This figure shows, for every pair of 5743 https://doi.org/10.1021/acs.jcim.2c01079 J. Chem. Inf. Model. 2022, 62, 5738−5745 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article Figure 6. Distance of the seed from the center and cluster compactness (lines between points are drawn only to help visually distinguish the four data sets in each plot.) (A) rmsd between the seed element and the center of geometry of the elements in the cluster. (B) Radius of gyration of the cluster. Black circles (black solid line): RTC clustering with θr = 0.12 nm. Red circles (red solid line): QTC clustering with θd = 0.20 nm. Black squares (black dashed line): RTC clustering with θr = 0.17 nm. Red squares (red dashed line): QTC clustering with θd = 0.25 nm. elements (excluding the seed) of clusters CRTC1 (0.12) (panel A) illustrates that, except for cases in which a tight match between and CRTC 1 (0.17) (panel B), the relation between the sum of the rmsd distributions in Figures 4 and 5 exists, the RTC their rmsd distances to the seed and the angle between the clusters tend to be more compact (lower rog) than the QTC vectors associated with these distances. Thus, while there are clusters. pairs of points for which the sum of rmsds adds up indeed to When looking at the total number of clusters generated by 2θr, that is, 0.24 nm (panel A) and 0.34 nm (panel B), the the two algorithms, we see that they differ remarkably (see the angle between the corresponding vectors is in no case close to Supporting Information). Thus, while the number of clusters 180°, which precludes the maximum diameter from being generated by the RTC algorithm using θr values of 0.12 and reached. This figure also shows that, as expected (see panel A 0.17 nm is 1338 and 493, respectively, corresponding numbers in Figure 2), the point density is abruptly cut (solid wall at 0.24 for the QTC algorithm with θd values of 0.20 and 0.25 nm are nm) when using the 0.12 nm threshold (panel A), while a 599 and 276, respectively. However, focusing on the upper part threshold of 0.17 nm leads to a well-defined cluster (panel B). of the ranking list, we see that for RTC with θr = 0.12 nm the Figure 4 shows the distribution of rmsd values within each of number of clusters with 100 or more elements is 6 and these the first 10 clusters, for the RTC clustering with θr = 0.12 nm cover 19% of the overall population, while for QTC with θd = and the QTC clustering with θd = 0.20 nm. Seven of the 0.20 nm the number of clusters is also 6 and they cover 20% of clusters, including the first four, overlap almost perfectly. For the population. If we extend this to clusters with 10 or more three of the clusters, the QTC algorithm appears to have a elements, the numbers start to diverge, with RTC producing higher tendency to populate the far-right side of the 98 clusters that cover 50% of the population and QTC distribution. Figure 5 shows the corresponding distributions for the RTC producing 161 clusters covering 70% of the population. The clustering with θr = 0.17 nm and the QTC clustering with θd = trends are similar but the results less divergent for the 0.25 nm. It becomes here more apparent that the QTC comparison between RTC with θr = 0.17 nm and QTC with θd algorithm has a higher tendency to generate split distributions = 0.25 nm. Here, the RTC algorithm generates 8 clusters with and populate the far-right side of the distribution. For example, 100 or more elements covering 25% of the population and 149 it can be observed that an artificial cluster CQTC , containing clusters with 10 or more elements covering 82% of the 4 two different populations with similar weights, has displaced by population, while the QTC algorithm generates 7 clusters with one position in the ranking the QTC clusters that correspond 100 or more elements covering 25% of the population and 150 to CRTC 4 and CRTC 5 . clusters with 10 or more elements covering 91% of the We obtained for each cluster Cm, as shown in Figures 4 and population. Thus, while the upper part of the ranking looks 5, the center of geometry of the elements of the cluster, x(c) m , very much the same with the two algorithms, RTC produces and then calculated the rmsd between the seed element and many more small clusters at the end of the ranking. this point, as well as the radius of gyration of the cluster. Note How can we explain such large differences in the total that the radius of gyration of the mth cluster is here defined as number of clusters? The neighbor-search volume for the RTC algorithm is strictly spherical, while for the QTC algorithm it is 1/2 ij yz a volume within a sphere, with the diameter as the only shape j 1 2z rog m = jjj (x i m) z x (c) z zz restraint. This makes the QTC algorithm more flexible in jj |C | z terms of cluster shapes, allowing it to capture more of the k m xi Cm { (1) elements that would lay just outside a cluster boundary in the The results are shown in Figure 6. As anticipated, panel A RTC case. Thus, the RTC algorithm tends to generate many confirms that while the seeds of QTC clusters are in a majority more artificial clusters with orphan elements in the lower part of cases further from the center of geometry of the cluster of the cluster ranking. On the other hand, the greater shape elements than the seeds of RTC clusters, the latter are clearly flexibility of the QTC algorithm makes it have a higher off-center also. We therefore suggest to avoid, in either case, propensity to incorporate points from non-self neighbor the term central element to refer to the seed element. Panel B densities in a cluster, thus mixing in it different populations. 5744 https://doi.org/10.1021/acs.jcim.2c01079 J. Chem. Inf. Model. 2022, 62, 5738−5745 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article However, and as a conclusion, these differences between the two algorithms tend to be rather irrelevant in practice. ■ REFERENCES (1) Heyer, L. J.; Kruglyak, S.; Yooseph, S. Exploring Expression ■ DATA AND SOFTWARE AVAILABILITY The data used in this study was downloaded from https:// Data: Identification and Analysis of Coexpressed Genes. Genome Res. 1999, 9, 1106−1115. (2) González-Alemán, R.; Hernández-Castillo, D.; Caballero, J.; Montero-Cabrera, L. A. Quality Threshold Clustering of Molecular github.com/LQCT/BitQT/blob/master/examples/as indi- Dynamics: A Word of Caution. J. Chem. Inf. Model. 2020, 60, 467− cated under the Computational Details section. Although the 472. calculations shown here were performed with inhouse (3) Daura, X.; van Gunsteren, W. F.; Mark, A. E. Folding−unfolding software, these calculations can be performed with a variety thermodynamics of a β-heptapeptide from equilibrium simulations. of software packages implementing the RTC and QTC Proteins 1999, 34, 269−280. algorithms (see some of the potential choices in González- (4) Daura, X.; Gademann, K.; Jaun, B.; Seebach, D.; van Gunsteren, Alemán et al.2). We provide in the Supporting Information a W. F.; Mark, A. E. Peptide Folding: When Simulation Meets comparison between our results for the tau polypeptide and Experiment. Angew. Chem., Int. Ed. 1999, 38, 236−240. (5) Shea, J.-E.; Levine, Z. A. Studying the Early Stages of Protein results using one of the alternative software options in each Aggregation Using Replica Exchange Molecular Dynamics Simu- case. The results for the two software choices are exactly the lations. Methods Mol. Biol. 2016, 1345, 225−250. same, with small differences in the QTC case due to (6) McLachlan, A. D. Rapid comparison of protein structures. Acta implementation details that are explained in the SI document Cryst. A 1982, 38, 871−873. and conform in both cases with the QTC algorithm. (7) de Vries, S. J.; van Dijk, A. D. J.; Krzeminski, M.; van Dijk, M.; Thureau, A.; Hsu, V.; Wassenaar, T.; Bonvin, A. M. J. J. HADDOCK ■ ASSOCIATED CONTENT * Supporting Information sı versus HADDOCK: New features and performance of HAD- DOCK2.0 on the CAPRI targets. Proteins 2007, 69, 726−733. The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jcim.2c01079. Tau_results_diff_software.xlsx: listing of clusters, in- cluding seed element and number of elements, from the tau-polypeptide example, with θr values of 0.12 nm and 0.17 nm (RTC clustering) and θd values of 0.20 nm and 0.25 nm (QTC clustering), corresponding to the data presented here (inhouse software) and computed with alternative available software, as indicated under the Computational Details section (XLSX) ■ AUTHOR INFORMATION Corresponding Author Xavier Daura − Catalan Institution for Research and Advanced Studies (ICREA), Barcelona 08010, Spain; Institute of Biotechnology and Biomedicine, Universitat Autònoma de Barcelona, Cerdanyola del Vallès 08193, Spain; Centro de Investigación Biomédica en Red de Bioingeniería, Biomateriales y Nanomedicina, Instituto de Salud Carlos III, Cerdanyola del Vallès 08193, Spain; orcid.org/0000-0001-9235-6730; Email: xavier.daura@ uab.cat Author Oscar Conchillo-Solé − Institute of Biotechnology and Biomedicine and Department of Genetics and Microbiology, Universitat Autònoma de Barcelona, Cerdanyola del Vallès 08193, Spain; orcid.org/0000-0003-4266-246X Complete contact information is available at: https://pubs.acs.org/10.1021/acs.jcim.2c01079 Notes The authors declare no competing financial interest. ■ ACKNOWLEDGMENTS The authors thank Joan-Emma Shea and Zachary Levine for permission to use their tau-polypeptide simulation trajectories. This work received financial support from the Spanish Ministry for Science and Innovation (grant PID2019-111364RB-I00). 5745 https://doi.org/10.1021/acs.jcim.2c01079 J. Chem. Inf. Model. 2022, 62, 5738−5745