Toward Advanced Quantitative Prediction of Missense Variant Effects

1. Introduction

In their influential review, Shendure et al. [1] proposed that advancing genomic medicine requires “pivoting back to basics”—focusing on solving the genotype-to-phenotype equation. Central to their proposed roadmap (their Figure 4) is the development of algorithms to predict the functional impact of genetic variants, a prerequisite for both genetic diagnosis and pharmacogenomics. Because disentangling the full genotype-to-phenotype map—from DNA to disease—remains out of reach, a productive intermediate step is to target the molecular endophenotype: a quantitative, biochemical description of how sequence changes alter protein function [2]. This approach assumes that accurately mapping this molecular layer is a prerequisite for ultimately resolving clinical phenotypes.

The field has increasingly pursued this through data-driven machine learning, a trajectory fueled by the growing availability of large-scale functional data such as deep mutational scanning (DMS) assays [3,4]. This methodological progression has moved from early regression models [5,6] and generative evolutionary models like EVE [7], to protein language models [8,9] and structure-informed predictors like AlphaMissense [10], with the most recent efforts extending to genome-wide language models [11]. While this has expanded predictive scale, it has come at a cost. As models have grown more powerful, they have become opaque black boxes, resisting the mechanistic inspection required for clinical adoption [12].

Furthermore, translating these continuous molecular endophenotypes into clinical utility faces a distinct challenge. Clinical genetics still relies heavily on binary pathogenic/benign classifications, and the binary classifiers used today (e.g., REVEL [13], BayesDel [14]) have become established enough to serve as evidence under ACMG/AMP guidelines [15,16,17]. Yet their true predictive capacity is hard to establish: because these tools are widely used to generate the very clinical labels against which they are later benchmarked, a degree of circularity inflates and obscures their apparent performance [18]. Continuous, mechanistically grounded models offer an alternative that is biologically more realistic and conceptually easier to generalize to other biological layers. Their utility for clinical genetics, however, hinges on a prerequisite: demonstrating that these molecular scores can bridge the gap to clinical reality.

Here, we address this by proposing a framework to model the molecular endophenotype through a fundamental redesign of the Quantitative Assessment of Functional Impact (QAFI) architecture. Building on QAFI 1.0 [19], we reformulate variant-effect prediction as a two-stage problem (Figure 1). First, we estimate the mutation tolerance of a position—defined as the median functional effect across all possible substitutions—reflecting whether the local structural environment is inherently fragile or robust. Second, we model the deviation of individual variants from that positional baseline. This decomposition preserves interpretability at the biological level, even when utilizing non-linear learning via XGBoost.

This work is structured around two goals. First, we establish QAFIsplit as a competitive, mechanistically grounded predictor of the molecular endophenotype. Second, we ask how far our endophenotype predictions carry into clinical classification — a test of the genotype-to-phenotype conjecture at the accuracy we currently reach.

In this work, we demonstrate that:

(i) The two-stage QAFIsplit reformulation provides a biophysically grounded architecture that outperforms general-purpose predictors on independent functional benchmarks—surpassing AlphaMissense in 86% of tested protein domains.

(ii) Carried into clinical variant classification, this molecular-level performance yields an AUC of 0.90 (on par with AlphaMissense), nearly reaching the limit set by its feature representation despite never being trained on clinical data, and it integrates into high-performing meta-predictors without loss of performance, enabling hybrid strategies with a stronger biophysical grounding.

(iii) Under fully blind conditions, both the stand-alone predictor and the meta-predictor strategies showed strong performance across four independent CAGI 7 challenges (ATP7B, LPL, ARSA, and TSC2) — supporting their generalization to entirely unseen proteins.

2. Results

In the following sections, we first describe in more detail our new prediction tool, and then we go through the different training and validation steps followed to establish its validity.

Table S1 and Table S2 provide a concise overview of all the models described and the relationships among them.

• A two-stage framework for quantitative variant effect prediction

In our new model of the impact of missense variants on protein function, we have given a central role to mutation tolerance. Operationally, this property is defined as the median functional effect across all 19 possible amino acid substitutions at a specific position in the protein sequence. It reflects the inherent sensitivity of a position to mutation, and its predicted value may be a relevant component for estimating the impact of specific mutations.

At the technical level, this view translates into a reformulation of the prediction problem into two distinct steps (Figure 1): first, estimating the mutation tolerance of a protein position, and second, using this information, combined with other features, to predict the impact of individual variants at that position.

Both steps employ the same ensemble learning strategy, where predictions from multiple protein-specific models (PSPs) are aggregated to produce the final output. The nature of these PSPs, trained using protein-specific DMS data from the literature, varies depending on the predictive target (mutation tolerance or variant effect). The resulting methodology keeps the name of our original predictor (QAFI, here referred to as QAFI 1.0) [19], although a ‘split’ suffix has been added to the name, QAFIsplit, to indicate the two main predictive steps.

The technical implementation of this design has yielded three versions of QAFIsplit—QAFIsplit 1.0, QAFIsplit 2.0, and QAFIsplit 3.0—which differ in how we solved the variant effect step. A graphical overview of the development process is shown in Figure S1.

In the following sections, we first describe the cross-validated performance of the protein-specific models (PSPs) for each of the two predictive tasks separately — mutation tolerance (Section 2.1) and variant effect (Section 2.2). We then describe the ensemble strategy and its performance as a unified framework (Section 3), followed by an extensive external validation and benchmarking of the resulting models (Section 4).

• The building blocks of our methodology: PSPs

In the following subsections, we describe how PSPs for mutation tolerance (mtPSPs) and variant effect (vePSPs) prediction performed. PSPs form the foundation of the QAFI framework: their within-protein performance defines the upper bound for the ensemble models described in subsequent sections.

2.1. Predicting Mutation Tolerance (mtPSPs)

