Redundant protein experiments fused into a ranking model for variant effects

Deep mutational scanning has produced enormous catalogs of variant effects, but each experiment speaks its own language. Assays differ in readout, dynamic range, cellular context, and scoring conventions, so effect magnitudes are often not directly comparable even when two studies test the same variants. A team at Molecular Systems Biology now reports that this partial redundancy, usually treated as noise, can be turned into a resource. By fusing overlapping experiments through reciprocal rank fusion, they derived a consensus ordinal score they call variant soundness, spanning more than two million missense variants, and used it to train ESMRank, a sequence-based learning-to-rank predictor that outperforms widely used tools on stability, abundance, and folding benchmarks.

The starting point was a survey of the MAVEdb repository, whose November 2023 snapshot contained 1122 score sets covering 2,176,052 mutations across 596 proteins, predominantly human. More than half of these datasets, 581 score sets and 1,143,050 mutations, partially or completely overlap in the variants they assay. Rather than averaging raw scores, which fails when assays overlap only partially, the researchers represented each experiment as a node in a graph, connected experiments sharing at least 20 percent of assayed variants, and fused the resulting rank lists. Variants that consistently appear among the most deleterious or most tolerated across experiments are reinforced, while assay-specific scaling differences are suppressed.

The integrated landscape showed pronounced biophysical structure. Buried residues were substantially less tolerant of substitution than surface-exposed positions, and nonsense mutations were uniformly deleterious across environments. Hydrophobic-to-polar or charged substitutions were strongly deleterious in protein cores, consistent with stability constraints, whereas surface residues displayed selective sensitivity to charged and polar changes, consistent with disruption of interaction interfaces. Tolerance increased with solvent exposure and decreased with side-chain volume perturbation, with large volume changes markedly more damaging at buried residues. The authors note that the prominence of these stability-associated gradients across heterogeneous assays is consistent with folding defects acting as a shared bottleneck across diverse functional readouts.

Because the consensus signal is inherently relative within proteins, the team formulated variant effect prediction as a learning-to-rank problem rather than regression. ESMRank was trained on 1,033,002 soundness-normalized variants using protein-level stratified cross-validation to prevent information leakage across proteins and experiments. The model employs LambdaMART, a gradient-boosted decision tree implementation of pairwise learning-to-rank, to optimize discrimination between more and less deleterious substitutions within each protein. Variants are encoded with 18 descriptors: three deep features derived from the ESM2 protein language model, including embedding-space differences, attention-derived contact information, and masked marginal residue-probability shifts, plus 15 shallow biophysical, structural, and positional descriptors such as melting temperature and instability index.

Benchmarking began with the Human Domainome dataset, which quantifies effects of 563,534 missense variants across 522 protein domains using an abundance protein fragment complementation assay. After excluding domains sharing more than 80 percent sequence identity with any training protein, ESMRank achieved a median Spearman correlation of 0.62, substantially outperforming ThermoMPNN at 0.46 and other widely used stability and fitness predictors. Performance remained stable under increasingly stringent homology filtering below 50 and 25 percent identity, and gains were consistent across conserved and variable positions, buried and surface residues, and zinc-finger domains, a class in which stability prediction is particularly challenging due to metal coordination requirements.

On curated ProteinGym stability assays in the zero-shot setting, after excluding proteins present in training and restricting evaluation to proteins sharing less than 80 percent identity with any training protein, ESMRank achieved the highest mean Spearman correlation among all competing approaches, at 0.63 across 14 independent stability assays. It ranked first in four of six stratified evaluation subsets based on conservation, disorder, and residue exposure, including variable and intrinsically disordered regions where many predictors degrade. However, the advantage did not generalize uniformly: in activity, binding, expression, and organismal-fitness assays, ESMRank showed lower relative performance and did not consistently approach the top-ranked methods. The authors interpret this as evidence that the model captures a stability-enriched axis of mutational tolerance rather than a general predictor of all functional readouts.

