Exploring structural diversity across the protein universe with The Encyclopedia of Domains

1The AlphaFold Protein Structure Database (AFDB) is a groundbreaking initiative which significantly broadened the protein structure universe by expanding 3D representation to over 200 million UniProt sequences. The implications of the AFDB have been profound, not only in academic research in the life sciences but also in the commercial sphere, where the integration of novel data and the advanced technologies used to generate accurate structures is being explored for next-generation structure-based drug discovery.

2Notwithstanding its revolutionary impact, the AFDB is not without its limitations and presents a new set of challenges. The sheer scale of the data makes many traditional tools and pipelines, originally designed for considerably smaller datasets, inadequate for navigating the extensive number of structures and sequences within. This necessitates an evolved strategy, warranting a new perspective on how best to represent and traverse the data and complex relationships within such an expansive database, as well as requiring synergy between new algorithmic methods and computational hardware. Recent studies explored AFDB by partitioning full-length AFDB models into clusters of structurally similar proteins, as well as characterising their functions. At a more granular level, Bordin et al. and Schaeffer et al. interrogated the composition of specific proteomes (the initial AFDB 21 model organism dataset, and 48 proteomes, respectively), cataloguing domains under the CATH and ECOD frameworks.

3Domain discovery is possible via sequence- and structure-based approaches, with Pfam and Gene3D as prime examples of the former. The Pfam database describes collections of protein families, each represented by a multiple sequence alignment (MSA) and a hidden Markov model (HMM) profile. In Gene3D, sequences of existing CATH superfamilies, assigned via structure, are leveraged to discover new domains in sequence space via representative profile HMMs. Sequence-based discovery allows for greater coverage but is limited by HMMs detection capabilities, often failing on remote relatives. Structure-based assignments allow for higher quality domain boundaries and can reveal very remote relatives but have been limited by the low numbers of experimental structures. The task is further complicated by conceptual difficulties when performing structural comparisons, such as the lack of a statistical null model (compared to the use of random sequences for sequence comparisons).

4An essential perspective yet to be explored in the context of the AFDB is the comprehensive mapping and analysis of protein domains across various branches of the Tree of Life using structural data. The relationships between protein folds and domains have been extensively highlighted by structural databases such as CATH, ECOD, SCOP and SCOPe, with differences in criteria for which regions are assigned to domains in proteins. Although early comparisons across domain databases showed agreement for fold level assignments, significant differences exist for the definition of a protein domain. CATH recognises that some structures can be decomposed into further structural units, each with their own repertoire of observed variations, while SCOP takes into account the idea of fold recurrence - whether a subunit has been observed as reoccuring in another family or only as a single domain.

5Examination of the AFDB through the lens of the CATH framework holds the potential to reveal unprecedented insights that can illuminate subtle yet important connections between structure and function across large swathes of organisms. As such, a structure-based mapping of domains within the AFDB promises not only to massively facilitate the exploration of these relationships, but also provide a foundation for further annotation.

6In this study, we present a comprehensive analysis of domain composition for the entirety of the AFDB (version 4). This description encompasses over 364 million putative domains, derived from more than 214 million protein sequences across more than 1 million taxa. The identification of these domain structures is made via the consensus of three automated parsing methodologies: Merizo, Chainsaw and UniDoc. Further, by employing structural comparison methods such as Foldseek and an in-house deep learning method called Merizo-search, more than 251 million domains can be placed on the CATH hierarchy.

7Our workflow combines three state-of-the-art domain parsing methods together with structure classification algorithms to identify known domain folds (Fig. 1a-b; Methods) within the AFDB. Using this workflow, we identified a total of 364 million ‘TED’ domains across the AFDB (Fig. S2, Table S1) - 100 million more domains than found via sequence-based methods (Fig. 1c). TED-100 describes a roughly 42:55 division between single and multidomain proteins (Fig. 1d-i), with the latter up to 20 domains in composition (Fig. S3). Only 2.8% of targets (5.3 million) in TED-100 lack identifiable domains, compared to 33.9% (64.1 million targets) in Gene3D and 26.2% in Pfam (49.4 million targets) (Fig. S4). In TED, these targets either consist entirely of non-domain residues (NDR) or lack any consensus among the three domain segmentation methods employed. Across superkingdoms, the fraction of NDRs varies, with approximately 10% identified across archaeal and bacterial lineages and up to 30% in eukaryotes (Fig. S5).

