Characteristics of predicated pathogenic missense variants in human genes for Mendelian disorders

Missense variants represent a major class of genetic variation and contribute substantially to the cause of Mendelian disorders. However, the genetic and protein properties governing their pathogenicity remain incompletely characterized. Existing variant pathogenicity databases, such as ClinVar, are enriched for variants in well-studied genes and therefore may not provide a representative view of pathogenic missense variants across the genome. Meanwhile, recent advances in computational pathogenicity prediction have significantly improved their predictive accuracy, enabling large-scale genome-wide analysis. In this study, we investigated the genetic and protein properties associated with the predicted pathogenicity of missense variants in Mendelian disease genes. We collected all possible missense variants in Mendelian genes curated from the Online Mendelian Inheritance in Man (OMIM). Each variant was annotated with pathogenicity predictions from gMVP and a range of properties, including evolutionary conservation, physicochemical substitution characteristics, and functional and structural genetic and protein features. Using multiple linear regression models, we investigated the relationships between the various properties and the predicted pathogenicity of the missense variants. Interaction effects among those properties were also examined. We first validated the quality of gMVP by comparing the effects of each property in the linear regression, fitting to either gMVP pathogenicity predictions or ClinVar pathogenicity labels, on the ClinVar dataset. Then, we extended the analyses to the entire missense dataset to fit gMVP pathogenicity scores. In addition, gene set and robustness analyses were further performed to validate the study findings. In 655,859 missense variants from ClinVar, 7.1% were labeled as pathogenic/likely pathogenic compared with 6.4% predicted as pathogenic using a gMVP threshold of 0.9. Regression analyses revealed a high degree of concordance between ClinVar and gMVP pathogenicity annotations, with most predictors showing associations in the same direction and of similar magnitude. The genome-wide dataset includes 23,339,549 missense variants in 4,709 Mendelian disease genes, with 31.8% predicted to be pathogenic using a default gMVP threshold of 0.75. Evolutionary conservation, solvent accessibility, and intrinsic disorder were strongly associated with gMVP pathogenicity scores. Other features, including functional annotations and amino acid biochemical properties, also contributed to the predicted pathogenicity with moderate effect sizes. Feature interactions support an integrative model in which conservation and local protein structure jointly influence the variant impact. GO-stratified and other gene-set analyses indicated that these associations were broadly conserved across gene functional contexts and inheritance patterns. gMVP-predicted pathogenic missense variants are consistently enriched in conserved and structured protein core regions, regardless of gene function and inheritance pattern. These findings highlight the importance of incorporating genetic and protein annotations into variant interpretation frameworks.

Missense variants, which result in single amino acid substitutions, represent one of the most common classes of coding genetic variation in the human genome (Chong et al., 2015; Cooper and Shendure, 2011). While most missense variants are benign, a small subset disrupts protein function and leads to Mendelian diseases, making the accurate identification of the pathogenic variants a central challenge in human genetics and precision medicine (Richards et al., 2015). Despite substantial advances in variant pathogenicity prediction, interpreting the factors underlying these predictions remains a major challenge (Starita et al., 2017; Grimm et al., 2015).

Current pathogenicity prediction methods leverage diverse sources of information, including evolutionary conservation, population allele frequencies (AF), protein sequence and structural information, and predictions from other computational methods (Ng and Henikoff, 2001; Adzhubei et al., 2010; Kircher et al., 2014; Sundaram et al., 2018). Recent methods, including EVE (Frazer et al., 2021), AlphaMissense (Cheng et al., 2023), Evo 2 (Brixi et al., 2026), and protein language model-based predictors (Brandes et al., 2023), have demonstrated improved predictive performance across diverse datasets. However, they often operate as black boxes, producing accurate predictions without providing clear biological interpretation. As a result, the extent to which specific protein properties drive pathogenicity predictions remains unclear. This gap limits our understanding of how missense variants exert their effects and constrains the interpretability of computational models.

