Site-saturation mutagenesis of 500 human protein domains

1The human genome encodes more than 20,000 proteins. Missense variants in nearly 5,000 of these proteins cause Mendelian diseases, however the functional consequences of nearly all missense variants in nearly all proteins are unknown–. Given the current size of the human population, most variants compatible with life are present in someone currently alive–, making the large-scale experimental analysis of variant function a central challenge for human genetics–. However, experiments have so far mostly quantified variant effects in one or a few proteins,. Despite recent improvements–, computational variant effect predictors (VEPs) are not deemed to provide sufficient evidence to classify clinical variants as pathogenic or benign. They also do not identify the molecular mechanisms by which variants cause disease, information that is important for therapy development and clinical trial design. Whereas many disease variants are likely to destabilize proteins and reduce their abundance,,, others may affect specific molecular interactions or cause gain-of-function phenotypes,.

2Previous studies have established reduced abundance as a frequent causal mechanism for pathogenic variants in diverse proteins,–, but larger-scale studies of human disease variants across many disease genes to test the generality of these observations are lacking. Most human proteins contain multiple independently folding structural units called domains,. For example, the human genome encodes more than 200 homeodomains that bind DNA to control gene expression and more than 250 PDZ domains that mediate protein–protein interactions,. The small size of protein domains (median around 100 amino acids) and their independent folding make them a useful target for large-scale experimental measurement of variant effects.

3Here, using a highly validated assay that quantifies the effects of variants on protein abundance in cells, we perform large-scale mutagenesis of human protein domains. We report the effect of more than 500,000 missense variants on the stability of more than 500 different human domains. This dataset, ‘Human Domainome 1’, provides a large reference dataset for the interpretation of clinical genetic variants and for benchmarking and training computational methods for prediction of variant effects on stability. We use the dataset to quantify the contribution of stability changes to human genetic disease and how this varies across proteins and diseases. We also show how stability measurements can be combined with protein language models to annotate functional sites across proteins, and that measurements made on a small number of proteins can be used to accurately predict stability changes across entire protein families.

4We used microchip-based massive in parallel synthesis (mMPS) technology to construct a library of 1,230,584 amino acid variants in 1,248 structurally diverse protein domains (Supplementary Table 1). In this Human Domainome 1 library, every amino acid is mutated to all other 19 amino acids at every position in each domain (Fig. 1a). Sequencing the library revealed it to be of high quality, with coverage of 91% of designed amino acid substitutions (Supplementary Table 1, Extended Data Fig. 1a). To quantify the effect of these variants on domain stability, we used an abundance protein fragment complementation assay (aPCA),, an in-cell selection system. In aPCA, the protein domain of interest is expressed as a fusion with a fragment of an essential enzyme, and the concentration of this enzyme linearly determines the cellular growth rate over at least three orders of magnitude (Fig. 1b). The effects of variants on protein abundance are quantified using high-throughput sequencing to measure the change in variant frequencies between the input and output cell populations in selection experiments (Fig. 1a). This strategy thus enables pooled cloning, transformation and selection of hundreds of thousands of variants in diverse proteins in a single experiment (Fig. 1a).

5In total, we performed 27 transformation, selection and sequencing experiments (3 independent biological replicates of 9 sub-libraries). After filtering (Methods), the final dataset consists of cellular abundance measurements for 563,534 variants in 522 protein domains, of which 503 are from human proteins (Supplementary Table 2). Abundance measurements were highly reproducible (median Pearson’s correlation coefficient, r = 0.85 between replicates for all variants; Fig. 1c and Extended Data Fig. 1b). They also correlated well with independent in vitro measurements of protein fold stability, (median Spearman’s ρ = 0.73 with folding free-energy changes (∆∆G), n = 10 domains; Fig. 1g,h and Extended Data Fig. 1d,e). Moreover, they correlated well with high-throughput stability measurements from protease sensitivity assays (median ρ = 0.65, n = 13 domains; Fig. 1h and Extended Data Fig. 1f).

6The 522 domains are structurally diverse, comprising 195 in the all-alpha structural class, 127 all-beta domains, 48 mixed alpha and beta domains and 148 metal-binding zinc-finger domains (Fig. 1e). In total, they cover 127 different domain families, including 14 families with 10 or more homologous domains (Figs. 1f and 2), and 97 families with only 1 or 2 domains (Fig. 1f). Together they comprise 2.1% of all proteins, 1.2% of all domains and 2.0% of all unique domain families in the human proteome. Two-hundred and seventy five of the domains are encoded by human disease genes with 108 domains containing annotated pathogenic variants. Across the dataset, mutations in the buried cores of the domains are more detrimental than mutations on their surfaces (Fig. 1d), with mutations to polar amino acids having stronger destabilizing effects in cores and mutations to hydrophobic residues having stronger effects on surfaces (Extended Data Fig. 1g). Mutations to proline are the most detrimental overall (Fig. 1d), in both core and surface residues, with highly destabilizing effects in beta strands and helices and milder effects in coils (Extended Data Fig. 1g,h). A subset of deep mutational scans of domains belonging to the most abundant families is shown in Fig. 2, and the full mutagenesis dataset is available in Supplementary Fig. 1.

7Our dataset represents a nearly fivefold increase in stability measurements for mutations in human proteins (Fig. 3a) and so provides an opportunity to evaluate how well computational VEPs predict changes in stability. Several general VEPs provide reasonable prediction of abundance changes (Fig. 3b,c), including the protein language model ESM1v (median ρ = 0.48) and the deep generative model EVE (median ρ = 0.48). Amongst tested dedicated stability predictors, the graph neural network ThermoMPNN performs best (median ρ = 0.50; Fig. 3b,c), and is the best predictor overall. Notably, all tested stability predictors perform poorly on small zinc-finger domains that require metal binding for stability (Extended Data Fig. 2a). After excluding zinc-fingers, ThermoMPNN is still the best performing method overall (median ρ = 0.57; Extended Data Fig. 2a), even when evaluating protein domains with no homology to the Megascale dataset used for ThermoMPNN training (median ρ = 0.57; Extended Data Fig. 2b).