Independent kinetic measurements provided further support. Across 1273 mutations in 39 proteins from the VariBench folding dataset, ESMRank scores correlated with both folding rates, with a median Spearman correlation of 0.52, and unfolding rates, with a median of 0.59, outperforming ThermoMPNN, FoldX, and other predictors. Importantly, ESMRank outputs a within-protein ordinal measure of mutational tolerance, not a calibrated thermodynamic quantity such as ΔΔG or ΔTm, so these benchmarks evaluate rank agreement with folding phenotypes rather than direct thermodynamic prediction.

The clinical relevance of the learned signal was examined without any disease supervision. Within the Human Domainome, which includes 1786 clinically annotated missense substitutions across 114 protein domains, ESMRank produced a sharper separation between pathogenic and putatively benign population variants from gnomAD than ΔΔG-based approaches. In receiver-operating-characteristic analysis, the pathogenicity-specialized predictor AlphaMissense achieved the highest overall performance at an AUC of 0.92, followed by the ESM-1v language-model baseline at 0.87; among stability-oriented predictors, ESMRank was strongest at 0.78, exceeding ThermoMPNN at 0.61, FoldX at 0.63, and the experimental abundance assay itself at 0.58. In an expanded analysis spanning roughly 12,000 proteins and about 95,000 ClinVar missense variants alongside around 8000 common gnomAD variants, ESMRank maintained separation between pathogenic and benign variants at both buried and exposed positions, whereas the discriminatory power of ΔΔG-based predictors decreased substantially at exposed sites.

Aggregating scores per protein also revealed a mechanistic gradient. Analyzing 8081 pathogenic ClinVar missense variants across 836 genes classified into gain-of-function, dominant-negative, autosomal recessive, and haploinsufficiency categories, the team found that genes associated with gain-of-function disease exhibited the highest overall mutational tolerance, followed by dominant-negative and autosomal recessive genes, while haploinsufficient genes were the most constrained. This stratification was most pronounced among buried residues and emerged despite the absence of inheritance or phenotypic information during training, suggesting that constraints relevant to disease mechanisms are implicitly encoded in experimentally derived mutational patterns.

The clearest demonstration came from cystic fibrosis. Mapping per-residue ESMRank scores onto the CFTR structure revealed marked intolerance to substitution within transmembrane helices and interdomain interfaces, including regions with only moderate evolutionary conservation. Against a comprehensive theratype screen of 549 CFTR missense variants, ESMRank correlated strongly with basal transepithelial conductance, with a Spearman correlation of 0.65, and with maturation efficiency measured by the ratio of fully glycosylated to immature CFTR, at 0.558. It also correlated with pharmacological responsiveness: the ivacaftor potentiator effect at 0.52 and the dual corrector response to elexacaftor plus tezacaftor at 0.40. Concordant trends appeared in an independent dataset of 129 missense variants, where ESMRank correlated with baseline CFTR abundance and with corrector responses. Variants responsive to both modulators showed the highest median ESMRank scores, while low or non-responsive variants showed the lowest. In ROC analysis for classifying therapeutically responsive variants, ESMRank achieved an AUC of 0.83, on par with its underlying ESM2 model and ahead of AlphaMissense, ThermoMPNN, and a solvent-accessibility baseline. Notably, ESMRank was not trained on CFTR-specific data or clinical labels, indicating that clinically relevant structure-function relationships can emerge from general sequence-derived constraints. The authors caution that the model does not capture higher-order epistatic effects, does not explicitly incorporate dynamic conformational changes or ligand-dependent states, and that generalization to remote structural classes will depend on continued diversification of experimental data. Whether similar relationships extend to proteins dominated by non-stability-mediated mechanisms remains an open question.

Subject of Research: Computational prediction of missense variant effects by integrating overlapping deep mutational scanning datasets into a learning-to-rank model

Article Title: ESMRank: learning-to-rank prediction of missense variant effects from deep mutational scanning data

Article References: Arnese, R., & Gambardella, G. (2026). ESMRank: learning-to-rank prediction of missense variant effects from deep mutational scanning data. Molecular Systems Biology. https://doi.org/10.1038/s44320-026-00252-3

Image Credits: AI Generated

DOI: 10.1038/s44320-026-00252-3

Keywords: deep mutational scanning, variant soundness, ESMRank, learning-to-rank, protein stability, missense variants, MAVEdb, ProteinGym, CFTR, ClinVar, protein language model, pharmacological responsiveness