While prior studies have demonstrated that pathogenic missense variants tend to be enriched in conserved or structurally constrained regions (Samocha et al., 2017; MacGowan et al., 2024), these analyses have typically focused on individual features or limited feature sets. Protein domains, binding interfaces, regulatory motifs, post-translational modification (PTMs), and 3D structures are all critical for protein stability, molecular interactions, and biological activities (Chothia, 1992; Pawson and Nash, 2003; Neduva and Russell, 2005; Walsh et al., 2005). Variants affecting these critical regions are more likely to disrupt protein function than variants occurring in structurally flexible or functionally less important regions (Yue et al., 2006). In addition, amino acid physicochemical properties, such as hydropathy, charge, molecular weight, and solvent accessibility, may further modulate the functional impact of missense substitutions (Dill, 1990; Rost and Sander, 1994; Tokuriki and Tawfik, 2009; Stein et al., 2019). Systematic, genome-scale quantification of how multiple genetic and protein structural and functional features jointly shape missense pathogenicity across Mendelian disease genes remains limited. Moreover, interaction effects between these properties have rarely been explicitly modeled, and the consistency of these protein determinants across different functional contexts has not been comprehensively evaluated.

Here, we address these gaps by performing a comprehensive, genome-wide analysis of missense variants across Mendelian disease genes. We compiled all possible missense variants in the human reference genome and annotated them with pathogenicity predictions, genetic, and protein annotations. By jointly modeling these features, explicitly testing their interaction effects, and evaluating consistency across gene subgroups and random subsets of variants, we provide a systematic framework for understanding the genetic and protein determinants of missense variant pathogenicity.

The analysis pipeline consisted of variant collection, pathogenicity annotation, genetic and protein feature mapping, and statistical modeling (Figure 1).

Datasets

Mendelian disease genes were curated from the Online Mendelian Inheritance in Man (OMIM) database (Amberger et al., 2015). We annotated the pathogenicity of all possible missense variants in these genes using the Graphical Missense Variant Pathogenicity Predictor (gMVP) (Zhang et al., 2022). gMVP was developed using graph attention neural networks to capture both evolutionary and co-evolutionary features and is among the best-performing predictors (Livesey and Marsh, 2022). We used its rank scores for correlation and regression analyses, and its binary classifications to estimate the percentage of predicted pathogenic variants and to compare it with ClinVar pathogenicity labels. All variants across the entire genome were ranked and normalized to [0, 1], with one indicating the most likely pathogenicity. We used the default gMVP rank score greater than 0.75 to estimate the percentage of pathogenic variants. A set of amino acid-level annotations was also collected for each missense variant, based on the canonical transcript in UniProt (The Ensembl Variant Effect Predictor).

Physicochemical properties of both reference and alternate amino acids were annotated based on the amino acid physical property table from the Thermo Fisher Scientific Pierce Protein Methods resource (Kyte and Doolittle, 1982), including hydropathy (Kyte and Doolittle, 1982), charge (Grantham, 1974), and molecular weight (Zimmerman et al., 1968). Hydropathy was categorized as hydrophilic, moderately hydrophobic, or hydrophobic, and charge as positive, neutral, or negative. Hydropathy, charge, and weight changes were calculated as the difference between the alternate and reference amino acids and were used in the downstream analyses. Amino acid substitution similarity was measured using the BLOSUM62 substitution score (Henikoff and Henikoff, 1992).

Evolutionary conservation was assessed using Eigen scores (Davydov et al., 2010; Ionita-Laza et al., 2016). Population AFs were obtained from gnomAD3 (maximum population AF) and the UK Biobank (Chen et al., 2024; Bycroft et al., 2018). Missing population AFs were imputed as 10 –6. Negative log-scaled AFs were calculated to represent selective constraint and population-level rarity in the analyses.

Protein structural features, including secondary structure and solvent accessibility (Rost and Sander, 1994; Dictionary of protein secondary structure, 1990), were predicted by DSSP (Kabsch and Sander, 1983) based on AlphaFold2-predicted protein 3D structures. For secondary structure, the DSSP-produced eight categories were grouped into three: helix (α-helix, 3-helix, and 5-helix), strand (isolated β-bridge and extended β-strand participating in a β-ladder), and other conformations (hydrogen-bonded turn, bend, and unassigned residues). These grouped secondary-structure annotations were then one-hot encoded. During the regression analyses, the group of other conformations was excluded from the model because it was perfectly collinear with the linear combination of the helix and strand groups, resulting in a singular design matrix. Solvent accessibility was normalized into relative solvent accessibility (RSA) using the same procedure as previous studies (Matveev et al., 2023) within each protein using min-max scaling. Missing solvent accessibility values were imputed using the mean value within the protein.