mtPSP are protein-specific models that estimate positional mutation tolerance and are trained on median effects derived from deep mutational scanning (DMS) data. Table S1 summarizes the key developmental decisions.

Across the 60 proteins, predicted mutation tolerances showed good agreement with experimental values (Figure 2, Figure S2). Visual inspection of observed-versus-predicted plots (Figure 2A) shows that the data tend to align with the diagonal, although with protein-dependent patterns. For example, some mtPSPs tended to overestimate tolerance near the lower end of the observed range (Figure 2A).

The observed–predicted Pearson (Figure 2B) and Spearman (Figure S3) correlations confirm the visual analysis, providing a quantitative assessment of performance variability across proteins. A clear trend emerges: this variability is lower-bounded, as for more than two-thirds of the proteins both correlation measures exceed 0.4, with several reaching above 0.8. Together, these results support the predictive value of the selected features and the robustness of the mutation-tolerance models.

To illustrate performance beyond summary metrics, we examine how it reproduces observed position-specific patterns along the protein sequence, using PRKN as a high-performing example (Figure S4; r = 0.85). Predicted tolerances closely track the observed variations across the sequence, largely independent of annotated domain boundaries. Notably, our method accurately captures the plateau of high experimental mutation tolerance corresponding to a known disordered region in PRKN [20]. This behavior is most likely due to the inclusion, among the model’s predictive features, of AlphaFold pLDDT [21], whose low values are associated with protein disorder [22].

A comparable correspondence between predicted and experimental mutation tolerance profiles is observed across additional proteins in our dataset, with its extent depending on the correlation level (Figure S5). Overall, the model captures global trends and, in higher-correlation cases, also local variations in experimental positional medians along the protein sequence. Mutation tolerance peaks are usually well located, although their predicted size tends to be smaller than their actual size. Together, these results indicate that our predictions preserve meaningful position-specific signals across a broad range of proteins, extending beyond individual high-performing examples.

2.2. Predicting the Variant Effect (vePSPs)

The vePSP models estimate variant functional impact using two PSP model versions, PSPsplit 1.0 and PSPsplit 2.0, which share common input data and feature sets but differ in their training targets. As a reference, we use two baseline models, PSP 1.0 and PSP 2.0. The former corresponds to our original QAFI model, and the second is an upgrade using the full set of predictive features presented here. A summary of the training details and model specifications for all these tools is provided in Table S1.

Model performance was evaluated by comparing observed and predicted variant-level functional scores, illustrated first using VKORC1 as an example (Figure 3A). In this case, we see how the split-based models shift the main cluster of variants closer to the diagonal relative to the baseline approaches, resulting in higher correlation coefficients. For the set of 60 proteins, these trends were summarized using Pearson and Spearman correlations together with prediction errors (MAE and RMSE) (Figure 3B-3C, Figure S6).

Examination of these results reveals several clear trends. First, among the baseline models, PSP2.0 consistently outperformed PSP1.0, demonstrating the added value of the expanded feature set and the non-linear regression model (Figure 3B). Turning to the split-based approaches, PSPsplit1.0 achieved improved performance over PSP2.0 for most proteins, confirming that incorporating mutation-tolerance estimates can enhance variant-level prediction. By contrast, PSPsplit2.0 did not show a systematic advantage: for most proteins, both PSP2.0 and PSPsplit1.0 performed better, indicating that the PSPsplit2.0 strategy did not confer consistent improvement.

All models reduced prediction errors relative to PSP 1.0, with PSP 2.0, PSPsplit 1.0, and PSPsplit 2.0 showing comparable absolute error levels, as further confirmed by MAE and RMSE analyses across all 60 proteins (Figure S6). Together, these results indicate that the main performance gain derives from incorporating mutation tolerance as an intermediate feature.

2.3. Proteome-Wide Quantitative Variant Effect Prediction: The QAFI Ensemble Framework

It is a priori unclear whether the observed PSP-level gains translate to the ensemble level: PSPs may improve by capturing protein-specific effects that are not transferable across proteins, resulting in performance loss upon aggregation. In this section, we examine whether and how these gains map to QAFI ensemble performance.

2.3.1. Defining the QAFI Ensemble Framework

The framework rests on two choices: how many protein-specific predictors (PSPs) the ensemble contains (its size, N) and which proteins they correspond to. Both are determined jointly by a nested leave-one-protein-out procedure that selects, for a given N, the PSPs generalizing most consistently across proteins (Figure S7A; full procedure in Methods).

We estimated performance by cross-validation, holding out each protein in turn and predicting its variants from an ensemble built only from the others — so that every protein is scored as if it were new. Across the 60 proteins, we tested five ensemble sizes (N = 10, 15, 30, 45, and 60) and selected the one performing best in most proteins, N = 30 (Figure S7B-S7C); the corresponding proteins are listed in Table S2.

2.3.2. Cross-Validated Performance of the QAFI Ensemble Framework

The ensemble strategy described in Section 3.1 was applied separately to the different PSP families (see Table S3 for the mapping between QAFI and PSP models), resulting in three main QAFI models: QAFI 2.0 (based on PSP 2.0), QAFIsplit 1.0 (based on PSPsplit 1.0), and QAFIsplit 2.0 (based on PSPsplit 2.0). In addition, we introduced a fourth model, QAFIsplit 3.0, derived from PSPsplit 2.0 but using a correlation-weighted fusion scheme (Figure S1B) instead of the median for aggregation. For each PSP–QAFI pair (e.g. PSPsplit 1.0 vs QAFIsplit 1.0), we compared model performance separately for the 60 proteins in our dataset, to assess the performance deviation introduced by moving from within-protein to cross-protein prediction and to estimate the generalization capacity of the ensemble framework.

In all cases, data points clustered close to, but predominantly above, the diagonal (Figure 4A, Figure S8), indicating that ensemble models largely preserve the performance trends of their underlying PSPs, although a generalization gap is observed. SNCA (α-synuclein) emerged as a consistent outlier showing a marked PSP–QAFI discrepancy across models, consistent with its intrinsically disordered nature [23], which contrasts with the predominantly globular character of the proteins whose PSPs define the ensemble.