8Fold stability is one of many biophysical properties that contribute to protein function and protein sequence conservation during evolution. In particular, many proteins bind other proteins, nucleic acids and small molecules or catalyse enzymatic reactions, and these functions often trade off with stability. The extent to which selection on protein sequences is driven by changes in stability rather than other biophysical properties is an important open question.

9To address this, we compared experimentally quantified stability to evolutionary fitness quantified by ESM1v across more than 500,000 variants in more than 500 domains (Methods). We find that protein stability accounts for a median of 30% of the variance in protein fitness when considering all domains (Extended Data Fig. 2e). However, the contribution of stability to fitness varies across domain families, with stability making a larger contribution to the fitness of all-beta domains (40% of the variance) than to that of all-alpha (25%) and mixed domains (25%) (Extended Data Fig. 2f). This is consistent with the lower structural tolerance to mutations of beta sheets compared with alpha helices, and suggests that stability is a more important determinant of the fitness landscapes of all-beta domains. The overall contribution of stability to protein fitness also varies in domains with different molecular functions. For example, stability has a larger contribution to fitness in SH3, WW, PHD finger and Pointed protein–protein interaction domains, and a lower contribution in DNA-binding HMG-box domains, nuclear hormone-type C4 zinc-fingers and homeodomains (Extended Data Fig. 2g).

10Mutations in binding interfaces, active sites and allosteric control sites typically have larger effects on function than can be accounted for by changes in protein stability,,. For example, quantifying the effects of mutations on protein binding and abundance enables binding interfaces and allosteric sites to be comprehensively mapped,. We reasoned that a similar approach could be used to identify functional sites in hundreds of domains by combining our abundance measurements with evolutionary fitness quantified by ESM1v.

11For individual domains, protein abundance is non-linearly related to evolutionary fitness predicted by ESM1v (Fig. 3d and Extended Data Fig. 3a). We used sigmoidal curves to model this relationship (n = 426 domains, ESM1v fitness range > 10, wild-type aPCA fitness percentile < 30; Fig. 3d and Extended Data Fig. 3a). The residuals to these fits identify mutations with larger or smaller effects on evolutionary fitness than can be accounted for by changes in stability (Fig. 3d–f).

12This analysis identifies a total of 102,231 mutations with larger effects on evolutionary fitness than can be accounted for by changes in stability (24% of the total; two-tailed z-test false discovery rate (FDR) < 0.1, normalized aPCA residual > 0.3). These mutations are enriched in known functional sites annotated in the Conserved Domains Database (CDD) (odds ratio (OR) = 2.72, Fisher’s exact test P < 2.2 × 10−16, n = 3,104 functional sites in 2,800 residues, Fig. 3g and Extended Data Fig. 3b). Defining evolutionary functional sites as residues with a weighted mean residual of less than 0.3 identifies a total of 5,231 sites in 426 domains. These sites are strongly enriched in CDD-annotated sites (OR = 4.50, Fisher’s exact test P < 2.2 × 10−16; Fig. 3h) and identify many known DNA-, RNA- and protein-binding interfaces (Fig. 3i,j). However, these evolutionary functional sites also include 1,942 sites in 180 domains without existing CDD annotations, and 1,873 additional sites in domains with other CDD annotations (Supplementary Table 3).

13Of note, evolutionary functional sites without known annotations are located in closer proximity to annotated functional sites (median side chain heavy atom distance (d) = 3.62 Å) than other residues (median d = 6.93 Å, P < 2.2 × 10−16). Many of these therefore act as ‘second-shell’ residues contacting DNA-, RNA- and protein-binding interfaces (Fig. 3f,j and Extended Data Fig. 3c), where sequence changes may indirectly impact binding via energetic interactions with interface residues,.

14Human Domainome 1 contains 3,652 variants with clinical annotations. Of these, 621 are classified as pathogenic or likely pathogenic (hereafter referred to as pathogenic), 322 are classified as benign or likely benign (hereafter benign) and 2,709 are classified as variants of uncertain significance (VUS), with 114 domains containing at least one pathogenic variant. Pathogenic variants are unevenly distributed across domains, with 75% of pathogenic variants contained in 25% of domains, and 41 domains containing only a single pathogenic variant (Extended Data Fig. 4a).

15In total, 380 out of 621 pathogenic variants (61%) cause a detectable domain destabilization (one-tailed z-test, normalized aPCA fitness < 0, FDR < 0.1) and 303 out of 621 (48%) are strongly destabilizing (FDR < 0.1, normalized aPCA fitness < −0.3; Fig. 4a,b). This contrasts with 129 out of 322 (40%) and 50 out of 322 (16%) of benign variants, respectively. However, the association between pathogenicity and destabilization varies across domain families (Fig. 4a,b and Extended Data Fig. 4b). For example, many pathogenic mutations in β/γ crystallins that cause cataract disease are strongly destabilizing (13 out of 18, OR = 11.98 compared with benign variants, Fisher’s exact test P = 1.99 × 10−6; Fig. 4b). By contrast, a smaller proportion of pathogenic variants are strongly destabilizing in homeodomains (OR = 3.94, P = 1.45 × 10−7), HMG-box domains (OR = 2.10, P = 0.045) and CUT domains (OR = 1.65, P = 0.55) (Fig. 4b), all of which bind DNA, suggesting that many pathogenic variants in these domain families affect biophysical properties beyond stability, such as DNA binding. Consistent with this variable association of destabilization with pathogenicity, aPCA has an overall lower performance than general VEPs in the classification of clinical variants across all domains, as the readout is specific for protein stability (Extended Data Fig. 4c,d).