8Analysis of the average predicted local distance difference test (plDDT) scores for TED-100 reveals that the majority of domains fall within the "very high" to "high" bins, with only approximately 2% within the lowest bin (Fig. 1d-ii). Since our domain segmentation methods do not consider residue plDDT when determining domains, this suggests that our domain finding pipeline effectively identifies plausible domains within the well-folded areas of AF2 models.

9The 324 million domains in TED-100 were clustered by sequence using MMseqs2 and compared against CATH representative domains using fast structure searching methods (Fig. 1a-iii, Methods). This produced approximately 121 million clusters at 50% sequence identity and a minimum coverage of 90% (Fig. S6; Methods). The majority of these clusters comprise just singleton sequences (roughly 81 million), with the largest non-singleton cluster comprising 12,847 domains.

10Parallel to sequence clustering, we use a combination of Foldseek and Merizo-search (an in-house structure search method using domain coordinate embedding) to search all TED-100 domains against CATH SSG5 domains (Methods), allowing 194 million domains to be assigned with CATH superfamily (H) labels and 46 million at the topology (T) level (Methods). These labels were further validated by scanning domain sequences against an updated library of HMMs for CATH PDB domains. Approximately 171 million superfamily predictions by Foldseek could be confirmed with exact HMM superfamily matches (88.54%), with an additional 1.8 million domains (0.95%) confirmable at the fold level. Only 4.1 million of the 16 million Foldseek predictions for fold matches on CATH can be validated by HMM scans (25.8%), with 11.8 million fold predictions and 20.3 million superfamily predictions by Foldseek not confirmed by an HMM match, suggesting an expansion by 15.4% in CATH labelled domain coverage using AFDB structures over HMM-based sequence assignments.

11By identifying sequence clusters with any CATH label assigned domains, the clusters were partitioned into two categories: 78 million clusters (over 251 million domains, nearly 80% of all TED-100 domains) which contain at least one CATH-labelled member, spanning 148 million proteins, and 26 million proteins having no domain annotations in Pfam and 30 million having none in Gene3D (Fig. S7). The remaining 41 million clusters have no members with CATH labels (approximately 73 million domains; Fig. S6). The absence of similarity to CATH domains in the latter clusters could be attributed to them being novel folds, extremely divergent relatives of existing CATH domains or simply being incorrect models, making them unmatchable to known folds.

12To develop an understanding of how folds are distributed across the AFDB, we assessed the composition of TED using the CATH hierarchy. Fig. 2a shows the top 100 CATH superfamilies of each class (alpha, beta and alpha/beta), which are greatly enriched in TED-100 compared to baseline sequence hits in Gene3D.

13These folds and architectures are unevenly distributed across the Tree of Life, however the majority of folds (61%) are shared and reused across all superkingdoms, suggesting essential roles for cellular life. Some folds were found across two superkingdoms (18.5%), while others are more exclusive, with 0.5%, 9% and 11% of CATH folds found only in archaea, eukarya and bacteria, respectively.

14The most abundant superfamilies in TED-100 and CATH are shown in Fig. S8 and Fig. S9. Directly comparing the top superfamilies assigned in TED compared to CATH, in terms of raw domain counts, sees the promotion of several superfamilies into the top 5 of each class, including the MFS general substrate transporter-like domain, translation factors and the FAD/NAD(P)-binding domain (Fig. S8). The set of superfamilies highly enriched in TED include those associated with the archetypal multi-drug efflux pump AcrB. AcrB is a Resistance-Nodulation-Division (RND) transporter and forms part of the AcrAB-TolC efflux pump in bacteria where it is responsible for the export of harmful substances such as antibiotics, contributing to antibiotic resistance. The constituent domains of AcrB, including the pore domain, transmembrane domain and TolC docking domain, are found greatly enriched in TED compared to the PDB, providing up to nearly a thousand-fold increase in representation for these biologically important superfamilies.