Protein functional and other structural annotations were obtained from UniProt and included protein domains, binding sites, PTMs, sequence motifs, and intrinsically disordered residues (The UniProt Consortium, 2025). All annotations were manually curated, except for intrinsic disorder, which was predicted using MobiDB (Necci et al., 2017), an ensemble method that combines outputs from multiple individual predictors. All variants were annotated with the above features, which were normalized by their absolute maximum value. This unified dataset served as the basis for our main statistical analyses.

Quality check

We validated the quality of gMVP predictions by comparing feature contributions to both clinically curated and gMVP-predicted pathogenicity on the ClinVar dataset, including variants classified as pathogenic/likely pathogenic or benign/likely benign (Landrum et al., 2016). To match the clinical variant classification standard, which requires a precision greater than 0.9 for a likely pathogenic variant, we used a gMVP rank score greater than 0.9 as the cutoff for comparison with ClinVar pathogenicity labels.

Statistical analysis

Correlation analysis

gMVP rank scores were utilized in all downstream analyses. Spearman’s correlation analysis was used to assess monotonic associations among gMVP rank scores and various genetic and protein features (Haigh and Conover, 1981; Mukaka, 2012). This approach also identifies pairwise correlation among the features we collected at the genome scale.

Regression analysis

Regression analysis (Gelman and Hill, 2007) was used to model missense variant pathogenicity using genetic and protein features and to estimate their contributions. In addition, we incorporate interaction terms between all pairs of features to capture potential non-additive effects. To assess gene functional context-specific effects, separate regression analyses, with and without interaction terms, were performed within three GO categories (The Gene Ontology Consortium, 2026): Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). To determine whether the findings generalized based on inheritance patterns, we repeated the gMVP-based analyses in genes with dominant inheritance. The inheritance patterns were obtained from the OMIM database. Given the large number of variants in the genome-wide analysis, we additionally evaluated the robustness of the regression results in subsets of variants. We repeated the regression analyses 10 times using randomly sampled 10,000 variants, summarizing the mean and standard deviation of the regression coefficients and the consistency of statistical significance (P value <0.05) across the 10 repeats.

Dataset

Distribution of gMVP scores

After applying reference version control, transcript harmonization, and variant-type filters, a total of 23,339,549 missense variants across 4,709 Mendelian disease genes in OMIM were collected. The distribution of gMVP rank scores is shown in Figure 2A. The median gMVP rank score was 0.584 (interquartile range, 0.337–0.809), with a mean of 0.566 (standard deviation (SD) = 0.277) and 31.8% as predicted pathogenic variants with gMVP rank scores above 0.75. The proportion of predicted pathogenic variants varied substantially across genes, from nearly 0% to nearly 100%, as shown in Figure 2B. Most genes contained a relatively small proportion of predicted pathogenic variants, whereas only a minority of genes exhibited a high predicted pathogenic variant burden, resulting in a unimodal distribution with a pronounced right tail. gMVP pathogenicity predictions based on a rank score above 0.75 are expected to classify 25% of missense variants as pathogenic. However, many genes exhibited a moderate to high proportion of predicted pathogenic variants, indicating that many Mendelian genes are conserved and that their variants are more likely to be pathogenic. The high percentage of predicted pathogenic variants does not imply a correspondingly high burden of disease-causing variants. In recessive diseases, the disease occurs only when both alleles carry pathogenic variants.

To determine whether this gene-level heterogeneity of pathogenic variant burden differed by inheritance pattern, we next compared proportions of predicted pathogenic variants between autosomal dominant (AD) and autosomal recessive (AR) genes (Figure 2C). This analysis included 1,172 AD genes and 3,537 AR genes. AD genes showed a slightly higher proportion of predicted pathogenic variants than AR genes, with a mean of 35.8% compared with 32.5%. Despite this moderate shift, the AD and AR distributions were highly similar, indicating that the inheritance pattern alone does not fully explain gene-level differences in predicted pathogenic variant burden. Instead, the broad variability within both inheritance groups suggests that additional genetic and protein properties likely contribute to differences in susceptibility to missense variants across Mendelian disease genes.