16We next quantified the relationship between stability and pathogenicity in all individual domains with at least 20 annotated clinical variants (n = 17). Stability is the major contributor to pathogenicity in some domains. In the LIM domain 2 of FHL1, stability is an excellent classifier of pathogenic variants that cause reducing body myopathy (Matthew’s correlation coefficient (MCC) = 1; Fig. 4c), a disease caused by the accumulation of FHL1 aggregates. Similarly, the dominant ankyloblepharon–ectodermal defects–clefting syndrome is caused by mutations in the SAM domain of TP63 that lead to TP63 aggregation. Accordingly, most pathogenic variants in the SAM domain of TP63 are destabilizing (MCC = 0.83; Fig. 4c).

17By contrast, stability changes are a poorer predictor of pathogenic variants in other domains. A large proportion of mutations in the methyl-binding domain (MBD) of MECP2 that cause the dominant Rett syndrome are not destabilizing (MCC = 0.4; Fig. 4c). This suggests that many of these haploinsufficient variants interfere with the methylated-DNA-binding function of MECP2– without affecting the overall stability of the domain. Indeed, pathogenic mutations in MECP2 MBD that are not strongly destabilizing are concentrated in its DNA-binding interface (OR = 5.00, P = 0.09; Fig. 4d), in second-shell (OR = 4.09, P = 0.15) and in positively charged surface residues (OR = 4.35, P = 0.019), probably leading to a loss of binding affinity. Similarly, multiple mutations in the CRX homeodomain that cause inherited retinal dystrophies are not strongly destabilizing (MCC = 0.52; Fig. 4c). These variants are also enriched in its DNA-binding interface (OR = 6.33, P = 0.14; Fig. 4d) and in positively charged sites (OR = 3.58, P = 0.11). Notably, the mode of inheritance of mutations in CRX correlates with their stability effects: whereas the recessive R90W is strongly destabilizing (Fig. 4c), the dominant K88N and E80A are stable (Fig. 4c), consistent with their described gain-of-function mechanisms. This again suggests that destabilization is the major disease mechanism for some proteins and diseases but is much less important for others.

18Comparing across all domains with at least 20 clinical variants, there is a substantial correlation between how well stability explains pathogenicity and how well stability explains evolutionary fitness, as quantified by ESM1v (Pearson’s r = 0.78; Fig. 4f). We therefore used the correlation between stability and evolutionary fitness to rank all 114 domains with at least one known pathogenic variant (Extended Data Fig. 4f). Domains for which stability is highly predictive of fitness include many PHD finger and PDZ domains (Extended Data Fig. 4f). In these domains, we expect stability changes to be the major driver of pathogenicity. By contrast, stability only poorly predicts the fitness of homeodomains, HMG-box domains and nuclear hormone receptor-type zinc-finger domains (Extended Data Fig. 4f), suggesting that other molecular mechanisms will more frequently cause pathogenicity.

19The contribution of stability to protein fitness also varies among genes with different modes of inheritance and disease mechanisms (Fig. 4g). Recessive diseases are strongly associated with loss of function, whereas dominant diseases can also be caused by additional mechanisms such as gain-of-function and dominant-negative effects or toxic aggregation,. A median 44% of the variance in protein fitness is accounted for by stability changes in proteins mutated in recessive diseases, in contrast to only 26% in dominant disorders (P = 1.1 × 10−3, Wilcoxon rank sum test), suggesting that protein variants more frequently affect biophysical properties other than stability in dominant disorders. Indeed, loss-of-function and aggregation diseases are better explained by protein destabilization than ‘altered-function’ diseases associated with gain-of-function or dominant-negative mechanisms (Fig. 4g). Despite this overall association, however, we also find that within loss-of-function diseases, the extent to which destabilization explains protein fitness is variable (Fig. 4g). In summary, mutagenesis of more than 500 domains suggests that stability changes are an important cause of pathogenicity, but that this varies across proteins, with changes in stability being particularly important in recessive diseases.

20An important goal of establishing Human Domainome 1 is to quantify the extent to which mutational effects are conserved in structurally homologous proteins. Mutations in homologous sites might be expected to have very similar effects. However, mutations can also interact energetically, resulting in changes in mutational effects that depend on the sequence context, a phenomenon known as epistasis. If epistasis is prevalent, then mutational effects will be poorly conserved in divergent homologous proteins, limiting the extent to which experimental data from some proteins can be used to predict stability changes in homologues. Human Domainome 1 contains saturation mutagenesis data for 5 or more domains for 26 different domain families, enabling us to quantify the importance of epistasis for protein stability in structurally diverse protein folds.

21To quantify the conservation of mutational effects, we fitted a thermodynamic model based on the Boltzmann partition function to all of the data for a domain family (Fig. 5a). The model assumes that mutations cause the same change in folding energy (∆∆G) in all homologous domains, and that the energetic effects of mutations combine additively with no specific epistasis (Fig. 5a,b). We first fitted this model to homeodomains, the most abundant domain family in the dataset, with 36 human homologues. The model provides very good prediction of mutational effects across all 36 homeodomains (Pearson’s r = 0.78 by tenfold cross-validation; Fig. 5c,d). The inferred free-energy changes can be considered homologue-averaged mutational effects, providing an energy model describing the entire family (Fig. 5e,f). A linear model provides similarly good performance but with biased prediction residuals (r = 0.78; Extended Data Fig. 5a,b). We additionally evaluated the performance of the Boltzmann energy model by leaving out single homeodomains from the training dataset (n = 36 models) and found similarly good performance (median Pearson’s r = 0.74). Predictive performance was, as expected, better for domains with more similar sequence to the training dataset, but reasonable across a wide range of sequence divergence (Pearson’s r = − 0.75; Fig. 5g).

