A complete, multi-level conformational clustering of antibody complementarity-determining regions

Acquisition of antibody structure files

The three-dimensional coordinates of all antibody structures were downloaded from the PDB (Berman et al., 2000). Since the presence of antibody variable chains inside a PDB file is not annotated in a unique and systematic way, the advanced search tool of the database was used in order to apply composite search filters. The simple text search query of the database with the keywords “antibody” or “immunoglobulin” returns hundreds of unwanted PDB files, for example those that only contain a constant antibody fragment (Fc) or those that contain the keyword in their primary citation without any relevant structures in the file. Conversely, in several cases, antibody variable chains (Fv) are found in PDB files that do not contain the keywords “antibody” or “immunoglobulin” at all. In order to refine the obtained results, multiple queries were run using a variety of relevant keywords and their combinations with appropriate logical AND/OR/NOT connectors. The keywords employed typically included: “antibody”, “immunoglobulin”, “Fab”, “Fv”, “Fc”, “light chain”, “heavy chain”, “intact”, “complete”, “camelid”, “llama”, “VHH”, “light dimer” and “Bence -Jones”.

The final dataset comprised of exactly 1,351 PDB structure files, 8 of which contain variable chains from two different antibodies (5 were idiotype-anti-idiotype complexes), increasing the total number of antibody structures to 1,359. The total number of included CDRs is 12,712, 2,827 of which are unique in sequence. Table 1 contains a summary of the dataset contents. The dataset was locked on the 31st of December 2011 and should reflect the complete repertoire of antibody CDR structures up to that date. The set should be complete, given the proviso that there was a lack of specific tagging or annotation in the required PDB files.

Numbering of antibody variable chains and definition of CDR extents

All the antibody variable chain sequences in the dataset were structurally numbered in order to detect the beginning and end of each CDR, using regular expressions for the detection of the location of conserved sequence patterns. The initially adopted numbering scheme was the Chothia scheme (Chothia & Lesk, 1987) because it correctly places the insertion points in CDR-L1 and CDR-H1, but also because it is very frequently used in the CDR-related literature. The definitions used for the extents of CDRs-L1, -L2, -L3 and -H3 were also those established by Chothia & Lesk (1987) because they are most commonly used. However, for CDR-H1 and CDR-H2, the definitions adopted were those used in North, Lehmann & Dunbrack (2011). Based on previous experience from the visual examination of CDR-H1 structural superpositions, it was noted that the N-terminal portion of the loop where Kabat’s (Kabat et al., 1991) and Chothia’s CDR-H1 differ shows great variability both in sequence and structure. Thus, it was judged that this cluster analysis would be more revealing and useful if the CDR-H1 extent was considered as the entire length of the loop, namely residues H23–H35. As far as CDR-H2 was concerned, it was observed that the C-terminal portion of Kabat’s definition (i.e., residues H59–H65) remained relatively unchanged conformationally in most CDRs. Therefore, only the length of the symmetrical loop portion between residues H50–H58 was retained for the CDR-H2 definition.

CDR length and numbering scheme amendments

A number of antibodies contained a CDR with more residues than the current scheme could accept. The CDRs concerned were CDR-L2, -L3, -H1, -H2 and -H3. These CDRs, except for CDR-L2, already contained an insertion locus so the maximum allowed length was extended by adding more insertion positions (letters) to the numbering scheme. An insertion point was required in CDR-L2 for an 11-residue length. By superposing the new 11-residue loop (PDB code 3FFD) on a typical 7-residue one (1A4K), it was strongly suggested that the insertion point in CDR-L2 should be placed at position L51 (Fig. 1).

