Geometric description of self-interaction potential in symmetric protein complexes

Proteins can self-associate with copies of themselves to form symmetric complexes called homomers. Homomers are widespread in all kingdoms of life and allow for unique geometric and functional properties, as reflected in viral capsids or allostery. Once a protein forms a homomer, however, its internal symmetry can compound the effect of point mutations and trigger uncontrolled self-assembly into high-order structures. We identified mutation hot spots for supramolecular assembly, which are predictable by geometry. Here, we present a dataset of descriptors that characterize these hot spot positions both geometrically and chemically, as well as computer scripts allowing the calculation and visualization of these properties for homomers of choice. Since the biological relevance of homomers is not readily available from their X-ray crystallographic structure, we also provide reliability estimates obtained by methods we recently developed. These data have implications in the study of disease-causing mutations, protein evolution and can be exploited in the design of biomaterials.


Background & Summary
The controlled association of proteins into functional complexes is central to the myriad of biochemical processes required to maintain cellular functions 1,2 . The symmetry of protein complexes enables unique biological properties: compact genetic encoding of large assemblies such as viral capsids, cytoskeleton tubules and filaments, or cooperative, switch-like transitions involving allostery. However, we recently observed that the repetition of subunits within homomers can exacerbate the effect of point mutations, resulting in the homomer's uncontrolled self-assembly 3 .
For a new mode of protein assembly to take place, a new interaction must be created. Previous work showed that the chemical composition of protein interfaces, although distinct from surfaces, is relatively close. Indeed, two amino-acid substitutions are sufficient, on average, to shift the chemical composition of a protein surface patch into that of an interface 4 , suggesting that point mutations may frequently trigger new interactions, as in the sickle-cell disease 5 .
Here, we need to distinguish homotypic interactions, where two identical parts of the structure are in contact, from heterotypic interaction, where two distinct structural parts are in contact. Homotypic interactions are more frequently sampled by chance than heterotypic interactions are [6][7][8] . When occurring at the surface of a monomer or at the surface of a cyclic complex, a new homotypic interaction will likely yield a finite dimerization event [9][10][11] . However, among homomers with dihedral symmetry, the emergence of a new self-interaction necessarily triggers an infinite (open) self-assembly 9 .
In our previous work 3 , we introduced point mutations solely designed to increase surface hydrophobicity into 12 dihedral homomers from Escherichia coli. Remarkably, these mutations triggered new self-interactions resulting in all complexes forming high-order supramolecular assemblies both in vitro and in vivo upon heterologous expression in Saccharomyces cerevisiae. Structural examination of these mutants allowed us to identify a novel descriptor: the normal distance to the closest bounding plane (nDp) of a symmetric oligomer, which describes a residue's position on the global quaternary structure. The lower the nDp, the closer the amino acid is to the apex 2 Scientific Data | (2019) 6:64 | https://doi.org/10.1038/s41597-019-0058-x www.nature.com/scientificdata www.nature.com/scientificdata/ or "tip" of the assembly along a symmetry axis, and the more its mutation has the potential to trigger the formation of a high-order assembly 3 .
Accordingly, we then showed that the greater potential of these geometric hot spots to trigger assemblies was counterbalanced chemically by an enrichment in hydrophilic amino acids 3 . We measured the interaction propensity of surface regions on 1,990 dihedral homomers of known structure using the 'stickiness' scale introduced by Levy et al. 3,12 and detailed below. Our results indicated that surface regions with high potential to trigger supramolecular assemblies upon mutation (i.e. low nDp) counterbalanced this risk by residues with low interaction propensity, or stickiness 3 .
Here, we present a dataset of descriptors that characterize these geometric hot spot positions and buffering effects on 165,916 proposed biological assemblies from the Protein Data Bank (PDB) 13,14 , together with the workflow and computer scripts used to compute these descriptors 15 .
These data serve multiple uses: (i) they will be important to consider in future studies predicting the molecular consequences of mutations, including single nucleotide polymorphisms, (ii) from an evolutionary standpoint, they describe molecular phenotypes that may constrain amino acid changes and thereby, could be considered in phylogenetic models of sequence evolution, and (iii) in the field of bio-materials design, these data facilitate the application of our simple strategy to program protein self-assembly at length scales up to several micrometers either in vitro or in vivo, using the PDB as a source of natural "building blocks".