22Extending this analysis to all 26 families with at least 5 homologues in Human Domainome 1 (Fig. 5h) results in similarly accurate Boltzmann energy models for all families (median Pearson’s r = 0.80, median percentage of explainable variance = 80.7%; Fig. 5i and Extended Data Fig. 6a,b). The performance on individual left-out homologues is also very good (median Pearson’s r = 0.66, median percentage of explainable variance = 73.5%, Fig. 5j). For most domains, predictions were best for the domains most similar to the training dataset but were also good for domains with higher sequence divergence (Extended Data Fig. 7). Predictions using these energy models were better than those made with ThermoMPNN, the top-performing stability predictor on our dataset (Fig. 5j), and also had a good performance on stability deep mutagenesis scans generated using in vitro proteolysis selections (n = 38 homologues, median Pearson’s r = 0.65; Fig. 5k).

23The excellent performance of these additive energy models is both useful and important: it demonstrates that epistasis makes only a small contribution to protein stability across these levels of sequence divergence. Combinatorial mutagenesis of individual proteins suggests a similar conclusion. The decay of predictive performance with sequence divergence does, however, suggest a role for epistasis in the evolution of protein stability. Indeed, we identify 25,410 mutations with evidence of epistasis as variants with large residuals to the energy model fits (FDR < 0.1 two-tailed z-test, |residual| > 0.05 h−1; Extended Data Fig. 8a,b). These epistatic variants are enriched in the buried cores of protein domains (OR = 2.71, Fisher’s exact test P = 1.05 × 10−5) and depleted from protein surfaces (OR = 0.19, Fisher’s exact test P = 9.79 × 10−14; Extended Data Fig. 8c) across the full range of measured stabilities (Extended Data Fig. 8d) and protein families (Extended Data Fig. 8e). This data suggests that genetic interactions are more important for the evolution of protein cores, consistent with these residues having a larger number of structural contacts than solvent-exposed residues (Extended Data Fig. 8f).

24Finally, the good performance of the Boltzmann energy models across all 26 families suggests that we can use them to provide proteome-wide stability predictions for these domain families. Using the models we made predictions for an additional 4,107,436 variants in 7,271 domains, including an additional 13,878 clinical variants, of which 1,310 are pathogenic, 951 are benign and 11,617 are variants of uncertain significance (Extended Data Fig. 9a–c). Consistent with the results for clinical variants with experimentally measured stability changes, 686 (52%) of these pathogenic variants are predicted to reduce protein stability (one-tailed z-test, scaled ∆∆G > 0, FDR < 0.1, OR = 2.95 compared with benign), and 452 (34%) are predicted to have strongly destabilizing effects (scaled ∆∆G > 0.3, FDR < 0.1, OR = 5.95 compared with benign; Extended Data Fig. 9d,e). Similarly to the fitness-level data, our energy models outperform stability predictors in pathogenicity prediction, and have a lower performance than general VEPs (Extended Data Fig. 9f). These high-quality stability predictions provide a resource for the mechanistic interpretation of clinical variants for entire protein families and suggest a strategy for expanding Human Domainome 1 proteome-wide by experimentally mutagenizing representative examples for all families.

25We present here a large-scale experimental analysis of mutational effects in human proteins. Human Domainome 1 demonstrates the feasibility of quantifying human protein variant effects at scale and provides high-confidence measurements of changes in protein abundance for 563,534 variants in 522 structurally diverse protein domains. The dataset increases the number of stability measurements for human variants about fivefold and serves as a large and standardized reference dataset for the interpretation of clinical variants and for benchmarking and training computational methods.

26An important goal of this project was to evaluate the contribution of protein stability to genetic disease and sequence evolution. We found that approximately 60% of pathogenic missense variants reduce protein stability, but that this varies quite extensively across domains and diseases, with stability changes being particularly important in recessive disorders. Similarly, we found that the contribution of stability to evolutionary fitness varies across protein families, with a larger contribution in all-beta protein domains compared with other fold types. How this relates to protein evolvability and the importance of trade-offs between stability and other molecular functions will be interesting avenues for future research. A further goal of establishing Human Domainome 1 was to quantify the conservation of mutational effects in structurally homologous proteins. Fitting additive energy models to the data for domain families revealed that mutational effects on stability are largely conserved in homologous domains, with a small contribution from epistasis that increases with sequence divergence. This energetic additivity enables proteome-wide prediction of stability changes for entire protein families and suggests an efficient strategy to complete a first draft of the human domainome by mutagenizing representative examples for every protein family.

27To maximize diversity and experimental efficiency, Human Domainome 1 focused on structurally diverse small protein domains. An important caveat of this approach is the use of domains isolated from their native sequence context, and the extent to which mutational effects differ in full-length proteins is an important question for future work. Advances in DNA synthesis, assembly and mutagenesis, should facilitate this. Moreover, many domains did not yield sufficient signal in aPCA, indicating low stability or solubility, and other domains had mutational effects that were incompatible with two-state folding. In future work, it will be important to understand the behaviour of these domains in the assay and to test methods to increase the signal-to-noise such as altering domain boundaries, expression as full-length proteins, mutagenesis of exposed hydrophobic residues and the use of different reporters or solubilization tags. In this study, we have reported the effects of Human Domainome 1 variants on cellular protein abundance. However, the same libraries can be re-used in future work to quantify variant effects on additional molecular traits, for example abundance in different cellular contexts,,, in vitro fold stability, protein–protein, and protein–nucleic acid interactions, protein localization, aggregation and allostery,,. It will also be important to include more complex genetic variants such as insertions and deletions and to extend to proteome-wide coverage, including extracellular and transmembrane proteins.