15The pore domain, which forms part of the pore selectivity filter in RND transporters, is principally found in bacterial species and in a small minority of archaeal and eukaryotic lineages based on Gene3D hits. However, structure-based searching via TED expands the coverage of the pore domain superfamily into an additional 18 archaeal, 1315 bacterial and 284 eukaryotic lineages unique to TED, which evaded even HMM searches. This broader coverage of organisms, revealed by structural comparisons, may reveal potential evolutionary events, such as lateral gene transfer between bacteria and eukaryotes. This expansion is exemplified in Fig. 2c, where we show that over 1000 CATH superfamilies, described by Gene3D sequence matches as occurring solely in a single superkingdom, could actually be mapped to other branches when structure was taken into account in TED.

16TED also enables us to better study the distribution of folds across different branches of the tree of life. Among the 193 million TED domains with superfamily labels, we observed folds that were exclusively localised to specific superkingdoms. These findings are shown in Fig. 2d, where folds are summarised at the CATH topology level, visualised by principal component analysis (PCA), and where points at the vertices represent exclusive occurrences of a topology and all its superfamilies, in a particular superkingdom.

17From our TED workflow (Fig. 1), we identified 41 million sequence clusters which could not be linked to CATH superfamilies. Representatives of these clusters were subjected to a workflow aimed at identifying novel domain folds (Fig. 1a-iv; Methods).

18While reviewing these clusters, it was apparent that we would need to treat repeat architectures with high internal symmetry separately. A good example of domains in this class are the various WD40 beta propellers, which are considered distinct domain architectures in their own right, but clearly comprise repeats of domain-like units. To identify similar domains in our workflow, we calculated Z-scores using the SymD program, sequestering any cluster representatives with a symmetry Z-score of greater than 9, into a new category of 6,433 highly symmetric novel fold clusters (Fig. S10).

19Within these clusters, we find new architectures such as an 11-bladed beta-propeller, a closed alpha ring-like 11-helix propeller and 6-helix propeller which have not been seen before (Fig. 3). Various other propeller arrangements including the “beta-flower” domain shown by Durairaj et al., were also identified, with some visually striking examples shown in Fig. 3.

20More curiously, we find a broad category of architectures which are composed of cyclic repeats, extruded along an axis to form highly repetitive and symmetric structures which we call “extruded repeats”. Several studies and databases such as the Database of Structural Repeats in Proteins (DbStRiPs) have curated repeat architectures from the PDB, including alpha- and beta-solenoids, and horseshoes formed of alpha-repeats such as HEAT, ankyrin and armadillo. Fig. 3 showcases some complex examples of these structures, with some featuring highly varied and unstructured loops between repeating units (Fig. S11). Many of these domains architectures resemble other solenoid folds that were recently reported in a systematic study of beta-solenoid fold space.

21The remaining low-symmetry clusters were assessed on domain quality, using a variant of the Foldclass network trained to identify poor quality domain choppings (Supp. Methods), and novelty by further matching against known structure libraries. The final output of our workflow produced 7,427 clusters of putative domains which appear to be well-folded but dissimilar to any known domain fold. Although there is no exact boundary between novel and just highly divergent examples of known folds, by applying a density-based anomaly detection algorithm (Methods), we could at least rank these domains by novelty relative to known domains.

22Fig. 4 showcases examples of novel domains identified by our workflow. Over a quarter of clusters corresponded to singletons at the sequence cluster level (1930 domains), which are uniformly distributed in terms of the Foldclass novelty score. Most novel domains identified are from bacterial proteins, which is unsurprising given that the latter comprise most proteins in the AFDB and in our TED-100 dataset. Overall, these domains are distributed evenly across different phyla, but several bacterial phyla were found overrepresented when compared to baseline counts across all domains in TED-100, specifically in the PVC group, Myxococcota, Spirochaetota, Bdellovibrionota, Nitrospinae/Tectomicrobia group and Calditrichota (Fig. S12), suggesting that species within these phyla may be underrepresented in terms of domain coverage in the current PDB.