Association between pathogenic variant proportion and protein size

We examined the relationship between the proportion of predicted pathogenic variants and protein size. The distribution is shown in Figure 3A. The exponential decay model estimated a negative slope of −0.027 per 100 amino acids. Similarly, we observed a weak negative correlation between the proportion of predicted pathogenic variants and the total number of missense variants per protein (Spearman ρ = −0.126) in Figure 3B. The small effect size indicates that protein length or variant count does not account for the substantial heterogeneity of predicted pathogenic variant burden observed across genes. Therefore, the proportion of predicted pathogenic variants likely reflects intrinsic gene-level properties.

Proportion of pathogenic variants across gene Ontology categories

To evaluate whether predicted pathogenic variant burden differs across broad functional contexts, we compared the portion of predicted pathogenic variants in the three primary GO categories: BP, CC, and MF. The number of proteins and missense variants represented in each GO term-specific subset is summarized in Supplementary Table S1. Figure 4 shows modest variability in the proportion of predicted pathogenic variants across GO categories. These findings suggest that gene-level predicted pathogenic variant burden is not primarily determined by broad GO functional categories, but rather by finer-grained structural and functional features that shape the landscape of missense pathogenicity.

Consistency between ClinVar and gMVP-predicted pathogenicity

The ClinVar dataset includes 655,859 unique missense variants across 2,862 genes, of which 46,549 (7.1%) are pathogenic/likely pathogenic. gMVP predicts 42,194 variants (6.4%) as pathogenic with a rank score greater than 0.9, comparable to the fraction identified by ClinVar labels.

Before investigating feature associations with gMVP pathogenicity at the genome scale, we first assessed whether gMVP pathogenicity was concordant with clinically curated pathogenicity from ClinVar. Figure 5 shows that the genetic and protein features were consistently correlated with ClinVar and gMVP pathogenicity in the ClinVar dataset. Eigen scores, solvent accessibility, motif, and binding sites were consistently associated with pathogenicity defined by the two methods. Other functional annotations, including BLOSUM scores, secondary structure, AF, and changes in biophysical properties, had little effect on the pathogenicity. These analyses indicate that gMVP provides robust pathogenicity predictions.

Pairwise correlations among gMVP pathogenicity and genetic and protein features

Figure 6A shows the Spearman’s correlation for all pairwise monotonic relationships among the gMVP rank score and various annotated features. Among features correlated with gMVP rank score, Eigen scores showed the strongest positive association (ρ = 0.403), whereas RSA showed the strongest negative association (ρ = −0.396). Intrinsic disorder was also negatively correlated with gMVP rank score (ρ = −0.278). This pattern was not detected in the ClinVar dataset because the number of strictly defined pathogenic variants in intrinsically disordered regions was too small to provide sufficient statistical power. Supplementary Figure S1 shows that variants within and outside intrinsically disordered regions have distinct distributions of gMVP rank scores. Other features correlated with gMVP rank scores include protein domain (ρ = 0.161), strand (ρ = 0.157), and helix (ρ = 0.106).

Previous studies have shown that the AlphaFold2 prediction confidence measure, pLDDT, was associated with conformationally flexible regions (Piovesan et al., 2022). Our study also shows a correlation between intrinsic disorder and unstructured secondary structures (ρ = 0.364). Intrinsic disorder was positively correlated with RSA (ρ = 0.319). RSA was negatively correlated with Eigen scores (ρ = −0.289), strand (ρ = −0.293), helix (ρ = −0.196), and domain (ρ = −0.191), indicating that residues buried inside proteins were more often found in constrained, structured, and functional domains.

Associations between genetic and protein features and predicted missense variant pathogenicity using linear regression