Figure 4. Cross-validated performance of QAFI ensemble models. (A) Per-protein Pearson correlations of PSP 2.0 (within-protein, y-axis) versus QAFIsplit 3.0 (cross-protein ensemble, x-axis) across the 60 proteins. Each point represents one protein, colored by its PSP 2.0 Pearson correlation (color scale, right). The dashed diagonal indicates equal performance; points above the diagonal reflect the within-protein advantage of PSP-level models over the ensemble. SNCA (α-synuclein) is highlighted as a consistent outlier; the inset shows its three-dimensional structure, illustrating its intrinsically disordered character, which likely accounts for its poor generalization under the ensemble framework. (B) Radial plot of per-protein Pearson correlations for all five QAFI ensemble models across the 60 proteins. Protein names are arranged along the perimeter; concentric rings indicate correlation values from 0.25 (inner) to 0.85 (outer). The shaded area corresponds to QAFI 1.0 (baseline); colored lines represent QAFI 2.0 (blue), QAFIsplit 1.0 (pink), QAFIsplit 2.0 (purple), and QAFIsplit 3.0 (orange).

Figure 4. Cross-validated performance of QAFI ensemble models. (A) Per-protein Pearson correlations of PSP 2.0 (within-protein, y-axis) versus QAFIsplit 3.0 (cross-protein ensemble, x-axis) across the 60 proteins. Each point represents one protein, colored by its PSP 2.0 Pearson correlation (color scale, right). The dashed diagonal indicates equal performance; points above the diagonal reflect the within-protein advantage of PSP-level models over the ensemble. SNCA (α-synuclein) is highlighted as a consistent outlier; the inset shows its three-dimensional structure, illustrating its intrinsically disordered character, which likely accounts for its poor generalization under the ensemble framework. (B) Radial plot of per-protein Pearson correlations for all five QAFI ensemble models across the 60 proteins. Protein names are arranged along the perimeter; concentric rings indicate correlation values from 0.25 (inner) to 0.85 (outer). The shaded area corresponds to QAFI 1.0 (baseline); colored lines represent QAFI 2.0 (blue), QAFIsplit 1.0 (pink), QAFIsplit 2.0 (purple), and QAFIsplit 3.0 (orange).

When compared with the original QAFI 1.0, all new ensemble models showed an overall performance gain (Figure 4B, Figure S9), demonstrating that the improvements observed at the PSP level translated effectively to the ensemble setting. Differences among QAFIsplit models and QAFI 2.0 were generally minor, with no model showing a clear systematic advantage over the others.

2.4. External Validation and Benchmarking of QAFI Models

The cross-validated results presented in the previous section provide an internal assessment of model performance. Here we go further, evaluating the QAFI framework on entirely independent data not involved in any stage of model development. This external validation comprises three complementary analyses: a large-scale benchmark across 402 human protein domains, an evaluation on clinically annotated variants, and performance assessment in four CAGI 7 challenges (TSC2, ATP7B, ARSA, and LPL). Each analysis also serves as an opportunity to benchmark QAFI against state-of-the-art predictors in the field.

2.4.1. Large-Scale Validation Across 402 Human Protein Domains

The large-scale site-saturation mutagenesis dataset recently published by Beltrán et al. [24] provides experimental measurements across an independent collection of human protein domains. The size and technical homogeneity of this dataset, together with the diversity of protein domains covered, make it an important validation resource.

After removing any potential overlap with the development set, the benchmark comprised 492,096 missense substitutions across 402 protein domains. Using this benchmark to compare QAFI family members, we observe that QAFIsplit 3.0 outperforms the remaining members of the QAFI family in the majority of domains, ranking above QAFI 1.0 (87%), QAFI 2.0 (71%), QAFIsplit 1.0 (67%), and QAFIsplit 2.0 (61%) of cases (Figure S10).

We next benchmarked QAFIsplit 3.0 against a representative panel of external predictors (Figure 5). Overall, QAFIsplit 3.0 outperformed all general-purpose methods and achieved comparable performance levels, without any adaptation or modification, to the specialized versions of EVE and ESM1v, which were applied taking into account the characteristics of the domains in Beltran et al.’s experiment. For thermoMPNN, given the substantial overlap between its training data and the Beltrán benchmark, the only dedicated stability predictor in the panel, we restricted the comparison to the homologue-free subset used by the original authors. In this setting, QAFIsplit 3.0 matched thermoMPNN’s performance (51% of domains).

An interesting feature of QAFIsplit 3.0 predictions is that they tend to compensate for poor performance of other methods: there is a tail extending towards the right of the distribution, corresponding to proteins where the alternative method performed poorly while QAFIsplit 3.0 did not. This feature supports the idea that, alongside its potential as a stand-alone predictor, QAFIsplit 3.0 could also be used to complement existing in silico tools.

The variant-level success of QAFIsplit 3.0 is underpinned by the ensemble estimation of positional mutation tolerance (QAFImt). Figure S11 illustrates the performance of this first stage across the independent Beltrán et al. dataset. The radial plot (right panel) maps the per-protein Pearson correlations between predicted and observed positional medians for all 402 domains, revealing a heterogeneous landscape. Less than 20% of domains exhibit poor agreement (r < 0.4), and more than 50% achieve good to excellent accuracy (r > 0.6). The representative sequence profiles (left panels) exemplify domains from this upper tier, demonstrating that when the baseline is well-estimated, QAFImt successfully captures the major trends and local variations in biophysical constraints. While the model struggles with certain protein contexts, this distribution confirms that the split architecture is nonetheless driven by a broadly reliable positional baseline for the majority of the proteome.

2.4.2. Benchmarking QAFI on Clinically Annotated Variants