28Human Domainome 1 is an important step in the comprehensive experimental analysis of human protein variants. It forms part of an ongoing global effort to determine the consequences of every mutation in every human protein and to produce reference atlases for the mechanistic interpretation of clinical variants. Beyond human genetics, it is also part of a broader effort to produce large, well-calibrated datasets that quantify how changes in sequence alter the biophysical properties of proteins. We believe that such multimodal biophysical measurements for millions of proteins and variants will enable machine learning approaches to be effectively brought to bear on the generative functions of molecular biology, allowing accurate prediction and engineering from sequence.

29LB: 10 g l−1 Bacto-tryptone, 5 g l−1 yeast extract, 10 g l−1 NaCl; autoclaved 20 min at 120 °C. YPDA: 20 g l−1 glucose, 20 g l−1 peptone, 10 g l−1 yeast extract, 40 mg l−1 adenine sulfate; autoclaved 20 min at 120 °C. SORB: 1 M sorbitol, 100 mM lithium acetate, 10 mM Tris pH 8.0, 1 mM EDTA; filter-sterilized (0.2 mm nylon membrane, ThermoScientific). Plate mixture: 40% PEG3350, 100 mM lithium acetate, 10 mM Tris-HCl pH 8.0, 1 mM EDTA pH 8.0; filter-sterilized. Recovery medium: YPD (20 g l−1 glucose, 20 g l−1 peptone, 10 g l−1 yeast extract) + 0.5 M sorbitol; filter-sterilized. SC-URA: 6.7 g l−1 yeast nitrogen base without amino acid, 20 g l−1 glucose, 0.77 g l−1 complete supplement mixture drop-out without uracil; filter-sterilized. SC-URA/ADE: 6.7 g l−1 yeast nitrogen base without amino acid, 20 g l−1 glucose, 0.76 g l−1 complete supplement mixture drop-out without uracil, adenine and methionine; filter-sterilized. MTX competition medium: SD–URA/ADE + 200 μg ml−1 methotrexate (BioShop Canada), 2% DMSO. DNA extraction buffer: 2% Triton X-100, 1% SDS, 100 mM NaCl, 10 mM Tris-HCl pH 8, 1 mM EDTA pH 8.

30We sampled PFAM-annotated human protein domains (version 34.0) of intracellular human proteins (not defined as ‘extracellular’, ‘transmembrane’ or ‘secreted’ in UniProt), prioritizing proteins that had at least one annotated pathogenic variant in ClinVar. We additionally included protein domains with in vitro stability measurements to benchmark the assay. We designed the library in two rounds: we first designed and selected two libraries (A1 and B3) containing single mutants of a total of 485 domains that had been previously tested to grow as wild types. To design a second set of libraries (C1 to C7), based on the results of A1 and B3, we excluded domains without a well-defined hydrophobic core (not having at least 10% of residues with rSASA <25%) and disordered domains defined as having an average AlphaFold2 pLDDT < 50, indicative of protein disorder. This second set contained 631 domains that had been previously tested to grow as wild type, and 132 additional domains not previously tested. The domain sequences were codon optimized with emboss backtranseq using the Saccharomyces cerevisiae codon usage table, and excluding the GCT alanine codon to prevent the appearance of HindIII restriction sites. The domain sequences were fused to 5′ gggctgctctagaatggctagc and 3′ aagcttggcggtggcgggtctg constant regions containing NheI and HindIII restriction sites for cloning.

31Libraries were synthesized by GCATbio using mMPS technology. The mMPS platform employs a streamlined four-module process—identification, sorting, parallel synthesis, and recycling—to enable an efficient and flexible workflow using traditional phosphoramidite chemistry for oligonucleotide synthesis. Silicon carbide wafers were used as the substrate, with each microchip uniquely engraved with a QR code, serving as both a chip identifier and a record of the predetermined sequence information. This allowed for precise sequence identification and assembly post-synthesis. Standard and degenerate bases (N/K, representing a mixture of A, T, C and G, or T and G, respectively) were employed for synthesis. The degenerate bases were pre-mixed in an optimized ratio to ensure equal representation of each nucleotide at the targeted mutation sites. Library A1, containing sequences shorter than 128 nt, was synthesized using single-stranded oligonucleotides directly cleaved from the microchip. By contrast, the remaining eight libraries (B3, C1 to C7), composed of longer sequences (sequence lengths ranging from 179 to 341 bases), were constructed via polymerase cycling assembly (PCA)-PCR. Initially, the sequences were divided into 378,896 sub-sequences, each ranging from 60 to 80 nt, for synthesis using the mMPS system. The full-length DNA fragments of each library were then assembled from these sub-fragments using PCA-PCR.

32Libraries were cloned at a variant coverage of ~100× or greater by restriction digestion and ligation into pGJJ162 (Supplementary Table 6), an aPCA assay plasmid where the DHFR3 fragment is fused at the C terminus of the target protein domains and the fusion is driven by the CYC promoter, and the DHFR1,2 fragment is expressed at high levels as is driven by the GPD promoter.

33Variant libraries were transformed in triplicate with a coverage of 20× or greater (with a mean coverage across libraries ~170×). For each transformation, we grew a 1 l YPDA culture of late log phase S. cerevisiae BY4741 cells (OD600 ~ 0.8–1), collected the cells by centrifugation (5 min, 4,000g), resuspended in 43 ml SORB medium and incubated for 30 min on a shaker at room temperature. Then, 875 μl of 10 mg ml−1 previously boiled (5 min, 100 °C) single-stranded DNA was added to the cells, followed by 17.5 μg library plasmid DNA and 175 ml plate mixture. The mix was incubated for 30 min on a shaker at room temperature. 17.5 ml DMSO were then added, and the cells were split in 50 ml Falcon tubes for 20 min heat shock at 42 °C. Following incubation, cells were pooled and collected by centrifugation, the supernatant was discarded with a pump, and cells were resuspended in 250 ml recovery medium and incubated at 30 °C for 1 h. Cells were then centrifuged for 3 min at 3,000g and transferred into 1 l SC-URA. Ten microlitres of this culture were immediately plated onto SC-URA selective plates to monitor transformation efficiency. The rest of the culture was incubated for one or two overnights at 30 °C.