To collectively characterize feature associations with the predicted pathogenicity while accounting for other variables, we performed linear regression using all features, including evolutionary conservation, structural, functional, AF, and physicochemical predictors. A forest plot for feature contributions was shown in Figure 6B. After accounting for all other variables, Eigen scores still showed the strongest positive association with the gMVP rank score (β = 0.601), and RSA showed the strongest negative association (β = −0.267), consistent with our correlation analysis. Intrinsic disorder remained negatively associated with predicted pathogenicity (β = −0.143). These three major associations were robust in the subsampling analysis, remaining statistically significant in all 10 independent runs of 10,000 variants (Supplementary Table S2). Several additional features were also associated with gMVP rank scores with smaller effect sizes, including motif (β = 0.071), PTM (β = 0.055), molecular weight change (β = −0.051), and protein domain (β = 0.040). Analyses using gMVP rank scores at the genome level yielded larger effect sizes than those using gMVP binary predictions in the ClinVar dataset, highlighting the power of quantitative measures and a larger dataset. Together, these results indicate that predicted missense pathogenicity is associated with multiple layers of genetic and protein information, including evolutionary conservation, residue exposure, local structure, and functional annotations.

Contributions of feature interactions to missense variant pathogenicity

To systematically evaluate feature associations and interaction effects, we constructed regression models that included interaction terms among all genetic and protein features. Supplementary Table S3 shows that Eigen scores remained the strongest positive predictor (β = 1.051), indicating that conservation is a dominant predictor of pathogenicity even in the presence of feature interaction terms. Several interaction terms showed strong associations with the gMVP rank score, including the interactions between intrinsic disorder and RSA (β = 0.308), which remained statistically significant in all 10 independent subsampling runs of 10,000 variants (Supplementary Table S3). Even though both intrinsic disorder and RSA are negatively correlated with pathogenicity, this positive interaction term suggests that a variant in the structured protein core region does not increase pathogenicity. The same applies to Eigen scores and RSA (β = 0.239); this interaction term counteracts the additive effects of conserved and buried protein regions. Together, these results indicate that predicted missense pathogenicity is shaped not only by strong associations among individual features but also by interactions among structural, functional, evolutionary conservation, and physicochemical features.

Subgroup analyses across gene Ontology categories

To assess whether the determinants of missense pathogenicity depend on gene functional context, we conducted regression analyses stratified by GO categories across MF (Supplementary Table S4), BP (Supplementary Table S5), and CC (Supplementary Table S6). The regression models fitted for each gene subset included the same covariates as the genome-wide interaction model and were used to evaluate whether the major feature associations were preserved across GO functional categories.

Eigen scores remained the strongest and most consistent positive predictor of gMVP rank score, consistent with the genome-wide models. Intrinsic disorder is also a stable determinant of predicted missense pathogenicity regardless of GO categories. In addition, several interaction terms were consistently observed across GO-stratified models. The disorder-solvent accessibility interaction was positive across the MF, BP, and CC subsets. Together, these results indicate that the major determinants of predicted missense pathogenicity are broadly preserved across GO functional contexts.

Targeted subset analyses

Having established that the evolutionary and structural signatures were preserved across GO functional categories, we next examined whether these findings also held true within other gene subsets.

Subgroup analysis in intrinsically disordered genes

Intrinsic disorder and the disorder-solvent accessibility interaction were among the major structural associations identified in the genome-wide regression analyses. We next evaluated whether these effects depended on the inclusion of genes lacking annotated disordered regions. Genes lacking disorder annotations may influence estimates of disorder-associated effects in the genome dataset. To address this possibility, we repeated the analysis in a disorder-restricted dataset comprising 2,189 proteins containing at least one annotated disordered region, thereby excluding proteins that are entirely structured.

Figure 7A shows that intrinsic disorder remained a strong negative predictor (β = −0.123, P -value <10 –300) in the dataset containing only disordered genes. In addition, other predictors show effects similar to those observed in the genome dataset. These results confirm the robustness of our findings. The interaction model also indicated that the disorder-solvent accessibility interaction remained strongly positive (β = 0.316, P -value <10 –300, shown in Supplementary Table S7). The consistency of both the single-variable effects and interaction terms across the genome and disorder-restricted datasets demonstrates that the observed relationships are not driven by completely ordered genes.

Subgroup analyses in dominant genes