23Ranked first by novelty, is a curious archaeal domain found as a sequence singleton from Candiatus Poseidoniales archaeon (TED: A0A7M3WA57_TED05; novelty score: 77.2; Fig. 4c-i). This protein is not documented in InterPro and no GO terms are available. The structure is composed of paired beta-strands, in a closed, twisted hairpin with both termini adjacent to one another. Reviewing the context of the full-chain from the AFDB shows that the domain forms part of an extended loop, protruding from the middle of a immunoglobulin-like domain (TED: A0A7M3WA57_TED04). The hairpin topology of A0A7M3WA57_TED05 is mirrored by another identified novel domain represented by E1Z635_TED02 shown in Fig. 4d-ii, but is alpha-helical in nature. This domain is found only in eukaryotes, primarily in species belonging to the Viridiplantae phyla (97 species) but also in a minority of Opisthokonta (14 species) and Amoebozoa (1 species), suggesting an evolutionary link between these lineages.

24To see if any functions could be assigned to the domains with novel folds and repeats, we use a sequence-based deep learning model (Supp. Methods) to predict Gene Ontology (GO) terms. This analysis shows that 1321/7427 (18%) of the domains in the putative novel fold set, and 1419/6433 (22%) of the repeat set can be assigned high confidence (p < 10-4) Molecular Function GO term labels. The top 20 GO terms predicted for the two sets of domains are shown in Table S4 and Table S5. Manual inspection of the domains predicted to have zinc binding and nucleic acid binding functions reveals that many of the domains contain plausible zinc binding sites, most containing 2 Cys and 2 His residues arranged in tetrahedral fashion, including as part of zinc finger-like supersecondary structure motifs. Fig. 4e-i shows one example containing two zinc binding sites and a possible nucleic acid binding alpha-helix, but which lacks the canonical zinc finger supersecondary motif. We also consider the set of domains predicted to have heme binding properties and find that most of these contain the canonical heme c binding motif (CXXCH). Inspection of the three-dimensional structures reveals that each of these domains has one or more heme binding sites in a plausible conformation; one example is shown in Fig. 4e-ii, with the residues of the heme c binding motif highlighted. The His residue that binds the heme iron is in a conformation compatible with placement of the heme group in the pocket, which is primarily hydrophobic. The presence of clear sequence motifs and structural features consistent with the assigned functions suggests that these novel domains may indeed have the predicted functions, and further work is needed to validate the remainder of the predicted functions for these novel domains.

25Unlike sequence-based domain annotations such as Gene3D, the availability of full-chain AF2 models for multidomain proteins allows us the unique opportunity to interrogate and compare domain pair packing interactions in TED and CATH. TED contains a total of 27,280,057 instances of interacting domains categorised into 13,771 Interacting Superfamily Pairs (ISPs). In contrast, the set of interacting domains in CATH consists of just 196,234 instances across 5,111 ISPs. The relative enrichment in the number of instances of ISPs common to CATH and TED is shown in Fig. 5a, which shows that most of these ISPs have many more members in TED. Assessing the diversity of interaction geometry for ISPs common to TED and CATH (Methods and Fig. 5b-i-ii), we find that for most ISPs, the diversity of interaction geometries in TED is consistent with that seen in CATH, indicating that on average, AF2 tends to recapitulate inter-domain geometries already seen in the PDB.

26A small proportion (5.4%) of these ISPs show enhanced diversity in interaction geometries in TED, as measured by an increase in Conservation of Interaction Orientation (CIO) score (Methods) of 0.3 or more. A smaller proportion of ISPs (2.3%) are more diverse in CATH. We find no strong dependence of difference in CIO value on whether the ISP in question is homotypic (2 copies of domains from the same superfamily) or heterotypic (Welch two-sample t-test, p=0.141). Most of the interacting superfamily pairs (ISPs) in the TED set (10,701 out of 13,771) are unique to TED, i.e. they are not observed in CATH, and CATH contains 2,041 ISPs not seen in TED (Fig. S13).

27That some interactions in CATH are not captured in TED is not entirely surprising, firstly due to the possibility of domain parsing and classification errors. Secondly, TED excludes any sequence not modelled in the AFDB; these include sequences from viruses (around 5,304,757 sequences), and longer sequences, which we estimate contain between 50 and 61 million domains (Supplementary Methods). The AFDB will also not contain experimental constructs comprising domain combinations not observed in natural sequences. Lastly, instances of ISPs in TED are filtered based on whether the constituent domains are in contact in the full-length AF2 structure, and whether they have a favourable inter-domain Predicted Aligned Error (PAE; Methods), both of which depend on the availability of homologous sequences and template structures that ideally span both domains. For these reasons, we do not expect that TED ISPs will be a strict superset of CATH ones. Nevertheless, we find that TED essentially doubles the set of known domain interactions at the superfamily level (10,701 novel TED interactions versus 5,111 known interactions in CATH).