Methods
Normal distance to the closest bounding plane (nDp) calculation. To study the effects of point mutations on symmetric homomers, we defined a novel structural descriptor based on quaternary structure geometry. We called this descriptor the "normal distance to the closest bounding plane", or nDp. These methods are expanded versions of descriptions in our related work 3 .
We reasoned that for point mutations to act synergistically in the creation of novel self-interacting interfaces, the affected residues at the surface of one copy of the homomer must be altogether accessible to the surface of other copies of this oligomer. Bounding planes, which are orthogonal to symmetry axes, capture such information. The nDp measure thus describes the distance of a residue from the closest apex of a quaternary structure along a symmetry axis (Fig. 1). The lower a residue's nDp, the higher its potential to mediate interactions with another copy of the homomer, and the more likely it is to trigger a novel self-interacting interface upon mutation 3 . www.nature.com/scientificdata www.nature.com/scientificdata/ To calculate the nDp, a symmetry axis is considered as a unit (1 Å) vector s originating from the center of mass of the assembly. Similarly, the Cα of each residue i defines a vector r i originating from the center of mass. For each symmetry axis a of the assembly, two bounding planes parallel to one another are defined. They are orthogonal to the symmetry axis considered, and intersect at the maximal (d a,max ) and minimal (d a,min ) values of the dot product s · r i , considering all residues i of the quaternary structure. The measure nDp for a given residue i is calculated with respect to a particular axis as the minimal distance to either of its bounding planes a, as follows: nDp a,i = min (d a,max − s · r i , s · r i − d a,min ) 7 .
Among cyclic complexes, which have a single axis of symmetry, there is no ambiguity to calculate nDp with the formula above. However, homomers with dihedral symmetry have multiple axes of symmetry, so multiple nDp values can be computed for each residue (one for each symmetry axis). Here, we consider three cases: (i) nDp relative to bounding planes originating from 2-fold axes, where each residue is assigned the lowest nDp value relative to all 2-fold axes (i.e. nDp-low-fold or nDp-2-fold, Fig. 1c), (ii) nDp relative to bounding planes originating from the n-fold axis (i.e. nDp-high-fold or nDp-n-fold, Fig. 1d), and (iii) nDp relative to all bounding planes originating from all axes, whereby each residue is assigned the lowest nDp value relative to all axes (i.e. nDp, Fig. 1e). In our previous study 3 , we employed this definition.
Importantly, D2 homomers have three 2-fold axes and so it is not possible to distinguish between axes' folds. Thus, for those we only employ nDp definition number 3. environment stickiness calculation. In our previous work, we observed that regions with high geometric potential to trigger self-assembly counterbalanced that potential by negative design consisting of a lower than average chemical potential for self-assembly. We measured the chemical potential for self-assembly of a given surface patch by the "stickiness" of amino acids it contains, introduced in our previous work 12 and described in detail below.
The stickiness of an amino acid is defined as the log-ratio of its frequency at protein-protein interfaces relative to solvent-exposed surfaces (Fig. 2a). The stickiness scale thus quantifies the trade-off between the probabilities of finding a given amino acid involved in an interaction with another protein versus being in a solvated environment (Fig. 2a) 12 . Its calculation is based on a set of 397 non-redundant protein structures from E. coli. Surface and interface protein regions were defined using the residues relative accessible solvent area in the complexed and unbound www.nature.com/scientificdata www.nature.com/scientificdata/ states (rASAc and rASAu, respectively) 4,12 . If a residue has a rASAc value superior to 25% and the delta between rASAc and rASAu is null, then this residue is assigned to the surface (ΔrASA = 0 & rASAc > 25%). Interface residues were defined as those belonging to the interface core (ΔrASA > 0 & rASAc < 25% & rASAu > 25%). The stickiness scale employed here is based on E. coli proteins, but it is robust to using different sets of proteins. For example, deriving stickiness scales based on proteins from S. cerevisiae and H. sapiens showed high correlation values (R E. coli−S. cere = 0.94, R E. coli−S. sapi = 0.97) 12 .
Next, the 'environment stickiness' of a residue of interest is calculated based on its surrounding surface residues, by averaging their stickiness values (Fig. 2b) 12 . The residue at the center of the patch is excluded since we focus on quantifying the buffering effects in the residue's vicinity. The reasoning behind this approach is that residues in more sticky environments are expected to have a higher probability of triggering protein-protein interfaces upon mutation to more sticky or more hydrophobic residues 12 . Surrounding surface residues are defined as those whose Cα is located within a 400 Å 2 patch centered on the Cα of the residue of interest (i.e. a maximum Cα-Cα distance of 11.28 Å). The surface region defined for the environment stickiness calculation are associated to a rASAc > 25%, without considering any delta between rASAc and rASAu. All buried residues (rASAc < 25%) are ignored and no stickiness is computed for those.
Biological relevance of homomers. The biologically relevant quaternary structure (QS) of a protein is not readily available from its X-ray crystallographic structure, which provides the atomic coordinates of the asymmetric unit (ASU) only. Indeed, the QS may be formed by parts of several ASUs or be a sub-part of one ASU. The challenge is, therefore, to distinguish fortuitous crystal contacts from biological ones forming the QS 16,17 . Numerous approaches such as PISA 18 and EPPIC 19 have been developed to predict QS information from X-ray crystallographic structures. In this dataset we provide predictions based on the integration of PISA and EPPIC approaches together with novel ones we recently developed, named QSalign/anti-QSalign and QSbio 20 . These methods are summarized from descriptions in our related work 20 .
QSalign employs evolutionary conservation of quaternary structure geometry as evidence of biological significance 20 . Quaternary structure conservation is inferred following the structural superposition of full homomers using Kpax 21 and is quantified by a multichain version of the TM-score 22 . Anti-QSalign takes a complementary approach where the absence of QS of homologues is predictive of a monomeric state.
Lastly, QSbio scores the relevance of a QS based on the predictions from three methods (PISA 18 , EPPIC 19 , QSalign/anti-QSalign 20 ) and provides a confidence estimate per assembly in the form of a probability for the QS to be incorrect 20 . Those probabilities are estimated based on a benchmark (Fig. 3), and are given in the table of assemblies descriptors (protein_assemblies_description.csv.tar.gz 15 ).
Other descriptors acquisition. Assemblies descriptors were retrieved from the 3DComplex database 23 : number of subunits, molecular weight, resolution, symmetry types, symmetry axes and Uniprot 24 accession codes (protein_assemblies_description.csv.tar.gz 15 ). Regarding residue descriptors, absolute and relative accessible surface area (ASA) were calculated using CCP4 25 Areaimol 26,27 . Relative ASA values initially superior to 100 were corrected to 100. For convenience, stickiness scale values from Levy et al. 12 were also included for each residue entry.
Datasets construction. As a starting point to build the datasets of assemblies we present in this paper, we interrogated the 3DComplex database 23 to retrieve assemblies that: (i) do not break into separated sub-structures when ignoring subunit-subunit contacts of less than 5 residues per chain on average, (ii) have at least one domain defined in either SCOP 28 , Pfam 28,29 or ECOD 30 , (iii) do not contain superposed chains, and (iv) do not exclusively contain Cα information (low resolution structures). This process allowed us to retrieve 165,916 proposed biological assemblies from the PDB 13 for which all descriptors cited in this study are provided 15 .  (Table 1). All data and scripts are available on figshare at: https://doi.org/10.6084/m9.figshare.6586958.v2 15 . We provide residue descriptors as 4 tables that regroup either symmetric or asymmetric protein structures ('sym' and 'asym' indicators, respectively) for either all 165,916 assemblies retrieved from the PDB 13 or their non-redundant subset of 40,109 assemblies ('all' and 'h80' indicators, respectively). To facilitate data loading and manipulation, table 'residues_h80_sym_ protein_assemblies' is a non-redundant subset of table 'residues_all_sym_protein_assemblies' and table 'resi-dues_h80_asym_protein_assemblies' is a non-redundant subset of table 'residues_all_asym_protein_assemblies' . Non-redundant subsets were derived from 3DComplex 23 . This process eliminates proteins that share the same domain architecture as defined in SCOP 28 , Pfam 29 or ECOD 30 and more than 80% sequence identity. Importantly, the quaternary structure is taken into account when filtering redundant structures, so different quaternary structures sharing the same sequence are kept. assembly descriptors. Table 'protein_assemblies_description' stores general assemblies descriptors (Table 2). Each line corresponds to one of the 165,916 assemblies of the complete dataset. Fields 'h_80' , 'h_90' and 'tv_sticky_discard' are binary values (0/1) indicating, respectively, whether the assembly belongs to a non-redundant subset using either a 80% sequence identity threshold, a 90% sequence identity threshold, and whether it was ignored to perform technical validation (see section "Technical Validation", Fig. 4c).
assembly symmetry axes coordinates. Information on symmetry axes coordinates is stored in table protein_assemblies_symmetry_axes.csv.tar.gz 15 (Table 3). Each line corresponds to one symmetry axis that belongs to one of the 69,922 symmetric assemblies in the datasets, minus 731 assemblies for which symmetry axes were incorrect. Tables 'residues_all_sym_protein_assemblies' , 'residues_all_asym_protein_assem-blies' , 'residues_h80_sym_protein_assemblies' and 'residues_h80_asym_protein_assemblies' store residue descriptors (Table 4) 15 . Each line corresponds to one unique residue of a structure's assembly. Please note that  www.nature.com/scientificdata www.nature.com/scientificdata/ descriptor 'nDp' is defined for monomers in tables 'residues_all_asym_protein_assemblies' and 'residues_h80_ asym_protein_assemblies' only because we use it as a control after having generated one random symmetry axis per monomer (see section "Technical Validation"). Otherwise, calculating any nDp version on asymmetrical structures is irrelevant because it directly depends on the symmetry axes of an assembly.