Here we ask how far our endophenotype models — our route to the genotype-to-phenotype problem — carry, at their present accuracy, into the clinic. As a test, we evaluated QAFIsplit 3.0 against a spectrum of established predictors on a common set of 12,787 clinically annotated variants (Figure 6).

To interpret QAFIsplit 3.0’s clinical coherence, we organized these methods into three reference groups. The first consists of established tools developed for clinical variant interpretation, many of which have been calibrated for use under ACMG/AMP guidelines (see Methods). This clinical ecosystem spans a remarkably broad performance range — from GERP++ (0.68), MPC (0.70) and PrimateAI (0.74) at the lower end, through PhyloP and FATHMM (0.77–0.78), SIFT (0.81), and CADD and PolyPhen (0.87), up to MutPred2 (0.93), REVEL and VEST4 (0.94), ClinPred (0.95) and BayesDel (0.96) at the top. Against this frame, QAFIsplit 3.0 reached an AUC of 0.90: without ever seeing a clinical label, it outperforms most of these tools and sits well above the median of the ecosystem.

The second reference indicates how much clinically relevant signal our feature representation can capture when optimized for classification. We estimate it with QAFIclas (0.92), an XGBoost classifier built on the same features as QAFIsplit 3.0 but trained directly on clinical labels. QAFIsplit 3.0 approaches this reference even though it never sees a clinical label.

The third group consists of advanced machine-learning methods — such as CPT-1 (0.93), AlphaMissense (0.91) and EVE (0.89) — which, like QAFI, are not trained on clinical labels, though some have since been calibrated for clinical use [17]. QAFIsplit 3.0 falls within the range of this group, which delimits what can currently be achieved by predicting variant impact from molecular data alone.

Taken together, these comparisons place the endophenotype route within reach of the tools built directly on clinical data. Applicability was further supported by combining QAFIsplit 3.0 with AlphaMissense, popEVE, CPT-1, REVEL and BayesDel by median aggregation (QAFImeta), which gave an AUC of 0.95 against 0.94 for the same combination without it: a continuous, biophysically grounded score integrates into established meta-predictors without penalty, forming the basis for the hybrid strategies we fielded in CAGI 7 (Section 4.3).

2.4.3. External Blind Validation Through the CAGI 7 Community Experiment

To extend the independent validation of our methodology, we participated in the 7th edition of the Critical Assessment of Genome Interpretation (CAGI 7) community experiment (December 2025). CAGI provides an unbiased evaluation framework in which participants submit predictions without access to experimental outcomes, subsequently assessed by independent expert panels [26]. We focused on four protein-centric challenges involving missense variants measured through high-throughput functional or stability assays: ATP7B (Wilson disease, n=15,795), LPL (familial lipoprotein lipase deficiency, n=19,138), ARSA (metachromatic leukodystrophy, n=8,867), and TSC2 (tuberous sclerosis complex, n=8,891). Since our methodology is designed exclusively for missense variants, predictions were restricted to this subset in all four challenges. For all these challenges, we took advantage of CAGI guidelines that allow participants to submit up to six predictors, identifying one of them as their method of choice.

At the time of writing, the final challenge rankings have not been formally published; however, communications from the CAGI 7 organizers confirmed that our preferred models consistently ranked among the leading methods. Pending formal publication of these results, we describe here how we selected, among members of the QAFI family, our six candidates per challenge. The selection process was primarily evidence-based but with a qualitative integration step, as limited protein-specific data were available for fully automated ranking. It combined the performance of our models on the Beltrán et al. [24] dataset with a ROC analysis based on ClinVar-annotated variants for each challenge protein (Figure 7). The six models selected per challenge are shown in Figure 7, and a brief description of each is provided in Table S3.

Pre-submission evaluation of our models on ClinVar-annotated variants for the four challenge proteins confirmed strong discriminative performance across all challenges (Figure 7). AUC values were particularly high for LPL (0.978–0.998) and ATP7B (0.904–0.931), with solid performance also for ARSA (0.871–0.899) and TSC2 (0.868–0.922). These results are consistent with the clinical applicability demonstrated in Section 4.2, now extended to four fully independent proteins not involved in any stage of model development.

3. Discussion

Responding to the call to resolve the genotype-to-phenotype equation, a prerequisite for genomic medicine [1], we framed variant-effect prediction as the modeling of a molecular endophenotype. This framing established two objectives for our work: first, to develop a method capable of accurately resolving this molecular layer; and second, to ask how far the resulting predictions carry into clinical classification. The QAFIsplit framework addresses the first objective. By reformulating prediction as a two-stage decomposition — first estimating the mutation tolerance of a position, and then modeling the specific deviation introduced by a substitution — QAFIsplit yields a biophysically grounded, mechanistically intuitive predictor. Crucially, this architectural choice does not sacrifice quantitative accuracy; as demonstrated by our benchmarks, it achieves state-of-the-art quantitative performance on independent functional assays, surpassing black-box alternatives in the vast majority of protein domains.

The success of this two-stage decomposition lies in its alignment with the biophysical reality of proteins. A variant’s functional impact is not generated in a vacuum; it is heavily conditioned by the local structural and energetic environment of the residue. Our results show that explicitly modeling this environment as “mutation tolerance” provides a highly informative prior for variant-level prediction, in accordance with Yu et al. [27]. The ability of the mtPSPs to recapitulate position-specific tolerance profiles—such as accurately identifying disordered regions in PRKN via AlphaFold pLDDT scores—confirms that our expanded feature set effectively captures the determinants of positional constraint. By separating this first-order effect from the second-order perturbation caused by a specific amino acid change, the framework reduces the burden on the downstream regressor, leading to more robust predictions.