28A visual illustration of ISP sets can be seen in Fig. 5c-i, in which a path is drawn between two superfamilies if at least one interaction is observed between domains assigned to those superfamilies. This representation uses the CATH hierarchy as a guide to ‘bundle’ paths drawn between superfamilies in related parts of the CATH hierarchy (Methods; Fig. S13), and shows visually that a very large number of new interactions are seen, especially between superfamilies in CATH classes 2 and 3 (all-beta and alpha-beta classes, respectively). As mentioned above, the majority of the interactions seen in the TED set are unique to TED, and a visual comparison of the set of TED-unique interactions to all known interactions in CATH (Fig. S13) reveals a huge expansion in the set of known interactions. The network of superfamily interactions also allows us to identify superfamilies that can be considered hubs on account of their ability to interact with many other superfamilies. As shown in Fig. 5c-i-ii, a large number of superfamilies in CATH have their hub status ‘promoted’ due to new interactions seen in TED. In Fig. 5d-i-ii we show two examples of superfamilies that have been promoted in this way. Notably, the set of interactions observed in TED for these superfamilies show that they can contain only novel (in the case of 3.40.50.10850) or, in some cases, a mixture of novel and already-observed interaction modes (for 2.40.30.200).

29The interaction data shows that TED greatly expands the set of known interactions and interaction modes between domains. Future work will aim to delve deeper into the interaction data and determine the functional roles that these interactions (particularly those not seen before) might play in cellular processes, in particular, taking into account the multi-domain architecture of individual proteins and their evolutionary relationships.

30Of the 214 million structures in the AFDB, nearly 39 million are exact sequence duplicates of other proteins in the database (13 million unique sequences within this set). The 26 million redundant proteins comprise our TED-redundant set (Methods).

31Remarkably, the AFDB models for these sequence-redundant proteins often diverge from one another, with approximately 42% of clusters (5.6 million) showing a maximum cluster RMSD of greater than 1Å (Fig. S16). The very largest RMSDs in the distribution tend to relate to changes in domain packing, but even at the domain fold level, changes can be observed. Fig. S16 shows the distribution of maximum pairwise RMSD for each cluster of identical sequences. We found structural variation at the chain-level (Fig. S16a), as well as in the PAE maps generated by AF2 (Fig. S16a-ii).

32One explanation for our observations here is that we are looking at alternative conformers of the protein chains relating to different MSAs. However, in these cases, the input MSAs should be identical given the described modelling protocol. The larger differences are also far too large to attribute to the short relaxation step that follows the model prediction step. This is most evident in the example shown in Fig. S16a-ii, where two AFDB models for identical sequences deviate by nearly 65Å and show clear differences in the PAE map.

33To further investigate structural diversity within these sequence-redundant clusters, we subjected the TED-redundant targets to our domain parsing workflow, identifying approximately 40 million domains across the set. In addition, we also derived a consensus across all structures in each identical sequence cluster (Methods). This allowed us to investigate domain-level changes in conformation between identical sequences, and identified many cases where the consensus domains were dramatically different (Fig. S16b).

34What we have shown so far in developing TED is a way in which structural data in the AFDB can be augmented, by carefully breaking down structures into their component domains, allowing them to be classified through the CATH framework. This initiative not only drives forward the associations that we can make between structure and function, but as shown in our study, can be used to discover and reclaim the dark areas of fold space that are not accessible to sequence-based discovery.

35Comparing TED with a recent study on the 21 model organisms dataset (released prior to the 214 million release), already shows that the TED workflow identifies not only a greater number of domains, but that the domains are also of higher quality and capture many more remote homologies (Table S3). TED currently annotates domains for over 1 million taxa, of which 600,000 are currently mapped to CATH domains in TED-100, and such an extended mapping of domains to superfamilies will enable many more evolutionary discoveries. A good example of how such expansion of CATH superfamilies enhances understanding of evolutionary processes and aids inheritance of functional properties, even towards drug repurposing is shown in Fig. S17.

