Assessing Variants of Uncertain Significance Implicated in Hearing Loss Using a Comprehensive Deafness Proteome

Hearing loss is the leading sensory deficit, affecting ~ 5% of the population. It exhibits remarkable heterogeneity across 223 genes with 6,328 pathogenic missense variants, making deafness-specific expertise a prerequisite for ascribing phenotypic consequences to genetic variants. Deafness-implicated variants are curated in the Deafness Variation Database (DVD) after classification by a genetic hearing loss expert panel and thorough informatics pipeline. However, seventy percent of the 128,167 missense variants in the DVD are “variants of uncertain significance” (VUS) due to insufficient evidence for classification. Here, we use the deep learning protein prediction algorithm, AlphaFold2, to curate structures for all DVD genes. We refine these structures with global optimization and the AMOEBA force field and use DDGun3D to predict folding free energy differences (ΔΔGFold) for all DVD missense variants. We find that 5,772 VUSs have a large, destabilizing ΔΔGFold that is consistent with pathogenic variants. When also filtered for CADD scores (> 25.7), we determine 3,456 VUSs are likely pathogenic at a probability of 99.0%. These VUSs affect 119 patients (~ 3% of cases) sequenced by the OtoSCOPE targeted panel. Approximately half of these patients previously received an inconclusive report, and reclassification of these VUSs as pathogenic provides a new genetic diagnosis for six patients.