This mechanistic grounding translates into tangible benchmark advantages. On the large-scale Beltrán et al. [24] dataset, QAFIsplit 3.0 outperformed AlphaMissense in 86% of protein domains and matched the performance of specialized, assay-tailored versions of protein language models like ESM-1v. We attribute this to the explicit inclusion of structurally informed descriptors—such as the normalized contact energy (NCE-NR) and evolutionary sequence context (neco features)—which capture local physicochemical constraints that deep-learning representations, whether structure-based or LLM-derived, may implicitly encode but cannot explicitly leverage for out-of-distribution generalization. Notably, QAFIsplit 3.0 exhibited a “compensatory tail,” maintaining high accuracy on domains where other predictors failed. This suggests that the QAFI framework captures complementary aspects of variant biophysics, making it particularly valuable as a robust, stand-alone tool for newly characterized domains.

Translating quantitative functional scores into clinical utility presents a distinct challenge, as clinical genetics still relies heavily on binary pathogenic/benign classifications. Here the endophenotype route proves efficient, as postulated by Masica and Karchin [2]. Trained exclusively on continuous DMS data and never exposed to a clinical label, QAFIsplit 3.0 reached an AUC of 0.90 — outperforming most tools in this clinical ecosystem, and sitting well above the median of that ecosystem, and close to the reference set by its own feature representation (0.92 for QAFIclas, trained directly on clinical labels). Part of what separates it from the upper end of that range (0.94–0.96) is the current accuracy of endophenotype prediction, and that part should narrow as this resolution improves. Part, however, may not: clinical labels rest on evidence that lies outside the molecular layer — segregation and trio studies, clinical observation, and, for some tools, population allele frequencies — which a model confined to the endophenotype neither uses nor aims to. Consistent with this, adding QAFIsplit 3.0 to a meta-predictor of established tools (QAFImeta) left their performance intact (0.94 to 0.95), showing that a continuous, biophysically grounded score can join the clinical toolkit without diluting it.

This behavior was further examined prospectively, under blind conditions, in the CAGI 7 community experiment. QAFI-based models — both the stand-alone predictor and meta-predictor strategies — were submitted to four independent challenges on proteins entirely absent from model development (ATP7B, LPL, ARSA, and TSC2), where their strong pre-submission validation on clinical variants (Section 4.3) was borne out in a fully blind setting. Communications from the organizers indicated that our models were among the better-performing entries, with formal rankings to appear in the CAGI 7 special collection. This extends the value of continuous, mechanistically grounded scores — alone or combined with existing classifiers — to genuinely unseen prediction settings.

An interesting technical observation emerges from the internal validation: while incorporating mutation tolerance as a feature (PSPsplit 1.0) consistently improved predictions, explicitly modeling the residual deviation from that positional value (PSPsplit 2.0) did not. We hypothesize that predicting residuals amplifies noise in the DMS median estimates, whereas feeding the predicted tolerance directly into XGBoost allows the non-linear model to learn the interaction between that positional estimate and the physicochemical properties of the new amino acid. This highlights an important balance in interpretable machine learning: explicit biological decomposition is powerful, but it must be coupled with algorithms that can handle the inherent noise of experimental biological data.

Despite these advances, limitations remain. The generalization gap observed between protein-specific predictors (PSPs) and the ensemble framework indicates that transferring local mutational landscapes across proteins is non-trivial. The failure to accurately predict SNCA (α-synuclein) underscores this: intrinsically disordered proteins violate the globular packing assumptions that underpin many of our structural features. Furthermore, the reliance on AlphaFold models and high-quality multiple sequence alignments restricts proteome coverage for the split models to approximately 80% of variants. Future iterations of QAFI could address these limitations by integrating disorder-specific predictors or graph neural networks capable of modeling long-range allosteric effects and protein complexes.

In conclusion, the QAFI family of predictors demonstrates that high-resolution quantitative prediction and mechanistic interpretability are synergistic, not mutually exclusive. As the field moves beyond binary variant classification, frameworks like QAFIsplit will be essential for translating the full complexity of the genotype-to-phenotype map into actionable medical knowledge.

4. Materials and Methods

4.1. DMS Variant Datasets

We compiled 60 deep mutational scanning (DMS) datasets to train the protein-specific predictors (PSPs) that form the core of QAFI. Thirty were inherited from our previous work [19], and thirty were newly selected from the ProteinGym collection [28] using the following criteria: (i) the assay contained more than 1,000 single missense variants; (ii) when several assays existed for the same protein, the one with broader positional coverage was kept; (iii) mutation coverage exceeded 50% of residue positions within the assayed region; (iv) an AlphaFold model was available [21,29]; and (v) a multiple sequence alignment (MSA) could be built. These criteria ensure sufficient statistical power for regression modeling and the availability of complementary structural and evolutionary features. The full protein list and metadata are given in Table S4.

4.2. Score Normalization

Because raw DMS scores are assay-specific and not directly comparable across experiments, all datasets were standardized using the normalization scheme of Gray et al. [5], as in the original QAFI predictor [5]. Under this scheme, a score of 1 corresponds to variants with no detectable impact on protein function, scores below 1 indicate decreased activity, and scores above 1 enhanced activity relative to the wild type. Normalized variant-level scores for all 60 datasets are provided in Table S5.

4.3. Functional Validation Dataset

QAFI models were validated in the large-scale mutagenesis dataset of Beltran et al. [24]. Entries without valid UniProt identifiers or AlphaFold structures were removed, retaining only missense variants, and excluding proteins with residue-numbering mismatches between the UniProt sequence and the AlphaFold model. The final set contained 402 proteins and 492,096 missense variants.

4.4. Dataset of Clinically Annotated Variants

We assembled a set of clinically annotated missense variants from the January 2025 releases of ClinVar [30] and Humsavar [31]. Variants labeled Pathogenic, Likely pathogenic, or Pathogenic/Likely pathogenic were grouped as pathogenic, and those labeled Benign, Likely benign, or Benign/Likely benign as benign. We retained only variants with a review status of at least one star, excluding variants of uncertain significance and those with conflicting classifications. Variants deposited before 2021 were excluded to minimize overlap with the training sets of the predictors evaluated here. All proteins present in the DMS training set were excluded to avoid circularity. The resulting set comprised 42,173 variants across 2,012 human proteins (13,167 pathogenic, 29,006 benign).