36The coverage of TED dwarfs sequence-based assignment methods (namely Gene3D and Pfam), identifying over 100 million more domains in the AFDB compared to the latter. The proportion of proteins that each method can find domains within is shown in Fig. S4, which illustrates the advantage of considering structural domains (TED finds domains in 40-50 million proteins that Gene3D and Pfam cannot). Interestingly, prior studies based purely on sequence comparisons have suggested that between 40-65% of prokaryotic proteins are composed of multiple domains, with a higher proportion proposed in eukaryotes.

37These values are comparable to those seen in TED, where we find a roughly 42:55% split between single and multidomain proteins in the AFDB (Fig. S4). Compared to TED, the proportion of multidomain proteins is much lower in Gene3D (29%) and Pfam (24%) assignments (Fig. S4).

38Although most structures in the AFDB are undeniably of high quality, the sheer scale of the data means that errors and anomalies are inevitable, and these should be discussed in order to assess the overall robustness of the data. Some limitations have already been pointed out by other studies. One such idiosyncrasy that appeared during our development of TED was the observation that models of 100% identical sequences were sometimes dramatically different. To reframe the issue, the implication of this is that for a given protein, a user may be able to find a better model or domain within the AFDB (and one that AF2 clearly places higher confidence in) if sequence-redundant copies are considered. This may mean that some alternative structures could only be detected when redundant copies are available to compare against. Given that the vast majority of AFDB entries (175/214 million) do not have duplicates, the implication here is that an unknown (but probably large) number of low-quality structures or regions in AFDB might be improvable by resampling the input MSAs and reassessing the quality of the models. Even though this might be too heavy a task to do across the whole of AFDB, users of the database should certainly bear this in mind when looking at specific domains of interest.

39The obvious explanation for the existence of divergent models for a given sequence must be that the models were generated using significantly different MSA information, and we note that explicit pre-sampling of the MSAs given to AF2 has been explored recently as a way to persuade the network to generate alternate conformations. As the structure generation pipeline in AF2 includes possible random MSA resampling during the forward pass, it could be that a borderline MSA could result in very divergent predictions on different runs due to this resampling (fixed random seeds were not used in making the AFDB). Another possibility is that the sequence data banks used for different instances of the same target sequence were inadvertently updated to newer versions if they were predicted at later times. We were not able to investigate these possibilities further ourselves as there is no information on the MSAs provided in the AFDB data, but we do suggest that it is a topic worthy of further investigation by the AFDB developers.

40Fortunately, we find many examples that reiterate the notion that plDDT is an invaluable discriminator of the confidence that a user should place on a model. Fig. S18 shows an example whereby an AFDB entry has been generated within the very low confidence bin. This model features a significant number of structural defects not typical of AF2 models. “Re-folding” the sequence in ColabFold produces a number of visually striking structures across the five AF2 models, the variation of which, along with their unusually low plDDT scores, strongly indicates that AF2 has hallucinated these folds and that they should likely be disregarded.

41However, overall, the large proportion of domains mapping to CATH evolutionary families gives confidence in the quality of the models, capturing well the structural features of folds and preserving their distinctive structural characteristics. In this context, the performance of the domain segmentation algorithms deployed here has also been key, as early pilot work on the 21 model organisms suggested a much higher proportion of problematic AF2 models largely caused by the poor segmentation of full-length proteins into domains by the sequence-based methods. As well as confirming earlier hypotheses that the majority of domain structural families had already been characterised experimentally, our study has also revealed some intriguing and beautiful new domain architectures and folds, especially some highly symmetric repetitive structures. Future development of TED will aim to comprehensively analyse these new folds and incorporate them into the CATH hierarchy, which will require extensive manual curation. Additionally, we recognize the need to optimise the current TED workflow for better detection of repeat structures, such as the different extruded repeats highlighted in our work. This optimisation will necessitate the development and integration of new tools and analysis methods specifically designed to identify these types of structures.

42Throughout our study, we had to make a number of algorithmic decisions which were primarily motivated by the monumental amount of data we had to process in the AFDB. As such, we intend for TED to be an ongoing development, which will evolve as the data and the needs of its users do. Our aim is to provide the community with the most comprehensive summary and breakdown of the structures within the AFDB. We expect TED to be used as a starting point for a whole host of analyses, including providing a comprehensive dataset to train and test a new generation of deep learning based applications in structural biology.