technical Validation
Both manual inspection of individual examples as well as global analyses were performed to ensure the validity of the data. All residue descriptors were projected onto a few hundred assemblies for visual inspection using PyMol 31 . Measurement tools then allowed for manual validation of different nDp versions and environment stickiness calculations on several randomly selected residues. While recalculating data to provide updated datasets, we added a negative control by generating hypothetical nDp values on monomers. We used a single randomly oriented axis passing through the centroid of each monomer structure, and calculated nDp values as for any self-assembling high-order structure using this single axis (see Methods). This technical validation was performed using the non-redundant subsets of assemblies exclusively: tables 'residues_h80_sym_protein_assemblies' and 'residues_h80_asym_protein_assemblies' and ignoring assemblies for which the QS error probability calculated by QSbio was very high, i.e. above 50%.
As expected, the average and maximum nDp per assembly increase linearly with assemblies' molecular weight, regardless of symmetry types and including control (Fig. 4a). Since all structures are considered in terms of biological assemblies, the average and maximum nDp per assembly within a given range of molecular weight   www.nature.com/scientificdata www.nature.com/scientificdata/ is superior for control and cyclic complexes compared to dihedral complexes (Fig. 4a). This is due to dihedral complexes having at least 3 orthogonal symmetry axes (i.e. at least 6 bounding planes), whereas control and cyclic complexes only have 1 symmetry axis (i.e. 2 bounding planes). Considering the different nDp definitions at residue level, distributions are very similar regardless of protein symmetry types (Fig. 4b). Only the distributions of the nDp-2-fold and nDp-n-fold, which are exclusively calculated on high-order dihedral complexes (Dn, n > 2), are slightly more spread due to the high molecular weights of these assemblies (Fig. 4b). The distribution of the nDp-n-fold is also shifted towards higher values, since high-order dihedral complexes tend to be wider along their 2-fold axes, as for isoaspartyl dipeptidase (Fig. 1b).
Finally, we observed residues' environment stickiness as a function of nDp across different symmetry types (Fig. 4c) and validated our previous computational results: surface regions with high potential to trigger supramolecular assemblies upon mutation (i.e. low nDp) counterbalance this risk by residues with low interaction propensity (i.e. stickiness) 7 . Environment stickiness is tuned as a function of nDp in dihedral complexes, but not in cyclic complexes nor in control (Fig. 4c). To avoid biases due to membrane and viral proteins when analysing surface stickiness, we discarded all assemblies containing one of the following chains of characters in their title, description or function PDB fields: 'lipid' , 'transport' , 'rhodopsin' , 'membran' , 'virus' , 'viral' . We also excluded from the technical validation those assemblies with a high probability to be non-biological (>50%).

code availability
Perl scripts we used to calculate residues' environment stickiness and different nDp versions are provided (Scripts. tar.gz 15 ). Accordingly, the scripts archive contains two folders with demonstration input and output files for each script. A wrapper allows to run all calculations from a PDB file. The PyMol 31 script we used to visualize results on protein structures is also provided. All scripts were extensively commented and made easily readable to facilitate re-use and adaptation (Table 5).

Field
Description Type   Table 4. Residue descriptors records. Each line of tables 'residues_all_sym_protein_assemblies' , 'residues_ all_asym_protein_assemblies' , 'residues_h80_sym_protein_assemblies' and 'residues_h80_asym_protein_ assemblies' 15 corresponds to one unique residue. (*) Descriptor defined for monomers in tables 'residues_all_ asym_protein_assemblies' and 'residues_h80_asym_protein_assemblies' only because we use it as a control (see section "Technical Validation"). (**) Descriptors exclusively related to high-order dihedral complexes (Dn, n > 2) and present only in tables residues_all_sym_protein_assemblies' and 'residues_h80_sym_protein_ assemblies' .  Table 5. Overview of the scripts archive content. Demonstration input and output files are provided for each script.