This collection is used here to compare predictor performance. Because coverage differs across methods, all cross-method comparisons were carried out on the subset scored by every one of them: 12,787 variants (6,606 pathogenic, 6,181 benign).

4.5. Protein Structure Dataset

To parameterize the structure-based features, we curated a non-redundant set of high-resolution structures from PISCES [32] (X-ray structures at ≤50% sequence identity, resolution ≤1.8 Å, R-factor ≤0.25). For each protein we kept only PDB chains with a one-to-one UniProt match and 100% sequence identity after alignment, using chain A when multiple chains were present. Structures were processed with ARPEGGIO [33] to compute atomic contacts, and entries with duplicate residue indices or insufficient MSA depth were excluded. After filtering, 4,355 structures with valid contact annotations were retained.

4.6. Multiple Sequence Alignments

MSAs were built following Riera et al. [34]. For each target protein, homologs were retrieved from UniRef100 [35] with PsiBlast [36] (E-value 0.001, two iterations); sequences with less than 40% identity to the target were discarded; and the remaining sequences, together with the query, were aligned with Muscle [37].

4.7. Feature Representation

Each variant was represented by sequence- and structure-based features. Fourteen features were taken from our previous work [19], spanning evolutionary, physicochemical, and structural information (e.g., Blosum62 scores, Shannon entropy, PSSM, pLDDT, solvent accessibility, and Miyazawa–Jernigan contact-energy terms). In this study we extended the representation with 14 new descriptors of the native residue and its environment, parameterized from the curated set of 4,355 structures. They are organized into four conceptual groups—an energy-based descriptor (NCE-NR), sequence-context descriptors (the neco family), probability-density descriptors (PDFFs), and mutation tolerance—described below. Full definitions for every feature are listed in Table S6.

4.8. Normalized Contact Energy of the Native Residue (NCE-NR)

NCE-NR is a coarse-grained measure of the contact energy of the wild-type residue within its structural environment. Using the Miyazawa–Jernigan (MJ) potential [38], we identify the spatial neighbors of the native residue (side-chain centroids within 6.5 Å, or Cα for glycine) and sum their MJ interaction energies to obtain $E_{\mathrm{native}}$. This is normalized by the maximum contact energy observed for that amino acid type across the curated dataset,

$$\mathrm{NCE}-\mathrm{NR}={\displaystyle \frac{E_{\mathrm{native}}}{E_{\mathrm{native}\_\max }}}$$

(1)

which removes intrinsic contact biases across amino acid types and yields a relative measure of local structural stability.

4.9. Evolutionary Compatibility with Neighbor Sequence Context (neco)

To capture context-dependent substitution effects, we introduced the neco family of features, which quantify how well a substitution fits its surrounding sequence environment. Unlike position-wise conservation scores, neco features model dependencies between the target residue and conserved neighboring residues using log-odds scores derived from MSAs. We defined four variants—neco1, neco2, neco3, and neco12—representing sequence context at different distances.

The immediate-neighbor score, neco1, evaluates compatibility with the first neighbors at positions $i-1$ and $i+1$. For a native residue $A_i$ flanked by $X_1$ and $Y_1$ and substituted by $B$,

$$\mathrm{neco}1=\ln \left({\displaystyle \frac{P(X_1\hspace{0.17em}B\hspace{0.17em}{Y}_1\hspace{0.17em}|\hspace{0.17em}\mathrm{human}=X_1{A}_i{Y}_1)}{P(X_1\hspace{0.17em}\mathrm{*}\hspace{0.17em}{Y}_1\hspace{0.17em}|\hspace{0.17em}\mathrm{human}=X_1{A}_i{Y}_1)}}\right)$$

(2)

where the numerator is the probability of observing the mutant residue $B$ under the conserved flanking context and the denominator ($\mathrm{*}$ = any central residue) reflects the background frequency of that context in the MSA. The remaining variants follow the same log-odds form but use more distant or combined neighbors: neco2 and neco3 use the second ($i\pm 2$) and third ($i\pm 3$) neighbors, respectively, and neco12 jointly uses the first and second neighbors on both sides ($i\pm 1,i\pm 2$). Each neco feature is a 20×20 matrix in which higher values indicate substitutions more compatible with the local evolutionary context; their explicit forms are given in Table S7.

4.10. Probability Density Function Features (PDFF)

Building on previously defined residue-level descriptors (e.g., pLDDT, contact energy, conservation entropy), PDFFs quantify how the empirical distribution of a feature shifts when a native residue $A$ is replaced by $W$, relative to the background at all positions carrying $A$. For a substitution $A\to W$ and a feature value $x$.

$$r\left(x\right)=\ln \left({\displaystyle \frac{p(x\hspace{0.17em}|\hspace{0.17em}A\to W)}{p(x\hspace{0.17em}|\hspace{0.17em}A)}}\right)$$

(3)

where $p(x\hspace{0.17em}|\hspace{0.17em}A\to W)$ is the density of the feature among homologs in which $A$ is replaced by $W$, and $p(x\hspace{0.17em}|\hspace{0.17em}A)$ the background density at all positions carrying $A$, regardless of mutation status. Because empirical densities can be noisy in rare substitution contexts, the numerator and denominator distributions were each fitted with three-component Gaussian mixture models (and, for a small number of features with atypical distributions, alternative functional forms) before computing $r\left(x\right)$. The final PDFF of each feature is a 20×20 matrix of scoring functions $r\left(x\right)$ that can be applied to any substitution in any protein.

4.11. Mutation Tolerance