Two more cases required intervention in the numbering scheme. The first was in Light chain framework-3 (LFR3), where structure 1PW3 showed a 2-residue insertion. Superposition of this structure to the respective portion of a typical Light variable chain (1A4K) revealed that an insertion point should be introduced at position L67 (Fig. 2). The second case was raised by two anti-HIV antibodies observed in structures 3RPI and 3SE8, showing an insertion of 3 and 7 residues respectively in Heavy chain framework-3 (HFR3). Superposition of these frameworks onto a typical HFR3 (3MLY) suggested that an insertion point should be placed at residue H74 (Fig. 3). Table 2 summarises all the amendments brought to the initial numbering scheme in order to accommodate the special cases discovered in the dataset.

Clustering overview

In order to increase the usefulness of the clustering result in a way that meets the needs of a wider range of applications, a novel three-level nested cluster architecture was devised. At the parent-level, members of the same cluster share the least similarity in terms of Cα-atom Root-Mean-Square Deviation (RMSD), as the cluster is designed to include all the variants of a conformational theme within the limits of a statistical cluster validation. At the daughter-level, RMSD variance is successively reduced and members of the same cluster are increasingly similar. This stratified scheme could also be perceived as a variation of sensitivity to the potential natural flexibility of a CDR conformation (looser clusters), as well as a trade-off to the specificity of a particular shape (tighter clusters).

First-level clusters were formed by the use of a statistical clustering method, while second- and third-level clusters were defined using qualitative criteria. More specifically, the data was initially analysed by average- and complete-distance hierarchical clustering using RMSD distance matrices, and pruning of the resulting trees was performed with the Dynamic Tree Cut algorithm (Langfelder, Zhang & Horvath, 2007). RMSD distance matrices were obtained by performing all-by-all Cα-atom superpositions of the entire CDR loops, per individual CDR length. The result of hierarchical clustering was a set of level-1 structural classes, as traditionally produced by various methods in all previous CDR conformational studies, meaning that members of the same cluster were similar to a degree that is defined by the tree-pruning and clustering criteria.

Subsequently, ϕ/ψ angles were calculated for all CDR residues, each residue was attributed to a Ramachandran region and Ramachandran logos were formulated for each CDR. For practical and computational reasons, the boundaries of the different Ramachandran regions were based on the Ramachandran Plot subdivision used by North, Lehmann & Dunbrack (2011) (Fig. 4). Two types of Ramachandran logos are defined for each CDR, namely one where similar conformational regions were represented by the same letter (also suggested in North, Lehmann & Dunbrack, 2011), which will henceforth be called the reduced-Ramachandran Logo or r-RL, and one where every conformational region is represented individually, called the full-Ramachandran Logo or f-RL. For the formation of level-2 clusters, the members of any given parent level-1 cluster were regrouped by identical r-RL, meaning that members of the same cluster contain residues at each CDR position that belong to similar conformational regions. For the formation of level-3 clusters, the members of any given level-2 cluster were regrouped by identical f-RL, meaning that members of the same cluster contain residues at each CDR position that belong to the exact same conformational region. An example showing the layout of this nested cluster architecture can be seen in Fig. 5. Outliers/singletons were all given the tag ‘-O-‘ in their conformational logo, which created a common parent class that allowed the subsequent formation of 2nd- and 3rd-level clusters within outlier space, as well.

Clustering method