34SC-URA cultures were used to inoculate a 1 l culture of SC-URA-ADE at an OD600 = 0.2–0.4, which was grown overnight (input culture). Cells from this culture were inoculated in 1 l SC-URA/ADE + 200 μg ml−1 MTX to select stably expressed protein domain variants. The remaining input cells grown SC -URA/ADE were collected and frozen for DNA extraction. MTX cultures were left to grow overnight to an OD600 = 1.6–2.5, corresponding approximately to 5 generations, collected and frozen for DNA extraction (outputs).

35Total DNA was extracted from yeast pellets equivalent to 50 ml of cells at OD600 = 1.6 as described in our previous work,. Plasmid concentrations in the resulting samples were quantified by against a standard curve of known concentrations by qPCR, using oGJJ152 and oGJJ153 as qPCR primers that amplify in the origin of replication of the aPCA assay plasmid.

36To generate the sequencing libraries, we performed two rounds of PCR amplification. In the first round, we used primer pools (oTB595+ and oTB748+; Supplementary Table 7) flanking the inserts that introduce frame-shifting nucleotides between the Illumina adapters and the sequencing region of interest. To maintain variant representation, we carried out eight 100 μl PCR1 reactions per sample, each of which starting with 125 million plasmid molecules that we amplified for 8 cycles. The reactions were column-purified (QIAquick PCR purification kit, QIAGEN), and the purified product was amplified further using the standard i5 and i7 primers to add the remainder of the Illumina adapter sequences and the demultiplexing indices (dual indexing) unique to each sample. We carried out a total of 8 100 μl PCR2 reactions per sample, each starting with 20–40 ng of purified product, that was amplified for 8 more cycles. The resulting amplicons were run on a 2% agarose gel to quantify the samples before pooling them for joint purification, and to ensure the specificity of the amplification and check for any potential excess amplification problems. The final libraries were size selected by electrophoresis on a 2% agarose gel, and gel-purified (QIAEX II Gel Extraction Kit, QIAGEN). The amplicons were subjected to Illumina paired-end 2 × 150 sequencing on a NextSeq2000 instrument at the CRG Genomics facility.