To further assess whether the associations observed across the genome dataset were driven by inheritance patterns, we conducted an additional analysis restricted to dominant genes. This subset included 1,173 unique genes and 5,701,973 missense variants. In the dominant-gene regression model (Figure 7B), Eigen scores remained the strongest positively associated feature (β = 0.597, P -value <10 –300). RSA (β = −0.248, P -value <10 –300) and intrinsic disorder (β = −0.149, P -value <10 –300) showed strong negative associations with pathogenicity. The interaction model for the dominant genes further reinforced the relationship between intrinsic disorder and RSA (β = 0.302, P -value <10 –300), shown in Supplementary Table S8. These results indicate that the associations do not depend on the type of gene inheritance.

Case studies

To visualize the relationship between predicted pathogenicity, evolutionary conservation, and protein structural features, we used KCNJ5 as a representative protein (Figure 8A). This protein encodes a G protein-activated inward-rectifying potassium channel subunit (Choi et al., 2011). In KCNJ5, gMVP-predicted pathogenic variants are enriched in the central region of the protein, indicating position-dependent differences in predicted pathogenicity. Evolutionary conservation exhibits a similar pattern, with higher Eigen scores in the central region, suggesting that stronger evolutionary constraint contributes to increased pathogenicity. In contrast, solvent accessibility is negatively correlated with pathogenicity, with the protein center displaying the lowest accessibility scores. Variants located within the intrinsically disordered region consistently have the lowest gMVP rank score across KCNJ5. The disordered region also corresponds to the low-confidence segments in the AlphaFold-predicted 3D structure (Figure 8B), located at both termini of the protein. Overall, the KCNJ5 case study demonstrates that missense pathogenicity is spatially heterogeneous and reflects the combined influence of evolutionary constraint and local protein structural context.

In this study, we performed a genome-scale, multi-predictor analysis of missense variant pathogenicity predictions across Mendelian disease genes, leveraging gMVP alongside various structural, functional, evolutionary, and physicochemical annotations. Across more than 23 million missense variants in 4,709 Mendelian disease genes, our results demonstrate that evolutionary conservation and protein structural context are major determinants of missense pathogenicity.

In both pairwise correlation and regression analyses, the gMVP rank score was strongly associated with evolutionary conservation, intrinsic disorder, and solvent accessibility. Variants located in conserved, structured, and protein-core regions were more likely to be pathogenic, whereas substitutions in unconserved, disordered, or protein-surface regions were relatively tolerated. Among other predictors, PTM was also positively associated with pathogenicity, suggesting that functional elements also contribute to pathogenic burden. Importantly, these associations persisted after adjusting for other predictors, including AF, amino acid substitution properties, and functional annotations. These findings indicate that missense pathogenicity is potentially governed by a combination of evolutionary conservation, structural, and functional context.

Subgroup analyses, including the ClinVar dataset, GO-stratified gene sets, dominant genes, and disordered genes, showed a consistent association between gMVP-predicted missense variant pathogenicity and evolutionary conservation, protein structural measures, and their interactions. Overall, these subset analyses support the robustness and biological relevance of these features as determinants of missense variant pathogenicity.

Several limitations of this study should be acknowledged. First, although gMVP demonstrates strong predictive performance, its predictions should not be interpreted as equivalent to clinically established pathogenicity. Second, the evolutionary conservation overlap conceptually with information incorporated into gMVP, but it is the most predictive feature used by almost all variant pathogenicity methods. Third, several protein features, such as solvent accessibility and intrinsic disorder, were derived from computational predictions and are therefore subject to prediction error. Likewise, protein functional annotations from UniProt, such as domains, PTMs, and binding sites, remain incomplete and may contain omissions or inaccuracies. In the future, more accurate AF estimates can be obtained from newer databases with larger sample sizes. Additional analyses incorporating nonlinear modeling may capture higher-order relationships between missense variant pathogenicity and other factors.

This study provides a quantitative and integrative framework for understanding how gene and protein context shape missense pathogenicity across Mendelian disease genes. By combining correlation, regression, interaction modeling, GO-stratified analyses, gene-subset robustness analyses, and ClinVar-based concordance analyses, we show that missense pathogenicity, as assessed by the gMVP prediction surrogate, is fundamentally shaped by evolutionary conservation and protein structural context. These structural effects are robust across functional categories and consistent between computational predictions and clinically curated datasets. This framework provides a mechanistic basis for interpreting missense variation and for prioritizing potentially disease-causing variants.