The RMSD distance matrices produced for each CDR/length combination were used for hierarchical analysis in the statistical package RGui (GNU project, http://www.sciviews.org/_rgui/). The average-linkage and complete-linkage algorithms were preferred to single-linkage in order to avoid chaining effects in dense configurations of the dataset in conformational space, and were both explored for every CDR/length combination. Hierarchical trees (dendrograms) that gave a Cophenetic Correlation Coefficient (CPCC) lower than 0.6 were directly discarded as pointing to poor fitting of the data. In all cases at least one of the hierarchical methods achieved a CPCC score greater than 0.6. Both hierarchical trees were considered whenever the CPCC was acceptable and comparatively evaluated using the criteria below.

The Dynamic Hybrid Tree Cut method of the Dynamic Tree Cut statistical package in RGui was utilised for dendrogram pruning. The package has been previously successfully used for the detection of biologically meaningful clusters in a protein–protein interaction network in Drosophila (Dong & Horvath, 2007). The Dynamic Hybrid Tree Cut algorithm offers flexibility, by allowing the user to set the desired pruning parameters for cluster and outlier recognition. Specifically, the algorithm defines four cluster shape criteria: (1) the minimum number of cluster members (N₀, minClusterSize), (2) the maximum scatter of the pairwise distances between the lowest merged objects (CDR structures) in each cluster, called the cluster core (d_max, maxAbsCoreScatter), (3) the maximum joining height at which a cluster attaches to the rest of the dendrogram (h_max, cutHeight), and (4) the minimum distance between the core scatter and the joining height of a cluster to the dendrogram, called the cluster gap (g_min, minAbsGap). The core scatter is defined as the average of all pairwise dissimilarities between objects belonging to the core of the cluster. Consequently, a branch is considered a cluster when it contains a minimum number of members (N₀), its joining height is at most h_max, its core is tightly connected (d_max) and distinct from its neighbourhood (g_min). Specifically, the minimum cluster gap distance (g_min) can be perceived as the minimum allowance for the cluster to expand its diameter from its core until it reaches a neighbouring cluster.

Although these pruning parameters are explained in depth in the corresponding method paper (Dong & Horvath, 2007), an example of the application of pruning parameters to an actual dendrogram from this analysis can be seen in Fig. 6. The number of objects assigned to the core of a cluster is derived from the following implemented formula: (1) ${\displaystyle n_c=min\left\{\int \left(N_0/2+\sqrt{N-N_0/2}\right),N\right\}}$ with n_c the number of core objects, N₀ the defined minimum cluster size and N the total number of objects in the cluster. As a consequence, the core of small clusters can be as large as the whole cluster, while the core of large clusters remains a fraction of the lowest joined objects.

The algorithm examines the dendrogram in a bottom-up manner and attempts to perform three types of branch merges: a merge of two singletons which creates a new branch, the addition of a singleton to a branch, or a merge of two branches. In each step two branches are tested against the pruning criteria: if both considered branches satisfy the criteria then both are declared “closed” and no further objects are added in the current step. Otherwise, the branches are merged and this new group is reassessed for cluster conformity during the next merge with an adjacent branch. Objects too far from a cluster are left unlabelled as outliers. Once all possible object assignments are performed, the method allows a further optional ‘Partitioning Around Medoids-like’ step (PAM). During this step, unlabelled objects (outliers) are considered one-by-one and are assigned to existing clusters based on a user-defined maximum allowable distance, or when their distance is smaller than the cluster’s radius. There are two options available for the cluster radius definition (parameter: useMedoids[=FALSE/TRUE]). If average distances are being used (FALSE), then the radius of the cluster is defined as the maximum of the average distances between objects in the cluster. If instead medoids are used (TRUE), then the radius is defined as the maximum distance of the cluster’s medoid to the cluster’s objects.

In order to detect the pruning parameters that lead to the best clustering result, an R routine was created which cycles the pruning method through a range of h_max, then g_min, then d_max using 0.1 increment steps. In each step, the quality of the clusters was assessed by calculation of the average Silhouette Coefficient (SC) and a cut-off of 0.51 was defined as the minimum required coefficient value for a reasonable structure to be found. The minimum number of members per cluster (N₀) was set to 2, in order to make sure that true singletons that could not form a compact cluster core with sufficient separation from neighbouring clusters were left as outliers. The output of this routine returned the clustering parameters, the number of clusters and outliers, the average SC and an auxiliary index showing the ratio of outliers over clusters.

Multidimensional scaling was applied to all distance matrices and 2D maps were produced for visual inspection of the clusters. In addition, 3D maps were created and consulted through the visualisation tool GNUPLOT (Williams et al., 2007–2011), for better perception of the configuration of the global population of each CDR/length combination. The 2D/3D maps and the respective Silhouette Plots of pruning results with average SC greater than 0.51 and all positive individual Silhouette Widths (SW) were consulted in all cases in order to continually have a visual appreciation of the data configuration and clustering evolution, and to make informed decisions which allowed the final formalisation of the clustering procedure. Given that the desired clustering result would ideally produce as many well separated clusters and as few outliers as possible, the auxiliary index offered a quick composite comparison between pruning results, and was defined as: (2) ${\displaystyle a=\left(1+S\right)/C}$ where S is the number of outliers/singletons and C the number of clusters. The unit (1) was added to the index’s numerator in order to allow the comparison between pruning results with 0 outliers/singletons, but a different number of clusters.

Another index employed during the clustering procedure was that of the ideal maximum cluster diameter, which took into consideration the examined CDR length (l): (3) ${\displaystyle D_i=1+\left(\frac{l-9}{10}\right).}$ The rationale behind this formula was to define an ideal maximum diameter by adding or subtracting 0.1 Å per residue respectively above or below a length of 9. For a CDR with 9-residues, this diameter was set empirically at 1.0 Å, based on experience of manual 3D superpositions of CDR-L3/9-residues with the graphics program Swiss-PdbViewer (Spdbv; Guex & Peitsch, 1997). Observations suggested 1.0 Å to be an appropriate cut-off for significant visual conformational similarity for CDRs of this length. This auxiliary index played no further analytical role than to merely define a cut-off at which the possibility of cluster splitting was to be explored during the clustering procedure. In no case did it impose a diameter threshold for cluster formation. Conversely, cluster merging was explored between clusters that contained one or more members with greater affinity for the second cluster (revealed by its negative SW). If the merge resulted in a global average SC ≥ 0.51 then it was retained, otherwise the entire partition was discarded. In the end, the preferred clustering parameters were those that resulted in global average SC equal or higher than 0.51, all positive individual SWs and the lower auxiliary index α (Eq. (2)). If the number of outliers remained high, the optional PAM-stage was applied at the end of the tree cut procedure, but its results were only retained if all of the above partition quality criteria were satisfied.

When the optimal clustering result was obtained, the clusters’ cores, medoids, most distant members and their diameters were extracted for that CDR/length combination via a dedicated R routine. Clustering summaries were created with Java code, as well as lists and various post-analytical data that are detailed later.

Clustering results

Tables of results were constructed for 58 CDR/length combination, gathering information that describes each individual cluster, which can be consulted for quick reference (Tables 3–7 for CDR-L1/-L2/-L3/-H1/-H2 and a separate supplementary table for CDR-H3, Supplemental Information 4). A summary table with all clustered lengths is available in Table 8. Detailed membership assignments can be found in two forms: one where every CDR is shown in alphabetical PDB order with all available clustering and data-mined information (cis/trans peptides, structure resolution, crystal space group, sequence, Ramachandran logos, cluster core label, bound state, light isotype, heavy or light chain only) and one where the same information is given in cluster order (Supplemental Information 6 and Supplemental Information 5). The ω-angle cut-off for cis-peptide detection was set to ± 30°; absence of cis-content that satisfied these limits resulted in an all-trans (allT) label. Bound state was flagged based on a list of bound antibodies obtained from SAbDab (Dunbar et al., 2014). This list did not contain idiotype-anti-idiotype complexes, therefore the 5 such files in the dataset were additionally flagged as bound (entries 1CIC, 1DVF, 1IAI, 1PG7, 3BQU).

Comparison of clustering results

The level-1 clusters obtained in this work were compared to the clustering results of previous major CDR studies (Tables 9–13 for CDR-L1, -L2, -L3, -H1 and -H2, Supplemental Information 2 for CDR-H3). Specifically, comparisons were made with the clusters found in Martin & Thornton (1996) because it was the first five CDR clustering performed on a significant CDR dataset (57 antibody structures, 269 CDRs), presented most major conformational classes and for these reasons is regularly cited in research of this kind. Comparisons were also made with the clustering results in North, Lehmann & Dunbrack (2011) as this is the most recent relevant analysis, which used the largest CDR dataset (932 antibody structures before filtering, 1897 CDRs after filtering) until the present study. Also included were the results from Kuroda et al. (2009) for the comparisons in CDR-L3, as this recent dedicated analysis used an RMSD-based approach, as is the case in this work, while using a considerable number of CDR structures (212 CDR-L3 structures). For the first five CDRs, the present study comprised 1,359 antibody structures and 10,680 CDRs (and a total of 12,712 CDRs including CDR-H3). Commenting on these comparisons is made in the discussion section below.

Rogue clusters and sequences

Assigned as ‘rogue’ were two conformational clusters that contain one or more members with identical CDR sequences. This definition was first used for CDR conformations by Martin & Thornton (1996) with respect to their unpredictability by canonical sequence templates when all their key residues are overlapping. In this work there is an expansion of this notion with the term ‘rogue CDR sequences’. This refers specifically to those identical sequences that are found to exist with more than one distinct conformation. The extraction of such sequences allows for further investigation, which can reveal any particular circumstances or neighbouring sequence features that led to a different CDR conformation despite the identical sequence. For example, examination of antibody Fvs with rogue CDR sequences may reveal the influence of neighbouring main-chain atoms, a particular framework residue influencing the CDR conformation, a conformational switch due to interface interactions (e.g., with an antigen), intrusive crystal-packing interactions, or even suggest some experimental error.

All cluster populations were parsed for rogue CDR sequences and a list of CDRs, tagged by their cluster assignment, was created for future detailed analysis (Supplemental Information 1). Also in the same file, entries with completely identical Fvs which belong to different conformational clusters (full-chain rogues) are reported separately, while entries containing bound antibodies are flagged as such by an asterisk. Furthermore, cluster populations were compared in all CDR/length sets, and the minimum number of amino acid differences, position-by-position, was calculated between any two sequences of different clusters. This difference was termed the ‘minimum pairwise Sequence Distance between clusters’, or mSD (essentially a minimum Hamming distance between sequences). Matrices showing the mSD between all clusters were constructed for every CDR/length, and heatmaps were produced in order to allow a quick visual appreciation of the degree of sequence dissimilarity between clusters (Supplemental Information 3). The purpose of these heatmaps is to assist mutation studies by promptly directing the researcher to clusters/CDR sequences of interest, as well as sequence-to-structure studies by biologists or modellers.

Investigation of structure resolution in outlier space

As a preliminary layer of quality assessment for the outliers in the present clustering, the min, max, average and median resolutions were calculated in clustered and outlier spaces per CDR/length (-L1, -L3, -H1, -H2, being of the highest interest). These values were plotted as stock charts for comparison, in order to observe any global correlation between the outlier space content and possibly erroneous CDR structures due to poor resolution (Supplemental Information 7). In only four cases (CDR-H1/15-, CDR-H1/16-, CDR-L1/12- and CDR-L1/16-residues) was the median resolution of outlier space found to be more than 0.5 Å higher than the respective median in clustered space, and in only two cases (CDR-H1/15-residues, 3 outliers in total, and CDR-L1/12-residues, 1 outlier in total) was the outlier median resolution value above 2.8 Å. In conclusion, average structure resolution does not appear to be a determinant factor of the outlier content, although it remains possible that wrong structures due to poor resolution may exist between the outliers. In fact, as proposed throughout this work, any decisions on structure validity should be considerably easier to make during targeted analysis of the structures/clusters of interest, when using the results of the present clustering. The supplementary file (Supplemental Information 7) also contains complementary bar charts showing the percentages of bound content in outlier and clustered space.