37FASTQ files from paired-end sequencing of all aPCA selections were processed with DiMSum v1.2.11 (https://github.com/lehner-lab/DiMSum) using default settings with minor adjustments. The option “barcodeIdentityPath” was used to specify a variants file in order to recover designed variants only (NNK mutations in the protein domains present in each sublibrary), and starting and final culture optical densities and selection times were specified to infer absolute growth rates across libraries.

38We computed several quality control metrics on a per-domain basis. First, the wild-type position in the fitness distribution, as the difference between the wild-type and the 95th percentile of the distribution, divided by the difference between the 95th and 5th percentile in the distribution (fitness 90% range). Second, the Pearson’s correlation coefficient between biological replicates for missense variants. Third, the correlation between fitness and solvent accessibility for all variants in each domain. And fourth, the correlation between fitness and hydrophobicity of all variants in each domain (measured as the first principal component of a comprehensive table of amino acid properties). All structural calculations of the domains analysed are based on AlphaFold2 predictions obtained from https://alphafold.ebi.ac.uk/.

39Compared to folded domains (more than 10% core residues and pLDDT > 50), disordered and domains without a well-defined hydrophobic core (see ‘Library design’) behave in a very distinct fashion in aPCA, with the wild type located in the middle of the fitness distribution, hydrophobicity negatively correlated with protein abundance, and narrow dynamic ranges resulting from typically small effects of mutations (Extended Data Fig. 1c). We thus excluded these domains from further analysis in the main text. Folded domains showed larger dynamic ranges, with the wild-type among the fittest variants, a negative correlation between mutation sensitivity and solvent accessibility, and a lack of correlation to hydrophobicity (Extended Data Fig. 1c). For folded domains, the four quality metrics described above were highly correlated with each other, and we combined them using principal component analysis. The first principal component explained 64% of the variance, and was used as a single quality metric to rank all domains.

40We retained all domains ranked 600 or higher, with a missense Pearson’s r > 0.485, and with more than 50% of variants measured (with at least 10 counts in at least 1 replicate), resulting in 538 domains. We removed from this set 16 additional domains that were not compatible with 2-state folding, either showing many large effect abundance-increasing mutations (O75364_PF00046_64, O75956_PF09806_73, P10242_PF00249_89, P52952_PF00046_140, Q13263_PF00643_205, Q86TZ1_PF13181_60, Q8IX03_PF00397_1, Q8NDW8_PF13181_799, Q9Y2H9_PF17820_968, Q9Y6M9_PF05347_15), narrow fitness ranges (P35637_PF00641_421, EHEE-rd2-0005_1, HHH-rd2-0133_1), or clear correlations to hydrophobicity (E9PAV3_PF19026_2040, HEEH-rd3-0223_1, Q5VTD9_PF00096_193), resulting in a final set of 522 domains.

41The distributions of quality metrics for retained and discarded folded domains, and for disordered domains and folded domains without a well-defined core are in Extended Data Fig. 1c. Of the 478 discarded folded domains, 224 had slow growth rates as wild type (<0.075 h−1), 176 had an unfit wild type (wild-type position <0.4), 107 were negatively correlated with hydrophobicity (Pearson’s r < −0.25), and 158 were not correlated with solvent accessibility (Pearson’s r < 0) (we note that many of the discarded domains meet more than one of these criteria). Discarded domains are shorter (median length = 49 amino acids (aa)) than retained domains (length = 58 aa, p = 1.87 × 10−14, Wilcoxon rank sum test), and are enriched in zinc-finger domains (57% of zinc-finger domains discarded in contrast to 43%, 41%, and 39% for all-alpha, all-beta, and alpha+beta, respectively, P = 2.33 × 10−6, Fisher’s exact test). Multimerization status only affected success rate marginally, as 49% of domains from proteins that form complexes were discarded, compared to 53% of those that do not engage in protein–protein interactions (P = 0.2, Fisher’s exact test). Finally, domains with yeast orthologues were less likely to be retained (45%) than domains without yeast orthologues (53%, P = 0.08, Fisher’s exact test).

42We normalized the growth rates and growth rate errors within each protein domain by linearly scaling the data such that the wild type-normalized fitness equals zero, and the 2.5th percentile of the distribution of growth rates of all missense variants plus the wild type is equivalent to a normalized fitness of −1. In all box plots shown in the manuscript, the central line represents the median, the upper and lower hinges correspond to the first and third quartiles, with the upper whisker extending from the hinge to the largest value no further than 1.5 × IQR, and the lower whisker extending from the hinge to the smallest value no further than 1.5 × IQR.

43We generated a structural distance matrix based on Foldseek hit probabilities. We computed all pairwise alignments between the domains in the final retained set by first creating a Foldseek database ‘foldseek_domainome_db’ containing all domains, and then searching all domains against the database using foldseek easy-search *.pdb foldseek_domainome_db foldseek_easy_allvsall tmp --format-output “query,target,alntmscore,qtmscore,ttmscore,prob” --exhaustive-search TRUE -e inf.

44To visualize the similarity network, we imported the domains as nodes and the Foldseek probabilities as edges into Gephi, and applied two layout algorithms: first the Fruchterman Reingold algorithm to equilibrium, which resulted in a clear separation of the different SCOP classes, followed by a Force Atlas 2 layout algorithm (preventing node overlap), which accentuated the separation between the different domain family clusters within each SCOP class.

45To compare aPCA measurements with in vitro ∆∆G values, we used domains containing at least 10 variants in more than a single residue measured in vitro and in aPCA, and with a range of at least 2 kcal mol−1 for the in vitro measurements, and of 0.075 h−1 log growth rate units for the aPCA measurements (n = 10 domains). To compare aPCA measurements with ∆∆G values derived from high-throughput proteolysis deep mutational scans, we used protein domains with an overlap (overlapping length/total alignment length) of at least 80%, and an aPCA measurement range of 0.075 h−1 log growth rate units (n = 12).

46The total number of human missense variants measured in the Megascale dataset was retrieved from the high-confidence dataset excluding double mutants and domain duplicates (from different Protein Data Bank (PDB) entries and genetic backgrounds of the same PDB entry). We also quantified the total number of human missense variants measured in all abundance (VAMP-seq, aPCA) datasets available (as of 7 July 2024) in MaveDB total number of human missense variants in ProthermDB.

47We generated ESM1v predictions (https://github.com/facebookresearch/esm) using the domain sequences alone as input, and using sliding windows across the human proteome (size = 1,000 aa, step = 250 aa). For each domain, we used the predictions corresponding to the window in which the domain is most centred. We obtained precomputed AlphaMissense– predictions from https://console.cloud.google.com/storage/browser/dm_alphamissense. We obtained EVE, popEVE and Tranception precomputed scores from https://pop.evemodel.org/. ‘EVE domain’ scores computed on high-coverage alignments generated specifically for domainome domains were kindly shared by A. Kollasch and D. Marks. We used precomputed RaSP scores, and generated DDMut and ThermoMPNN predictions with the aid of AlphaFold2 structure predictions. FoldX predictions were obtained from https://ftp.ebi.ac.uk/pub/databases/ProtVar/predictions/stability/. We used Spearman’s ρ to quantify the relationship between the predictions and aPCA fitness. Domains with homology to protein domains in the Megascale dataset were defined using hmmer (http://hmmer.org/) hmmscan against PFAM, and predictors were evaluated on a homologue-free set to prevent leakage from ThermoMPNN training.

48To estimate the fraction of the variance in mutational effects on evolutionary fitness that are attributable to protein stability, we calculated the correlation between ESM1v fitness predictions on full-length protein sequences and aPCA scores. The correlation coefficient was adjusted by the measurement error of the aPCA scores according to the Spearman disattenuation formula:where Rxy is the disattenuated correlation coefficient, rxy is the observed correlation coefficient, rxx is the mean correlation coefficient between aPCA biological replicates. This procedure was applied to both linear (Pearson’s) and rank (Spearman’s) correlation coefficients. The two were highly correlated (r = 0.94) and we report the Pearson’s r-based version for ease of interpretation.

49We carried out this analysis for domains where (1) the wild type is above the 30th percentile in the fitness distribution; and (2) the range between the 5th and 95th percentile of the distribution of ESM1v-predicted fitness is greater than 10 (n = 426 domains). We used sigmoid curves to model the relationship between normalized aPCA fitness (stability) and predicted function (ESM1v) of all variants in each individual domain, with an upper bound of 0 (wild-type aPCA-normalized fitness) and a lower bound of −1:Where faPCA is the aPCA-normalized fitness, fESM1v is the ESM1v-predicted fitness, xmid is the midpoint of the sigmoid and scal is the steepness parameter. To prioritize fitting the low stability variants, we weighted the fit by the aPCA-normalized fitness:

50We used a two-tailed z-test to identify mutations whose effects on fitness cannot be accounted for by stability effects. We calculated z-scores as the aPCA residuals to the fit divided by the aPCA error, derived P values based on the normal distribution, and performed multiple testing correction using Benjamini–Hochberg’s FDR. We additionally calculated per-residue mean residuals weighted by the aPCA-normalized fitness error.

51We obtained residue-level functional site annotations corresponding to the CDD using the InterPro API. Second-shell residues were defined as those with a minimum heavy atom distance of 5 Å to functional site residues.

52We identified destabilizing variants using a one-tailed z-test. z-Scores were calculated as the normalized fitness divided by the normalized fitness error, a P value was derived on the basis of a normal distribution, and FDR multiple test correction was applied. Destabilizing variants were defined as variants with FDR < 0.1, and strongly destabilizing variants as FDR < 0.1 and a normalized aPCA fitness < − 0.3. We obtained clinical variant annotations from ClinVar (January 2024 version) and from the UniProt API. Variant annotations were highly consistent between the two sources, and were merged. We tested for enrichments of clinical variant classes in stability classes using a two-tailed Fisher’s exact test.

53To estimate to what extent destabilization explains pathogenic mutations in individual domains, we used MCC. To increase statistical power, we included gnomAD variants with an allele frequency >10−5 as benign. We estimated the errors in MCCs by resampling the dataset based on the mean and errors of aPCA scores, 10 times. To analyse the distribution of pathogenic mutations in the MECP2 methyl-binding domain and the CRX homeodomain, we generated AlphaFold3 predictions. The sequences of the domains and DNA strands used for the predictions are available as Supplementary Table 8. DNA-binding interface residues were defined using getcontacts, and second-shell residues were defined as those with a minimum heavy atom distance of 5 Å to binding interface residues.

54To extend the analysis to a larger number of domains with small numbers of pathogenic variants, we used the fraction of variance in ESM1v-predicted fitness explained by aPCA scores as described above. MCCs based on clinical and gnomAD variants correlated well with the fraction of evolutionary fitness explained by stability effects in domains with at least 20 clinical and gnomAD variants (r = 0.78), validating the approach. We estimated errors in MCCs by resampling the aPCA fitness dataset using the mean and error of each variant, and recalculating the MCCs on clinical variants (n = 10 samples). To further validate the estimates of the contribution of stability to fitness, we estimated the fraction of variance in ESM1v-predicted fitness explained by thermoMPNN predicted stabilities for non-zinc-finger domains with at least 1 pathogenic variant (n = 86). These correlated well with the original estimates derived using our abundance data (Pearson’s r = 0.81; Extended Data Fig. 2d).

55Modes of inheritance and mechanisms of disease information were obtained from OMIM (https://www.omim.org/). To analyse the contribution of stability changes to disease according to mode of inheritance and mechanism of action controlling for protein composition (Extended Data Fig. 4e), we modelled the contribution of secondary structure and the percentage of core residues to the scores using linear models. We then extracted the model residuals as the composition-corrected scores.

56We used the pROC R package to generate receiver operating characteristic (ROC) curves and calculate ROC area under the curve. We incorporated gnomAD variants with an allele frequency >10−5 as benign. To test the performance when combining structural and sequence features (secondary structure, rSASA, wild-type residue, mutated residue), ESM1v and aPCA scores, we used generalized linear models in R. In addition to training and evaluating on the full dataset, we trained the logistic regression models with 90% of the data and evaluated on the remaining 10% unseen data.

57We used MoCHI to fit 2-state thermodynamic models on a per-family basis (for families with at least 5 human homologues and variants with a mean count >29). We specified a neural network architecture consisting of a single additive trait layer for shared folding energies across the family and a shared linear transformation layer. We used a ‘TwoStateFractionFolded’ transformation derived from the Boltzmann distribution function that relates energies to proportions of folded protein. To map homologous positions across families, we used PFAM alignments of human domains. To input into MoCHI, we recoded each domain wild-type sequence as an indel sequence as long as the PFAM alignment, plus additional positions in the alignment with variants that encode the wild-type identity of each homologue. Mutations in each domain were encoded in the corresponding alignment position. This design allows MoCHI’s one-hot encoding of both the mutations and the wild-type identities simultaneously, and the joint fitting of mutation ∆∆G and starting ∆G of each homologue. We trained the model using a tenfold cross-validation approach and evaluated the performance on held-out data.

58We additionally fitted a linear model, and a ‘TwoStateFractionFolded’ model with shared energies across all homologues of a family but with domain-specific linear transformations that account for potential differences in solubility of the folded and unfolded states between domains of the same family. These models resulted in highly correlated inferred energies to the original Boltzmann model (Pearson’s r = 0.972 for the linear model, and r = 0.938 for the Boltzmann model with domain-specific linear transformations) and similar performances on held-out data across families. We discarded the linear model due to highly biased residuals, and chose the ‘TwoStateFractionFolded’ model without homologue-specific linear scaling as the simplest of the remaining models. To compare across protein families, energies were rescaled such that the 2.5th percentile of the distribution of energies is equivalent to a scaled ∆∆G = −1 and the wild type to a ∆∆G = 0.

59We further evaluated the performance of the TwoStateFractionFolded model in families with at least ten homologues by training the model leaving out a single domain at a time. The genetic distance (hamming distance and BLOSUM62 distance) of each left-out homologue to each of the homologues that went into training was calculated and averaged, to compare to model performance. The fraction of explainable variance in aPCA scores accounted for by the models was estimated as the R2 between observed and predicted fitness divided by the R2 between replicates.

60We identified epistatic variants with large residuals to the MoCHI fits using a two-tailed z-test. We calculated z-scores as the residuals to the fit divided by the aPCA fitness errors, derived P values on the basis of a normal distribution, and applied FDR multiple testing correction. Epistatic mutations were defined as those with |residual| > 0.05 and FDR < 0.1. Enrichments of alignment sites in epistatic variants were calculated, and significantly enriched sites were identified using a two-tailed Fisher’s exact test. Epistatic sites were defined as those with a log2(OR) > 1.5 and a Fisher’s Exact test FDR < 0.05. We classified alignment sites as core residues (rSASA <25% in at least 75% of homologues), surface residues (rSASA > 25% in at least 75% of homologues), or changing residues (the rest). We defined contacts using getcontacts.

61Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

62Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41586-024-08370-4.