To capture the overall sensitivity of a residue position to mutational effects, we explicitly model mutation tolerance, operationally defined as the median DMS score observed at each position; constrained positions show low tolerance, whereas permissive positions show high tolerance. Unlike the other descriptors, mutation tolerance is not a direct physical quantity but is predicted from local features. For each protein, the per-position median DMS score is paired with a vector of position-level features $x_p$ (PSSM, Shannon entropy, entropy of sequence neighbors, pLDDT, pLDDT bin, Colasi, fraction of conserved 3D neighbors, FANC, FBNC, and NCE-NR), and a multiple linear regression (MLR) is fitted across positions,

$${\tilde{y}}_p={\beta }_0+x_p^{\mathrm{\top }}\mathsf{\beta }+{\epsilon }_p$$

(4)

where ${\beta }_0$ is the intercept, $\mathsf{\beta }$ the regression coefficients, and ${\epsilon }_p$ the error term. The fitted model predicts the mutation tolerance at any position from its feature vector; this value is the input used by the split models below.

4.12. Feature Standardization

All features were standardized to zero mean and unit variance using z-score standardization,

$$x\mathrm{\text{'}}={\displaystyle \frac{x-{\mu }_{\mathrm{train}}}{{\sigma }_{\mathrm{train}}}}$$

(5)

where ${\mu }_{\mathrm{train}}$ and ${\sigma }_{\mathrm{train}}$ are computed on the training set only and then frozen and reused for the validation and test data, preventing information leakage from the held-out sets.

4.13. Baseline Variant Effect Models: PSP 1.0 and PSP 2.0

To quantify the improvements introduced by the two-stage design, we defined two baselines at the protein-specific (PSP) level. PSP 1.0 corresponds to the original QAFI protein-specific predictor [19]: a separate multiple-linear-regression model for each of 30 proteins, using 14 sequence- and structure-based features. PSP 2.0 is an updated version that extends both the feature set and the number of proteins — a separate model for each of the full set of 60 proteins, using the 28 features described above — and replaces multiple linear regression with XGBoost regression. Both predict variant effects directly from the full feature set, without explicitly separating position-level from mutation-specific effects.

4.14. Two-Step Position-Informed Regression

Because the functional impact of a mutation is strongly shaped by the structural and evolutionary constraints of its position, we developed a two-step regression framework that separates position-level effects from mutation-specific deviations. In both variants, the first step estimates the mutation tolerance ${\tilde{y}}_p$ at each position as described above; the second step uses this estimate, together with the mutation-level feature vector $x_i$, to predict the variant effect.

In PSPsplit 1.0, the predicted mutation tolerance is added as an extra input feature to an XGBoost regressor,

$${\widehat{y}}_i=f_{\mathrm{XGB}}\left(\right[x_i,\hspace{0.17em}{\tilde{y}}_p\left]\right)$$

(6)

In PSPsplit 2.0, the model instead predicts the residual between the observed score and the predicted tolerance,

$$r_i=y_i-{\tilde{y}}_p, {\widehat{r}}_i=f_{\mathrm{XGB}}\left(\right[x_i,\hspace{0.17em}{\tilde{y}}_p\left]\right), {\widehat{y}}_i={\tilde{y}}_p+{\widehat{r}}_i$$

(7)

where $y_i$ is the observed (normalized) DMS score of variant $i$, $r_i$ the residual relative to the predicted tolerance ${\tilde{y}}_p$, ${\widehat{r}}_i$ the residual predicted by the XGBoost regressor $f_{\mathrm{XGB}}$ from the mutation-level features $x_i$, and ${\widehat{y}}_i$ the final predicted variant effect. This formulation isolates the mutation-specific perturbation from the predicted tolerance. Both frameworks were evaluated by leave-one-position-out cross-validation within each protein [39], with performance assessed by Pearson ($r$) and Spearman ($\rho$) correlation, mean absolute error, and root mean squared error. All models were implemented with the XGBoost version in scikit-learn [40], using hyperparameters fixed across proteins.

4.15. Ensemble Architecture and Aggregation

The protein-specific predictors (PSPs) are the building blocks of the general QAFI models: a variant is scored by aggregating predictions from many PSPs, so that any variant in any protein can be predicted even when no DMS data exist for that protein. All QAFIsplit variants share the same three-step pipeline:

Step 1 — Training. A PSP is trained independently on each of the 60 proteins.

Step 2 — Selection. A fixed subset of N = 30 PSPs — those generalizing best across proteins — is selected using the nested leave-one-protein-out procedure described below (Ensemble construction and size selection).

Step 3 — Aggregation. The 30 selected PSPs are combined into a single prediction. The choice of PSP type and fusion rule defines the QAFIsplit versions, where ${\widehat{y}}_{i,m}$ denotes the prediction of model $m$ for variant $i$.

QAFIsplit 1.0 and QAFIsplit 2.0 (median fusion). These aggregate the top-30 PSPsplit 1.0 and PSPsplit 2.0 models, respectively, giving every model equal influence through the median,

$${\widehat{y}}_i^{\mathrm{median}}={\mathrm{median}}_{m=1,\dots ,30}\hspace{0.17em}{\widehat{y}}_{i,m}$$

(8)

This is robust but treats all 30 models as equally informative, ignoring that they differ in how reliable their predictions are.

QAFIsplit 3.0 (correlation-weighted fusion). This version uses the same 30 PSPsplit 2.0 candidates as QAFIsplit 2.0 but replaces the median by a weighted average that favors models agreeing with the ensemble consensus. Each candidate is applied to the validation proteins and the pairwise Pearson correlations among the 30 predictions are computed; the reliability of model $m$ is its mean correlation with the other 29 candidates,

$${\overline{r}}_m={\displaystyle \frac{1}{29}}\sum _{j\ne m}{C}_{mj}$$

(9)

where $C$ is the correlation matrix. The reliabilities are normalized into weights and used to combine predictions,

$$w_m={\displaystyle \frac{{\overline{r}}_m}{{\displaystyle \sum _{j=1}^{30}}{\overline{r}}_j}}, {\widehat{y}}_i^{\mathrm{fusion}}={\displaystyle \sum _{m=1}^{30}}{w}_m\hspace{0.17em}{\widehat{y}}_{i,m}$$

(10)

Candidates systematically out of step with the others are thus down-weighted, while consistent models contribute more strongly.

Baseline ensemble models. Two ensemble baselines frame these results. QAFI 2.0 applies the same pipeline to the PSP 2.0 predictors (top-30 selection, median aggregation), so that its gain over QAFI 1.0 reflects the expanded features and XGBoost regression alone, not the two-stage split. QAFI 1.0 is the published ensemble [19] — the median of its ten original PSP 1.0 predictors.

4.16. Size and Composition of QAFI Ensemble Models

A key feature of the QAFI framework is its ensemble strategy, in which predictions from many protein-specific models are aggregated into a single, generalizable variant-effect score. Two choices define the ensemble: its size (the number of PSPs, N) and the identity of its members. Both were determined with a nested leave-one-protein-out procedure (Figure S7A).

In the outer loop, each protein k is held out in turn (k = 1…60) and its PSP excluded. In the inner loop, each of the remaining 59 PSPs, i, is applied to the DMS data of the other 58 proteins (excluding both k and i’s own protein), yielding 58 Pearson correlations whose median is taken as the performance of PSP i. The 59 PSPs are then ranked by this median, and the top N are recorded. Repeating the outer loop over all 60 proteins yields 60 such ranked lists, and the N PSPs appearing most frequently across them define the deployed ensemble.

The same procedure provides an unbiased performance estimate: for each held-out protein k, the N PSPs selected from the others are applied as an ensemble to its DMS data, and the Pearson correlation between predicted and observed scores is computed. This yields 60 cross-validated correlations across the full dataset. The ensemble size was chosen by repeating the whole procedure for N = 10, 15, 30, 45, and 60 and, for each protein, identifying the N with the highest correlation; the value winning in most proteins, N = 30, was selected (Figure S7B–S7C).

4.17. Clinical Variant Classifier (QAFIclas)

As a reference for the classification performance attainable from our feature representation, we trained QAFIclas, an XGBoost classifier that uses the same features as QAFIsplit 3.0 but is trained directly on clinical labels rather than on continuous DMS scores. QAFIclas is intended only to indicate how much clinically relevant signal these features can capture when optimized for classification; it is not proposed as a variant classifier for practical use.

QAFIclas was fitted on the curated pathogenic/benign set described above and evaluated on the same common set of 12,787 variants used throughout Section 4.2, reaching an AUC of 0.92. The model used reg_alpha = 10, reg_lambda = 10, max_depth = 4, learning_rate = 0.1, subsample = 0.9, colsample_bytree = 0.8, n_estimators = 200, and log-loss as the evaluation metric; its output ranges from 0 (pathogenic) to 1 (benign).

4.18. Functional Benchmarking

QAFI was benchmarked against existing variant-effect predictors on the Beltran et al. dataset [24] described above. The dataset provides scores from twelve representative predictors—ESM1v_domain [9], RaSP [41], ddMut [42], FoldX [43], popEVE [44], EVE [7], Tranception [45], EVE_domain [7], rsASA (this is a residue accessible surface area computation), ThermoMPNN [46], AlphaMissense [10], and ESM1v_full-length [9]—used for comparison in quantitative functional prediction.

Whereas the predictors above were all evaluated under identical conditions, ThermoMPNN [46] required a separate treatment because it was trained on the mega-scale dataset of Tsuboyama et al. [47], which overlaps substantially with the Beltrán et al. benchmark [24]. To avoid the data leakage that would occur if it were evaluated on domains from its training set, we restricted its comparison to the homologue-free subset defined by Beltrán et al. [24]— benchmark domains with no Pfam homology to the Tsuboyama et al. training set. All other predictors were compared on the full 402-domain benchmark.

4.19. Clinical Benchmarking

QAFIsplit 3.0 was benchmarked against a panel of established variant-effect predictors on the curated pathogenic/benign set described above. Most predictor scores were retrieved through VEP (release 113) [48], which at the time of analysis drew on the dbNSFP v4.7c plugin [49]; scores for popEVE and CPT-1, which are not available through VEP, were obtained directly from their respective official sources (accessed January 2026). All predictors were evaluated on the common subset of variants scored by every method (n = 12,787), and scores were oriented, where necessary, so that higher values consistently indicate greater predicted pathogenicity. Discriminative performance was quantified as the area under the ROC curve (AUC).

For interpretation, the predictors were organized into three reference groups: (1) tools developed for clinical variant interpretation—several of which have been specifically calibrated for ACMG/AMP guideline evidence [16]. To provide a comprehensive view of the clinical ecosystem, this group includes a representative set of such classifiers available in the dbNSFP v4.7c plugin (e.g., SIFT [50], PolyPhen [51], CADD [52], GERP++ [53], PhyloP [54], FATHMM [55], MPC [56], PrimateAI [57], MutPred2 [58], REVEL [13], VEST4 [59], ClinPred [60], and BayesDel [14]); (2) QAFIclas, trained directly on clinical labels as a theoretical reference for the performance attainable from the feature representation (see above); and (3) advanced machine-learning methods not trained on clinical labels (e.g., AlphaMissense [10], EVE [7], CPT-1 [61], popEVE [44]). We further assessed whether QAFIsplit 3.0 contributes complementary signal by combining it with AlphaMissense, popEVE, CPT-1, REVEL and BayesDel through median aggregation (QAFImeta), evaluated on the same common subset.

4.20. Writing Support

The text of this article was reviewed for grammar and clarity using an AI language model. All suggested corrections were carefully checked to ensure they faithfully reflect our intended meaning and results.