Introduction
Hearing loss is the most prevalent sensory de cit, affecting approximately 5% of the world's population. In its evaluation, following an audiogram, genetic sequencing with a multi-gene panel is recommended as the most informative diagnostic test for infants and children with hearing loss ( In each patient screened, an average of 545 genetic variants is identi ed (Shearer et al. 2013). Ascribing a pathogenic consequence to these variants is challenging and requires deafness-speci c expertise. To help meet this challenge, we developed the Deafness Variation Database (Azaiez et al. 2018) (DVD). This resource includes 128,167 missense variants, which are classi ed by a genetic hearing loss expert panel and thorough informatics pipeline into one of ve categories: benign (B, n = 1,725), likely benign (LB, n = 27,907), likely pathogenic (LP, n = 2,441), pathogenic (P, n = 6,328), and variant of uncertain signi cance (VUS, n = 89,766). If a variant is classi ed as a VUS, a de nitive diagnosis cannot be made for patients affected by that variant.
For variant reclassi cation, additional studies are required and can include family segregation analysis, identi cation of the variant in a family member with hearing loss or an unrelated proband, or speci c wet lab based functional evidence (Richards et al. 2015). Given the disproportionate number of VUSs, making genotype-phenotype correlations from such evidence is infeasible. Therefore, we sought to apply deep learning-based protein structure prediction , atomic resolution simulation (Tollefson et al. variants classi ed as VUSs to determine whether it would be possible to reclassify some VUSs as P.
In 2019, protein structures of deafness-associated genes were known for fewer than 40% of all proteins and missense variants implicated in hearing loss (Tollefson et  accompanying prediction of protein misfolding, abrogated function and possible degradation can be done on a deafness proteome wide basis. However, computing ∆∆G Fold using protein structures from AlphaFold2 as input to rigorous molecular dynamics-based simulation for all 128,167 missense variants listed in the DVD is currently intractable due to computational expense. As an alternative, we use a high-throughput in silico tool to predict ∆∆G Fold (Guerois et Zhou and Zhou 2002) and identify VUSs most likely to induce signi cant protein misfolding (often ∆∆G Fold >2-3 kcal/mol), potentially allowing these variants to be classi ed as P. First, we use AlphaFold2 to curate full-length, isoform-speci c protein structures for all genes in the DVD (OtoProtein2). We then reduce biophysical inaccuracies (i.e., steric clashes and side-chain errors) in the OtoProtein2 structures by re ning them with an amino acid side-chain optimization algorithm (Tollefson et al. 2019)  Trained on experimentally known protein structures from the Protein Data Bank (PDB) (Berman et al. 2000), the AlphaFold2 neural network predicts protein structures from amino acid sequences to an accuracy comparable to experimental results using two modules ). The rst module develops a general hypothesis for the protein's structure in part from relationships between co-evolving amino acids associated with a multiple sequence alignment. The second module predicts the spatial relationships between subsequent amino acids to produce an explicit three-dimensional protein structure. By default, the two modules are generally applied in three iterative cycles to re ne the structure prediction; however based on prior work(Mirdita et al. 2022) we applied the modules in 15 cycles to achieve higher quality predictions.
Biophysical Re nement of the AlphaFold2 Deafness Proteome Using thermodynamic (see supplementary information) observations, we identi ed a ∆∆G Fold threshold to predict genetic variants that induce signi cant misfolding, loss of function and possibly protein degradation. We used classi ed DVD variants to determine the positive predictive value (PPV) of this ∆∆G Fold threshold. We applied this threshold to all P variants to determine which P variants are deleterious due to protein misfolding. We further applied this threshold to all VUSs in the DVD to determine which VUSs most likely impact protein misfolding and are therefore most likely to be P.

Integrating CADD Scores with ∆∆G Fold to Prioritize Variants
We combined the ∆∆G Fold predictions and threshold with CADD (Rentzsch et al. 2018) scores to prioritize VUSs most likely to be deleterious. Because variants with higher CADD scores are predicted to be more damaging(Rentzsch et al. 2018), we anticipated variants with both a large ∆∆G Fold and a high CADD score are more likely to be P. We set the CADD score threshold (25.7) to re ect a 99% PPV for classi ed DVD variants to be P when both ∆∆G Fold and CADD scores are combined. We then applied both the CADD threshold and the ∆∆G Fold threshold to identify VUSs that are deleterious with 99% certainty.

Curating Variant Features for Further Analysis
In addition to annotating ∆∆G Fold and CADD scores for each DVD variant, we aggregated features from the optimized structures to be used for variant analysis, prioritization, and deep learning. For each variant, we collected AlphaFold2's con dence in the protein structure at that variant's position, which can be used to prioritize analysis of variants in regions where protein structure is predicted with a high degree of con dence. Similarly, because amino acids buried within a protein domain are often intolerant of variation as compared to amino acids on the surface of a protein domain, we computed the percent of solvent accessible surface area (SASA) for each DVD variant. Finally, previous work has shown that minor allele frequency (MAF) can be used to classify common variants as LB in deafness-associated genes (Shearer et al. 2014); therefore, we included the MAF for each variant in the dataset of variant features.

Quality and Characteristics of Deafness Protein Structure Predictions
Using AlphaFold2, we developed complete protein structures for all genes and relevant isoforms in the Deafness Variation Database(Azaiez et al. 2018) (DVD, Fig. 1a, b). Called OtoProtein2, this dataset increases structural coverage of the deafness proteome from approximately 30% by experimental and homology protein structures curated during prior work(Tollefson et al. 2019) (i.e., called OtoProtein) to 100% (Fig. 1c). For each amino acid in a prediction, AlphaFold2 provides a unitless con dence score ranging from 1 to 100, with higher scores corresponding to higher con dence in the prediction. Model con dence is > 70 for 64% of wild-type amino acids and 60% of missense variant locations in the deafness proteome. The remaining amino acids and missense variants fall in regions that are predicted only with low con dence (i.e., con dence < 70).
[ Figure 1 Here] Approximately 41% of missense variants in the deafness proteome belong to a functional protein domain as characterized by InterPro (Apweiler et al. 2001;Blum et al. 2021), while 59% belong to exible termini, natively disordered regions, or uncharacterized domains (Table 1). InterPro characterized domains are enriched in high con dence protein structures, while natively disordered regions exist in lower con dence regions. Of the 128,167 missense variants in the deafness proteome, 34% belong to both a characterized domain and a high con dence structural region. Although missense variants are evenly distributed across InterPro characterizations (e.g., 41.3% and 41.4% of wild-type amino acids and missense variants are in a characterized domain, respectively), benign and likely benign variants favor lower con dence, uncharacterized regions while pathogenic and likely pathogenic variants favor higher con dence regions with functional protein domains (Tables 2, S1 and S2).  (Table 3).
The re nement procedure improved the dataset's mean MolProbity score from 2.86 to 0.97 (Fig. 3), making the OtoProtein2 structural quality equivalent to experimental structures at atomic resolution. Table 3 Average MolProbity re nement statistics for all deafness associated protein models in OtoProtein2 before and after optimization with Force Field X. A lower clash score, a lower percentage of poor rotamers, a higher percentage of favored backbone phi/psi angles, fewer backbone outliers and lower MolProbity score are each better. regions of a protein structure (i.e., often functional regions) have a higher mean ∆∆G Fold and a wider distribution of ∆∆G Fold than variants that fall within low con dence regions (i.e., often natively disordered protein regions).
[ Figure 3 Here] Using thermodynamics principles (see derivation in supplementary information), a ∆∆G Fold of > to be deleterious. Higher CADD scores are associated with P and LP variants (Fig. 4a). These variants also favor protein regions with high con dence (Fig. 4b) and consist primarily of domains and motifs that are intolerant to variation. Establishing a CADD threshold independently has a reasonable PPV (e.g., a CADD cutoff of 20 results in a PPV of 88.3%). We applied a CADD cutoff of 25.7 and combined this threshold with the ∆∆G Fold threshold, which resulted in a PPV of 99% and a speci city of 99.5%. While these stringent CADD and ∆∆G Fold thresholds limit prioritization to 3,456 destabilizing VUSs (Tables 4 and S4), these VUSs can be classi ed as LP due to protein misfolding (Fig. 4C).
[ Figure 4 Here] We found that P and LP variants are often in buried residues (i.e., solvent accessible surface area near zero percent) with con dent structure regions (Fig. 5ab). The prioritized dataset of 3,456 VUSs are consistently present in buried, con dent regions of the OtoProtein2 structures (Fig. 5c). Additionally, ∆∆G Fold , CADD scores, solvent accessible surface area, and structure con dence from the OtoProtein2 models for all variants in the DVD can be utilized for deep learning applications or for variant analysis.

Discussion
The classi cation of genetic variation in relationship to a disease phenotype is challenging. For hearing loss, the DVD uses an expert panel and rigorous informatics pipeline to classify changes in deafness-associated genes based on evidence of pathogenicity. This database includes over 128,167 missense variants, the majority of which (> 70%) are classi ed as VUSs due to insu cient evidence to classify as P or B. A VUS classi cation is problematic for both the healthcare provider and the patient as a de nitive diagnosis cannot be made. Using these ∆∆G Fold and CADD thresholds, we identi ed 3,456 VUSs that are LP due to protein misfolding.
Of these 3,456 prioritized VUSs, we have observed 79 across 119 patients who underwent comprehensive genetic testing using OtoSCOPE. Over half of these patients (60 patients) previously received an inconclusive genetic diagnosis. In ve patients with variants affecting autosomal recessive genes, the proband carried a second LP/P variant in the gene. Segregation analysis (SA) con rmed that the second LP/P variant occurs on the opposite allele in three of ve patients; in the remaining two patients, SA was not available. One patient carried a variant affecting an autosomal dominant gene. The work here delivers a de nitive genetic diagnosis for these six patients and directly impacts their subsequent healthcare (Table 5). For example, patient six carried a known P variant in TMPRSS3 in trans with a novel missense variant predicted to cause protein destabilization by this work (Fig. 6). The phenotype of the patient's hearing loss is highly speci c for TMPRSS3-related hearing loss (DFNB8/10). Reclassi cation of patient six's novel missense variant from VUS to LP results in a de nitive genetic diagnosis, ultimately directing subsequent medical care and recurrence risk calculations for offspring. Current guidelines established by the American College of Medical Genetics and Genomics (ACMG) for hearing loss do not incorporate in silico ∆∆G Fold calculations, however, our work demonstrates the utility of protein modeling for hearing loss diagnostics. Further work is indicated to guide incorporation of protein modeling into ACMG guidelines for hearing loss and deafness.
[ Figure 6 Here] Table 5 Patients with de nitive diagnoses from upgraded classi cation of priority VUSs. Segregation analysis con rms that the second variant occurs on the opposite allele in three probands.

Yes
The number of prioritized VUSs and impacted patients is greatly affected by adjustments to the ∆∆G Fold and CADD thresholds. We used a ∆∆G Fold threshold of 1.8 kcal/mol and a CADD threshold of 25.7 to reach a PPV of 99.0% (false positive rate < 0.5%), but by increasing the CADD threshold to 30.0, the PPV approaches 100%. These stringent thresholds leave negligible room for a false positive diagnosis but provide a prioritized dataset of only 419 VUSs that are LP. Seven of these 419 VUSs impact 18 OtoSCOPE patients. Alternatively, a more lenient PPV of 95% is reached by disregarding CADD scores and dropping the ∆∆G Fold cutoff to 1.0 kcal/mol. These parameters provide a substantially larger dataset of 12,585 VUSs that are LP, albeit with a 5.6% false positive rate and impact 775 OtoSCOPE patients.
Though we applied the ∆∆G Fold and CADD thresholds on a deafness-proteome-wide scale, these cutoffs can be tuned to better t a protein, domain, or amino acid speci c level. Biochemical, environmental, and structural differences contribute to a protein's ability to tolerate changes to its structure. For example, ACTG1 encodes gamma actin, a highly conserved cytoskeletal protein; even small ∆∆G Fold signi cantly disrupt gamma actin's structure and function. While no P ACTG1 variants from the DVD surpass both the ∆∆G Fold and CADD cutoffs for the proteome-wide scale, a smaller ∆∆G Fold threshold may detect subtle structural changes that will affect gamma actin's highly conserved structure. Similarly, different domains within an individual protein can bene t from domain-speci c ∆∆G Fold analysis. Cochlin, the protein product of the COCH gene, has one Limulus factor C (LCCL) domain and two Von Willebrand factor A (VWFA) domains. P variants in COCH are known to localize in the LCCL and second VWFA domains (Gallant et al. 2013). Known P variants aggregating in just one of cochlin's two VWFA domains demonstrate the need for domain-speci c analysis to identify which domains are more sensitive to amino acid variation and are intolerant of misfolding. Even individual amino acid characteristics such as the structural con dence of the wild-type amino acid, SASA, or number of hydrogen bonds can affect an amino acid's ability to tolerate a missense variant that disrupts the protein's structure. As approaches for ∆∆G Fold predictions are improved, context-dependent thresholds will be signi cant for variant interpretation.
The ∆∆G Fold and CADD thresholds used to identify VUSs that induce substantial protein destabilization can also provide an estimate of the number of deafness-causing genetic variants yet to be classi ed as P. Because ∆∆G Fold quanti es the disruption to protein folding induced by variants, ∆∆G Fold resolves only those VUSs that are P due to protein misfolding. Applying these thresholds to listed P and LP variants in the DVD allows us to identify that subset of missense variants that destabilize protein structure. Of the 6,328 known P predictions often require that monomeric proteins be predicted in segments. This memory limitation is only exacerbated by the prediction of protein complexes where memory limits are more easily reached. Nevertheless, attaining a comprehensive model of the deafness interactome and subsequent analysis of ∆∆G Bind will be the subject of future studies. The analysis of indels, non-coding variants, and other variants, are beyond the scope of our current work, however, prioritization and characterization of these variants should be considered in context with the VUSs prioritized herein. Regardless of the work remaining, the deafness proteome and ∆∆G Fold analysis we present has revealed trends for P variants and provides insight on VUSs that are LP due to protein misfolding.
In summary, by using ab initio protein structure prediction, optimization, and thermodynamic analysis, with 99% con dence, we have identi ed 3,456 VUSs that are LP in patients with hearing loss due to protein misfolding. The deafness protein structures developed here have been incorporated with the DVD to inform deafness-associated variant analysis. As atomic resolution protein structures and in silico variant analysis techniques progress, continued and re ned analysis of free energy differences for deafness-associated variants will inform pathogenicity classi cations and lead to enhanced patient diagnoses. All data accumulated during this project are available on Github (https://github.com/SchniedersLab/OtoProtein). Figure 1 Structures and quality of proteins implicated in deafness. AlphaFold2's novel predicted protein regions are color coded by con dence in the prediction. Gray domains represent homology or experimental structures curated in prior work for a) cochlin and b) stereocilin. a) The root-mean-square deviation (RMSD) of the LCCL and vWFA domains of cochlin (COCH) from AlphaFold2's domain predictions to the previous models are shown in parentheses. b) This work increased protein structural coverage of stereocilin (STRC) from 12% to 100%. c) Structural model coverage of wild-type amino acids and missense variants for the entire deafness proteome shows that this work increased coverage from <30% (gray, prior work) to 100% coverage. The stacked bars are color coded based on con dence in the protein structure. The wild-type amino acids and missense variants in the deafness proteome are present in similar proportions across all structural con dence ranges, indicating that speci c con dence regions are not enriched for the presence of missense variants MolProbity score histogram for the OtoProtein2 database. Before optimization (red), the mean MolProbity score of the models is 2.86 and after optimization (blue) the structures are consistent with atomic resolution at a mean MolProbity score of 0.97. MolProbity scores are calibrated to re ect the expected crystallographic resolution of the diffraction dataset employed to create a protein structural model (i.e., a MolProbity score of 1.0 indicates that the structure is consistent with 1.0 Å resolution X-ray diffraction data) Figure 3 The range of ∆∆G Fold predictions for missense variants in the Deafness Variation Database (DVD). a) Box plots are grouped based on DVD pathogenicity classi cation and bars are colored based on the structure con dence at the variant's amino acid position. Pathogenic variants and variants in con dent portions of protein models have a larger distribution of ∆∆G Fold than the benign and low con dence (e.g., usually solvent exposed) counterparts. The number of observations belonging to each box is printed below the box. b) A box plot for all VUSs in the DVD. Each outlier in the boxplot can represent multiple VUSs due to overlap in ∆∆G Fold .

Figures
Outliers colored in red are prioritized VUSs that have a large ∆∆G Fold (≥1.8 kcal/mol) and a high CADD score (>25.7). Unprioritized VUSs do not have a high CADD score. The number of prioritized VUSs belonging to each box is printed in red below the total number of observations belonging to the box Figure 4 Prioritizing variants of uncertain signi cance (VUSs) from thermodynamic data and CADD scores. Folding free energy differences (∆∆G Fold ) versus CADD score for all missense variants observed in the Deafness Variation Database (DVD) with points colored according to DVD classi cation (panels a and c) or model con dence at the variant's amino acid position (panel b). CADD score and ∆∆G Fold show a positive correlation. A high ∆∆G Fold and high CADD score in con dent regions of a protein model favor pathogenic variants; low ∆∆G Fold and low CADD score favor benign variants and exhibit greater variety in model con dence. Prioritized VUSs have both high ∆∆G Fold and high CADD scores shows a methionine (blue) at position 384, which interacts with three neighboring amino acids (orange). F) Magni cation of Panel E shows three hydrogen bonds between the methionine and neighboring amino acids.
G) The NP_076927.1:p.Met384Lys variant introduces a lysine (blue) in place of methionine, which interacts with four neighboring amino acids, only one of which remains the same as the wildtype interacting neighbors. H) Enlargement of the boxed region from Panel G shows four hydrogen bonds between the lysine and neighboring amino acids. While one hydrogen bond remains the same between the wildtype and variant structures, the NP_076927.