Improved pathogenicity prediction for rare human missense variants

1Fully realizing the clinical potential of human genome sequencing will require the ability to accurately determine the pathogenicity of individual sequence variants. Although genome-wide association (GWA) studies have limited power to associate rare missense variants with disease,, rare missense changes cannot be ignored: over 99% of the set of observed missense variants are rare (with a global minor allele frequency [MAF] below 0.5%), and 90% are extremely rare (with MAF < 10−6). Moreover, if genome sequences were known for all humans alive today, each single-nucleotide variant that is compatible with life would appear an average of ∼50 times. Because identifying the subset of rare missense variants that are damaging represents a major unmet challenge for personal genome interpretation, we focus here on rare and extremely rare missense changes.

2It has become possible to experimentally assess the functional impact of nearly all possible missense variants for a target protein, yielding “deep mutational scans” or “variant effect maps.”,,,,,,,,,, However, the MaveDB resource currently contains variant effect maps for fewer than 1% of the ∼4,000 human disease-associated proteins, and a high-quality fully-comprehensive experimental atlas of functional missense variation could be decades away.

3By contrast, computational methods to infer variant pathogenicity can already be applied at genome scale. However, performance of current methods,,,,,,,,,,,,,,,,,,,,,,,,,,, has not fully addressed the challenge of variant interpretation. For example, using thresholds such that 90% of pathogenicity predictions are correct, such methods identified only 10%–20% of pathogenic variants. Current guidelines for clinical variant interpretation recommend that all computational methods be (at best) treated as “weak evidence.”

4To improve computational predictors, we would like to increase the size of the training dataset while also ensuring that the quality of training data is high, i.e., the training examples are accurately labeled and representative of the variants for which pathogenicity inference is most needed (e.g., rare variants). However, there is a tension between these competing objectives. For example, although common variants can be more confidently annotated as truly benign and are therefore likely to offer higher labeling accuracy, there is evidence that common benign variants are not representative of rare benign variants. However, excluding common variants would substantially reduce the number of clinically annotated benign variants available for training computational models (Figure S1). While ClinVar requires clinical interpretations to be based on a rigorous and attributable process, we can be less sure of accuracy for other resources offering pathogenicity annotations. For example, HumsaVar offers annotations solely from literature reports and HGMD seeks to maximize recall while accepting the risk of including false positives. The variants reported in gnomAD have been used as “putatively-benign” training examples by many computational methods,, despite the potential for contamination with pathogenic variants, especially in genes associated with late-onset diseases. MaveDB provides experimental “variant effect maps” in which the functional impact of many variants has been measured in parallel, but the set of scored variants may not be representative of the spectrum of patient variation.

5Here, we describe the VARITY approach for pathogenicity prediction, which has been specifically optimized for rare missense variation. While VARITY uses a meta-prediction strategy, it limits the circularity that can arise in such approaches by excluding any feature that was informed by variant pathogenicity annotation. We also exclude features that may serve as proxies for protein identity because these may lead to predictions that are biased by the fraction of each protein’s variants in the training set that are annotated pathogenic (which may be inaccurate in future application settings). VARITY judiciously harnesses a larger set of training examples with uncertain accuracy and representativity and uses differential weighting strategies to ensure that training set expansion improves performance on a high-quality test set (Figure 1). For the task of identifying rare pathogenic variants, we find that VARITY outperforms all other computational approaches examined.

6To assemble a resource of labeled variants and associated features for human proteins, we first retrieved from dbNSFP (V4.0b2) all missense variants for ∼18,000 human proteins that had a “reviewed” status in UniProt [download date: 2019-04-08]). We next assembled variant features (summarized in Table S1) within four main categories. First, from dbNSFP we collected scores for existing predictors based on only conservation information learned via unsupervised learning methods, e.g., the SIFT, method. To avoid circularity in subsequent performance testing, we excluded as features those predictors that had made direct or indirect use of annotated pathogenic and benign variants in training,,,,,,,, e.g., the PolyPhen2, method. Where multiple scores for a missense variant were available from a given predictor due to different coding isoforms, the score from the canonical isoform (defined in UniProt) was chosen as the final score. Where multiple scores were available for different nucleotide-level changes yielding the same missense variant, the most deleterious score was chosen as the final score. Second, we defined a set of “delta values” by the difference in various quantitative physicochemical properties (e.g., molecular weight, polarity, or charge) between the reference and substituted amino acid. Third, we assembled features related to protein-protein interaction. For example, “maximum buried area”—defined by the largest solvent-accessible surface area around a given amino acid that is buried by any protein interaction partner—or the “maximum solvation energy change” —defined for a given amino acid by the greatest change in solvation energy provided by any interaction partner (both were estimated by PDBePISA). Fourth, we assembled structure-related features. These included a set of binary variables indicating secondary structure (obtained from PSIPRED) and accessible surface area (estimated by PDBePISA). In addition, we used substitution matrix (BLOSUM100) scores and whether the amino acid position is within or outside of a Pfam (download date: 2019-04-08) domain. All 43 features and their sources can be found in Table S1. The resource of features for all missense variants from ∼18,000 human proteins can be downloaded via the VARITY web portal (see web resources).

7For the purpose of training, we first identified a set of ∼4,000 human disease-related proteins from all ∼18,000 human proteins on the basis of whether any variant (not restricted to missense variants) had been annotated as “pathogenic” or “likely pathogenic” in the ClinVar database or “deleterious” in the HumsaVar database (download date: 2019-04-08). We further included a list of human proteins from multiplexed assays of variant effect (MAVE) studies,,,, (see Table S2). Only variants from the set of human disease-related proteins were labeled and considered in model training. We assigned every missense variant from this set a label of either positive (putatively pathogenic), negative (putatively benign), or unknown. At this stage, our confidence in these labels might be quite low, but we will later depend on the machine-learning strategy to weight these variants. The labeling criteria depended on the source of information about the variant. For the ClinVar database, regardless of the number of “review stars,” variants with “pathogenic,” “likely pathogenic,” or “pathogenic/likely pathogenic” annotation were labeled as positive. Similarly, variants with “benign,” “likely benign,” or “benign/likely benign” were labeled as negative for training purposes. All other variants were labeled as unknown. For the HumsaVar database, variants with the “polymorphism” annotation were labeled as negative, variants with “deleterious” annotation were labeled as positive, and all other variants were labeled as unknown. For gnomAD, variants for which at least one homozygous individual had been observed were labeled as negative, and other variants were labeled as unknown. For variants from MAVE studies,,,, (encompassing 12 proteins, see Table S2), we used functional impact scores rescaled by a previously developed pipeline, such that a score of 0 indicates fitness approximating that of a nonsense variant and a score of 1 indicates fitness approximating that of a synonymous variant. Here, variants with a score above 0.5 were labeled as negative, while those below 0.5 were labeled as Positive. For the HGMD database, all variants were labeled as positive. To handle the minority of cases (∼10% of all variants) where a variant appears in multiple databases, we assigned a primary source for that variant by using a somewhat arbitrary ranking of the sources from which to derive the variant label: ClinVar, HGMD, and HumsaVar pathogenicity annotations, fitness scores from variant effect maps, and presence of homozygotes in gnomAD.

8All missense variants labeled as positive or negative were considered as potential training examples. To train and test VARITY models, we assembled a core training set of variants that were of high quality (i.e., accuracy and representativity). We also used various “add-on” training sets that potentially have less predictive utility (e.g., because they are less accurately labeled or less representative of the core set), but these sets were subjected to further weighting. To obtain high-accuracy core sets, we used ClinVar missense variants. To ensure high representativity, we first separated ClinVar missense variants into three subsets: ClinVar_C, ClinVar_MR, and ClinVar_ER with labeled missense variants that were common (MAF > 0.5%), moderately rare (10−6 < MAF < 0.5%), and extremely rare (MAF < 10−6), respectively. Given evidence that common variants are not representative of rare variants, we excluded ClinVar_C from the core training set of the VARITY model optimized for rare missense variants (VARITY_R; see Figure 1) and treated it as an add-on set. The properties of moderately rare and common variation may differ from those of extremely rare variants. Therefore, to optimize a VARITY model (VARITY_ER; see Figure 1) for extremely rare missense variants, we used ClinVar_ER as the core set and used ClinVar_C and ClinVar_MR as add-on sets.

9All labeled missense variants outside of the core set were considered as add-on training examples and were separated into different add-on sets so that each set can be employed separately with distinct effects on training VARITY models. Add-on variants were first grouped into add-on training sets on the basis of source of annotation (the source where the variant label was determined) and label (positive or negative). Each negative-labeled add-on set (e.g., negative-labeled variants from gnomAD) was further split on the basis of rarity (MAF cutoff = 0.5%) except the negative set from MAVE studies in which most scored variants are rare. A full list of core set and add-on sets for both VARITY_R and VARITY_ER models can be found in Table S3.

10To evaluate performance of the VARITY models, we used the standard area under the receiver operating characteristic curve (AUROC). However, we also wished to use measures based on precision and recall, which often correspond more closely to the user's needs and intuition. Unfortunately, any measure using precision, such as the area under the precision recall curve (AUPRC), changes depending on the prior (the frequency of positively labeled examples). One strategy to compare AUPRC measures applied to different test sets with varying priors is to force each test set to be balanced, e.g., by down-sampling the most abundant label class. However, this approach destroys information by eliminating test data. Therefore, we developed two measures: area under the balanced precision recall curve (AUBPRC) and recall at 90% balanced precision (RB90P). For these measures, we simply calculate the precision for each test set and then derive the “balanced precision,” i.e., the precision that would have been expected had the prior been balanced (equal to 50%):

11Proof:

12Let be a random variable representing the true label of an event (e.g., if a variant is pathogenic) and represent the label predicted by a model (e.g., if the model predicts a variant to be pathogenic). Precision and prior of the model can then be written as

13The following equation holds according to Bayes Rule:

14Using Equation 2 and Equation 3, we can rewrite as:and then from Equation 4, we obtain the likelihood ratio

15Because we know the likelihood ratio is not affected by the prior (ratio of conditional probability conditioned on prior), therefore we can write down the following, given a balanced prior of 0.5:

16By combining Equation 5 and Equation 6, we can solve for balanced_precision and obtain Equation 1.

17A balanced precision versus recall curve (BPRC) can be drawn by calculation of the balanced precision value corresponding to every precision value, and AUBPRC can be calculated just as AUPRC would be. We note that there are different ways to calculate an AUPRC (or AUBPRC). Here, we calculated AUPRC (and AUBPRC) as the weighted mean of precision (or balanced precision) at different recall thresholds and used the increase in recall from the previous threshold as the weight, which is also called average precision. It can also be shown that AUBPRC can be calculated more directly from AUPRC via

18To derive RB90P, we need to first plot the BPRC and then locate (or interpolate) the recall value at which balanced precision is 90%.

19For VARITY we adopted a variant of the gradient boosted tree (GBT) machine-learning algorithm XGBoost (Python Package V0.90). Briefly, the first stage of the GBT method is to generate an initial prediction of the target variable for each training example. In successive stages, GBT iteratively generates a series of decision trees. After each tree, the current prediction for any given training example is the initial prediction added to the sum of the outputs for all trees generated thus far. Each tree in the series is optimized to yield incremental outputs that, when added to the sum of previous outputs, tend to minimize the loss function (prediction error). The incremental value is determined by approximating both the first and second order gradient of the loss function at the current prediction stage. Optimization of each tree is accomplished by successively selecting features that split training examples into smaller groups (or nodes). The splitting feature is selected such that a single incremental output value assigned to each subgroup can most closely approximate the set of best incremental output values. For each training set (core or add-on set), each of several quality-informative properties was used as input to a logistic function that assigned a weight to each training variant within the set. The parameters of each logistic function were treated as hyperparameters that were optimized for performance (see details in the Bayesian optimization of hyperparameters section below). Subsequently, the weights for all positive training examples were scaled by a factor such that the total of re-scaled weights for positive training examples was equal to the total weight of negative examples.

20To identify the quality-informative properties used as described above for variant weighting, candidate properties were evaluated by moving window analysis. Candidate informative properties included: allele frequency (for all training sets except the ones from MAVE studies,,,, ), “review stars” (for ClinVar training sets), “number of homozygotes” (for gnomAD training sets), “label confidence,” and “mutational accessibility” (the latter two are candidates for MAVE training sets,,,, ). The “label confidence” score for each variant was defined to be either the MAVE score itself (for scores < 0.5) or 1 − score (for scores > 0.5). Use of the mutational accessibility property was motivated by the fact that, although variants with more than one nucleotide change in a given codon can be generated by mutagenesis strategies used in variant effect mapping, these rarely occur in humans. The mutational accessibility of each possible amino acid change was the total probability of all possible single-nucleotide change events yielding that amino acid change, where the probability that a single nucleotide change event occurs in each codon was defined via the human codon usage frequency for that codon divided by nine (because nine single-nucleotide edits are possible for any given codon). Amino acid changes requiring more than one single-nucleotide change within the same codon received a mutational accessibility score of 0.

21For each proposed candidate property, we carried out moving window analysis to evaluate whether model performance depends on the value of the property. To this end, add-on set examples were first ordered by the property. We then defined NW moving windows, each containing a fraction FE of the total X add-on training examples and incrementing the position of successive windows every I = [X ∗ (1 − FE)/(NW − 1)] training examples. Where the number of examples is not evenly divisible by I, the final window may contain slightly fewer than X ∗ FE examples. Here, we used NW = 100 and FE = 0.5.

22To evaluate the data quality of each window, we estimated model performance (AUBPRC) on the core set by using 10-fold cross-validation via XGBoost (with default algorithm level hyperparameters) where the training examples in each fold were supplemented by add-on examples in that moving window. To assess whether a property was informative, we calculated correlation between moving window performance and window number for each candidate informative property and also for each of ten random orderings of add-on examples. From this, we calculated a Z score for the observed correlation relative to the distribution of randomly ordered correlations and considered a property to be informative if it exhibited a Z score above 1 or less than −1. The negative and positive sign of a Z score indicates whether there was a negative or positive correlation between performance of moving windows and the property. Although moving window analysis combined one or two add-on sets at a time with the core set, we note that add-on sets may provide a contribution that is greater or lower than that observed here when they are used in conjunction with many other training sets during hyperparameter tuning. Therefore, all training sets and their associated informative properties (Table S3) were included for the weighting process that was optimized during hyperparameter tuning, even for Z scores exhibiting modest departures from zero.

23Two types of hyperparameters were employed for VARITY models. First, there were algorithm-level hyperparameters such as number of trees, tree depth, and learning rate (see Table S4 for a full list). Second, there were hyperparameters in the logistic functions that were used for assigning weights to different training sets (see Table S3 for full list). There are three parameters in each logistic function, namely the maximum weight L, the midpoint x0, and the growth rate k. Each logistic function takes one quality informative property as input. For the training set that has more than one associated quality-informative property, the weights were assigned as the product of multiple logistic functions.

24All hyperparameters were optimized simultaneously by Bayesian optimization with HyperOpt,, (Python package V0.2.2). This process involves random initialization of hyperparameter values (with uniform prior) followed by iteratively repeated trials (here we used 300 trials) with alternative hyperparameter values suggested on the basis of an expected improvement measure calculated via a probability distribution estimated from the performance of previously tried hyperparameter settings. For each trial, the hyperparameter set was evaluated via model performance on the core training set estimated in 10-fold cross-validation where the core training examples in each fold were supplemented with add-on training sets examples and all training examples were weighted with hyperparameter values suggested on the current trial. The metric used here for model performance (the objective function for hyperparameter optimization) was the mean AUBPRC on validation sets (average over ten validation sets).

25To limit model overfitting that can arise from repeated interrogation of the evaluation set with different hyperparameter settings, we used the following procedure to automatically determine the final hyperparameter values: (1) re-order all trials from HyperOpt,, by mean AUBPRC on training sets (averaged over ten training sets) from low to high, (2) calculate a moving window (we used window size 30 for total 300 trials) average of mean AUBPRC on validation sets, (3) define an “early stopping” point at the first moving window (the “fittest” region) for which mean AUBPRC on validation sets begins to descend, and (4) select as final the hyperparameters from the trial within this “fittest” region that achieved the highest mean AUBPRC on validation sets.

26To obtain an independent test dataset based on de novo variation in neurodevelopmental disorders, we downloaded all de novo variants from denovo-db 1.61 (download date 2020-10-20) that corresponded to neurodevelopmental disorder case/control studies. We extracted the subset of rare missense variants (MAF < 0.5%) from 253 candidate neurodevelopmental disease genes for which variation was enriched in neurodevelopmental disease patients. For the purpose of performance evaluation, the variants that appeared in “control” and “case” populations were labeled as negative (putatively benign) and positive (putatively pathogenic), respectively. Because a subset of variants annotated in HGMD had been used to train computational methods against which we wished to compare our results, variants annotated in HGMD (version 2020) were removed, as were variants used in VARITY_R and VARITY_ER model training. There are in total 367 de novo missense variants (321 putatively pathogenic and 46 putatively benign) left for comparison. Considering only variants scored by each of the 23 computational methods to be compared, the final validation dataset contained 215 de novo missense variants (188 putatively pathogenic and 27 putatively benign). For the validation analysis with published variant effect maps, we assembled map scores for a subset of six proteins associated with human genetic disease: CALM1 (MIM: 114180), TPK1 (MIM: 606370), CBS (MIM: 613381), PTEN (MIM: 601728), BRCA1 (MIM: 113705), and VKORC1 (MIM: 608547).

27To evaluate the generalization performance of VARITY models more broadly, we adopted 10-fold nested cross-validation strategy. For each outer loop, the hyperparameters were determined via Bayesian optimization,, based on inner-loop 10-fold cross-validation. The optimized hyperparameters were then used for training based on all data in the outer-loop training set. Performance on VARITY models (and other computational approaches) were assessed on the ten held-out outer-loop validation sets (Figure S2).

28To compare the performance of VARITY models to previously developed pathogenicity predictors, we collected the scores from 25 predictors. For 22 of these, scores were collected from dbNSFP (V4.0b2), while scores for two predictors (MPC and EVMutation) were assembled from links provided by the corresponding papers and scores for DeepSequence were collected from an independent benchmarking study of 11 proteins: UBE2I (MIM: 601661), TPMT (MIM: 187680), TPK1 (MIM: 606370), TP53 (MIM: 191170), SUMO1 (MIM: 601912), PTEN (MIM: 601728), MAPK1 (MIM: 176948), HRAS (MIM: 190020), CALM1 (MIM: 114180), BRCA1 (MIM: 113705), and ADRB2 (MIM: 109690).

29For each variant, the VARITY model provides a log(odds) (“lod”) score (before subsequent transformation to the final score via a sigmoid transformation), which is more positive for variants that are inferred to be pathogenic and more negative for variants that are inferred to be benign. To provide intuition about which features of a given variant drove its score, we defined an “output contribution” score by using Shapley additive explanation values, retrieved from the XGBoost output (with option “predict_pred_interactions” turned on).

30It can be useful to estimate the contribution of individual feature groups or pairs of feature groups to a model’s performance (as opposed to a model’s output values). We therefore defined a “performance contribution” score for all features for each variant used in training. For a positive (putatively pathogenic) variant used in training, performance contribution was the output contribution multiplied by the weight given to the variant during hyperparameter tuning. For a negative (putatively benign) variant used in training, each feature performance contribution is obtained by negating the output contribution and again multiplying by the weight used for that variant. Thus, output contributions moving in the right direction are counted as beneficial to model performance.

31We first assembled a large resource of input features that are potentially informative about variant function for all possible missense variants in ∼18,000 human proteins. To limit the circularity that may inflate the model performance and to limit the effects of training set bias, we excluded features informed by variant annotation (e.g., scores predicted by supervised models) or protein identity (Table S1). For model training and performance evaluation, we next identified from ∼4,000 human disease proteins a “core” set of variants with high quality annotation from ClinVar. This set was limited to rare variants (MAF < 0.5%) for the VARITY_R model. To assess sensitivity of our results to this threshold, we also defined a core set of extremely rare variants (MAF < 10−6) to train and test a VARITY_ER model. For the core set, we labeled variants as “positive” (putatively pathogenic) or “negative” (putatively benign) on the basis of pathogenicity interpretations. We then collected various “add-on” variants, employing potentially less-reliable information about whether non-core variants should be labeled as positive or negative (see material and methods). Add-on variants were drawn from gnomAD, HGMD, HumsaVar, ClinVar (excluding the core set variants), and a handful of multiplexed assays of variant effect (MAVE) studies,,,, (Table S2). Add-on variants were grouped into add-on training sets based on source (e.g., gnomAD), label (e.g., “negative”), and rarity (e.g., “MAF > 0.5%”) so that each add-on set can be employed separately with distinct effects on training VARITY models (see material and methods; Table S3).

32Although VARITY employs a “workhorse” machine-learning method (the gradient boosted tree algorithm XGBoost) (see material and methods), it is distinct from previous machine-learning strategies applied to variant impact prediction. Distinct a priori weights were assigned to different training examples, and these weights were considered in the loss function used for tree optimization. For each training set (core or add-on set), one or more quality-informative properties were used as inputs to a logistic function that assigned weights to each variant. The parameters of each logistic function were optimized as hyperparameters (see material and methods; Figure 1).

33To identify properties related to the utility of different add-on sets for inferring pathogenicity, we ordered the examples in an add-on set (or logically related pairs of add-on sets) by candidate property and carried out a moving window analysis (see material and methods), To evaluate the predictive utility of each window, we estimated the model performance by using 10-fold cross-validation on the core training set where the training examples in each fold were supplemented by all of the add-on examples in that moving window. When a candidate property correlated with moving window predictive utility in an add-on set, this property was identified as an informative property for the add-on set. For example, a moving window analysis of putatively benign variants from gnomAD showed variants with both lower and higher allele frequency to be less useful in training VARITY_R (Figure 2). This agreed with our expectation that variants with lower allele frequency are more likely to be “contaminated” with damaging variants and that common variants might not be representative of the rare variants used for testing. Allele frequency was similarly identified as a quality-informative property for negative-labeled ClinVar44 and HumsaVar45 variants. For several positive-labeled (putatively pathogenic) add-on sets, variants with higher allele frequency were found to be less useful. We also identified both label confidence and mutational accessibility properties (see material and methods for definitions) as quality-informative properties for variants with MAVE scores (Figure S3; Table S3). Although the core training set was considered high quality, we considered ClinVar review stars as a quality-informative property to further weight these training variants.

34The parameters of the variant-weighting logistic function used for each training set were treated as hyperparameters subject to optimization, as were algorithm-level hyperparameters such as learning rate, number of trees, and the maximum depth of a single tree (Table S3 and Table S4). All hyperparameters were tuned simultaneously with Bayesian optimization,,, which consists of many consecutive trials, each making use of a different set of hyperparameter values suggested by the previous trials. For each trial, the hyperparameter set was evaluated via VARITY model performance on the core set examples in 10-fold cross-validation where the core training examples in each fold were supplemented with add-on training sets examples and all training examples were weighted with hyperparameter values suggested on the current trial. To avoid overfitting to the validation set during hyperparameter tuning, we selected a hyperparameter set by using a scheme akin to “early stopping” (see material and methods; Figure S4). With thus-optimized hyperparameters, all 34,087 core set variants were given high weight (≥70%). While nearly all (123,591 of 123,621) add-on set variants were given non-zero weights for VARITY_R, only 1.8% of add-on variants received ≥70% weight (see Figure 3 for the optimized weights of gnomAD47 add-on sets and Figure S5 for all training sets). The effective number of training examples (the sum of weights of all core and add-on set examples) for VARITY_R was 41,898 (24,842 positive examples and 17,056 negative examples) as compared with the 157,708 total of core and add-on examples. A subsequent reweighting step equalized the effective number of positive and negative examples.

35Although the threshold that we used to define rare variants (MAF < 0.5%) in the VARITY_R model is commonly used, we considered the possibility that moderately rare variants (10−6 < MAF < 0.5%) may not be representative of extremely rare variants (MAF < 10−6, for which no allele count was detected in gnomAD). We therefore developed an alternative VARITY_ER model by using only extremely rare variants from ClinVar as the core set (but including moderately rare ClinVar variants as add-on variants). Like VARITY_R, VARITY_ER modeling proceeded via moving window analysis and hyperparameter tuning, and all 18,792 core set variants and 133,848 out of 138,916 add-on set variants received non-zero weights. The effective number of training examples was 24,074 for VARITY_ER, well below the total of 157,708 training examples considered (Figure S6, Figure S7, and Figure S8).

36The final VARITY_R and VARITY_ER models were each trained with both weighted core and add-on sets with optimized hyperparameters and were each used for inference of functional impact for all possible missense variants for each of ∼18,000 human proteins. We observed that cross-validation performance for VARITY_R and VARITY_ER models was only slightly higher for the training set than for the validation sets using the optimized hyperparameter settings (See Figure S4 and Figure S7), suggesting that overfitting, if any, was modest. However, to avoid inflating performance estimates as a result of overfitting, we applied a leave-one-variant-out (LOO) strategy such that each variant used (i.e., having non-zero weight) in training was excluded in turn from the VARITY_R (or VARITY_ER) model used to score that variant, yielding a set of VARITY_R_LOO (or VARITY_ER_LOO) scores.

37To better understand the rationale for each VARITY score, we used Shapley additive explanation values, to assess the contributions made by features and their pairwise combinations. For each variant, contributions were estimated both to VARITY score and, if the variant was used in training, to model performance (see material and methods). For each feature group (Table S1), we estimated the total contribution (including independent and pairwise contribution) to model performance as weighted average of the contribution from all training examples by using weight of each training example that was optimized during hyperparameter tuning (Figure 4). Contributions to model performance at the level of individual features and feature pairs were also evaluated (Figure S9).

38For the VARITY_R model, we found that conservation-based features were the most important for model performance, which is perhaps unsurprising given that they are foundational to most computational predictors of variant function. Among these, Provean, SIFT,, EVMutation, and LRT provided the greatest contribution. The next most important feature group was “delta AA properties,” a set of features capturing differences in physicochemical properties between the missense and wild-type amino acid. Among these, the “delta cyclic property” (capturing whether a proline residue has been substituted for a non-proline residue or vice versa) was the most important. The “delta sulfur property” (capturing whether a cysteine or methionine residue has been substituted for a non-sulfur-containing residue or vice versa) and the “delta hydropathy index” were also very helpful. Next were features related to “surface-accessible area,” representing the extent to which the amino acid is exposed to solvent, and “PPI,” a group of features modeling the extent to which the environment of an amino acid changes upon protein-protein interaction (see material and methods).

39Features and feature groups will of course combine to contribute to VARITY models in a non-additive way. For example, the “IN/OUT Pfam domain” feature (capturing whether an amino acid position falls within a conserved protein domain), which is derived from multiple sequence alignments, contributed to the performance of the “conservation scores” feature group. It is therefore understandable that this feature group masks the predictive value of “IN/OUT Pfam domain” (as indicated by negative pairwise performance contribution values for VARITY_R in Figure 4). A negative pairwise performance contribution for individual features was also observed among nine conservation scores. Because these conservation scores are similar, partial redundancy with one another was to be expected (see Figure S9). We also have observed that the combination of “conservation scores” with other feature groups such as “delta AA properties” and “PPI” provides a positive (synergistically favorable) pairwise feature contribution to model performance. Features contributing to each VARITY_ER model prediction score were also identified, and the performance contributions of individual and grouped features was again analyzed with the weighted training examples, yielding similar results (Figure S10 and Figure S11).

40To assess VARITY, we examined missense variants from neurodevelopmental case/control studies, focusing on 253 genes enriched for de novo variation in neurodevelopmental disease patients., For performance comparison, we considered all rare de novo variants in these genes that had not been used in training by VARITY and were not annotated by HGMD (see material and methods). After retaining only variants that had MAF < 0.5% and that had been scored by all 23 of the methods to be compared, there were 188 variants in the case group and 27 in the control group.

41We wished to judge performance using both receiver operating characteristics (ROC) and precision versus recall analysis. Precision versus recall analysis has the advantage of being more intuitively related to the needs of most users. However, precision varies as a function of the prior probability, and this prior is by default determined by the fraction of test examples with a Positive (putatively pathogenic) label. Moreover, this default prior may not have any relationship with the prior probability that is most relevant to a clinical geneticist using our output. To enable greater intuition and fairer comparisons, we therefore transformed each empirical precision-recall curve to the curve corresponding to a balanced (50% probability of pathogenicity) prior. Thus, our performance measures were as follows: area under the ROC curve (AUROC), area under the balanced precision recall curve (AUBPRC), and recall at 90% balanced precision (R90BP; see material and methods for details). For this de novo variant application, we used the case/control status of the patient(s) in which each variant appears as a proxy for pathogenicity.

42VARITY_ER numerically outperformed all other methods in terms of AUBPRC. When performance was measured with the subset of variants that had been scored by all methods, AUBPRC performance improvement achieved statistical significance (p < 0.05) relative to every method except MPC (p = 0.125) and M-CAP (p = 0.099). Moreover, when we used larger variant sets that had been scored by MPC (285 variants in case group and 38 variants in control group) or M-CAP (320 variants in case group and 45 variants in control group), VARITY_ER’s performance improvement was statistically significant (p = 0.048 and 0.029 for MPC and M-CAP, respectively). When the score threshold for every method was tuned to achieve 90% balanced precision, VARITY_ER recovered 10% more (presumed pathogenic) variants from the case group than the closest non-VARITY method. VARITY_ER outperformed VARITY_R with a 2% greater AUBPRC (p = 0.122), which may be explained by the fact that most de novo missense variants are extremely rare. However, VARITY_R also numerically outperformed all other methods (see Figure 5 and Table S5 for AUBPRC comparison; Figure S12 and Table S6 for AUROC comparison; Table S7 for individual comparisons with MPC and M-CAP).

43As another form of independent validation, we wished to assess agreement between impact scores and systematic experimental assessments of variant function. We therefore evaluated the ability of scores from each computational method to predict the quantitative scores from each of six recently published missense variant effect maps on human disease related genes.,,,,, Because five out of six maps (all except VKORC1 map) had been indirectly used in VARITY training, we included leave-one-variant-out scores (VARITY_R_LOO and VARITY_ER_LOO) for comparison. VARITY_ER_LOO showed best correspondence in terms of Pearson correlation coefficient (PCC) with variant effect map scores (PCC = 0.405), significantly outperforming all other non-VARITY predictors (p < 0.05; Figure 6 and Table S8 for PCC results; Table S9 for Spearman’s rank correlation coefficient [SRC] results).

44To further compare generalization performance of VARITY models with previously developed computational predictors of missense variant effects, we used the high-quality core sets from ClinVar. To avoid possible circularity from model training and hyperparameter tuning, we adopted a nested cross-validation strategy (see material and methods).

45VARITY_R significantly outperformed all 23 methods on rare variants from ClinVar (MAF < 0.5%; VARITY_R core sets) by all criteria (i.e., AUBPRC, R90BP, and AUROC; p < 0.05 via one-sided paired t test), recovering 13% more pathogenic rare missense variants at 90% balanced precision than the next-best method, REVEL (Figure 7 and Table S10 for AUBPRC comparison; Figure S13 and Table S11 for ROC comparison). EVMutation and DeepSequence had few scores available for core set variants and, because we limited the overall comparison to variants assessed by all methods, their inclusion would have dramatically reduced the space of test variants. However, direct comparison of VARITY_R with each method with the intersection of test variants scored by each respective method showed that VARITY_R significantly outperformed both methods (Table S12). Methods designed to evaluate amino acid substitutions generally outperformed methods that were designed for nucleotide-level changes (see predictors marked with “(•)” in Figure 7). Indeed, among nucleotide-level methods, only CADD (which may have benefited from amino-acid-level PolyPhen2, scores as training features) was able to outperform an amino-acid-level method.

46To evaluate the extent to which VARITY’s success depended on the variant-weighting scheme, we retrained a “VARITY_R_unweighted” model, giving equal weight to all core and add-on variants. This reduced AUBPRC performance from 92% to 90.9% (p = 6.5 × 10−5; Student’s t test comparing the two sets of outer-loop validation performance values). We expected that performance for the no-variant-weighting version of VARITY would fall further as more add-on training examples with uncertain quality are included. To evaluate this, we included additional gnomAD variants that were originally excluded from the putatively-benign add-on training set because they had never been seen in a homozygous individual. After including these variants in training, the nested cross-validation AUBPRC performance of VARITY_R_unweighted was reduced further from 90.9% to 89.9% (p = 2.1 × 10−4).

47We also evaluated VARITY_ER performance by using the core set of extremely rare variants from ClinVar (MAF < 10−6; VARITY_ER core sets), again with 10-fold nested cross-validation. VARITY_ER outperformed all 23 methods by all criteria (p < 0.05 via one-sided paired t test), recovering 12% more pathogenic rare missense variants at 90% balanced precision than the next-best method REVEL (Figure S14 and Table S13 for AUBPRC comparison; Figure S15 and Table S14 for ROC comparison). Although EVMutation and DeepSequence were again excluded, VARITY_ER outperformed EVMutation in a pairwise comparison based on extremely rare core set variants with available EVMutation scores (Table S15). Only 16 extremely rare benign core set variants had DeepSequence scores, leaving too few variants in each held-out subset during 10-fold nested cross-validation to accurately evaluate performance.

48Evidence to evaluate variant pathogenicity is typically less available for rare as opposed to common variants. Detecting population-level correlation between a trait and a genotype has reduced power at lower minor allele frequencies and is futile for extremely rare variants. Experimental functional evidence is also less available for lower-allele-frequency variants. For these reasons, and because the majority of unique clinically observed variants are rare, there is a particularly strong motivation to improve computational methods for inferring the pathogenicity of rare variants.

49Here, we developed predictors of variant pathogenicity that were specifically optimized for performance on rare and extremely rare missense variants. Given previous evidence that common variants are not representative of rare variants in the context of training pathogenicity models, perhaps the primary advantage of VARITY models is that rare variants have been given greater weight in training. Other key advantages are that VARITY makes use of large sets of potentially biased or error-prone training data while optimizing the potential for these add-on variant sets to increase performance. VARITY models benefitted not only from the process of filtering and weighting add-on examples but also from the use of features that have, either individually or in combination, not been used previously. To limit protein-specific bias and circularity, VARITY models also excluded features that are informative about protein identity or that may have been informed by previous pathogenicity annotations.

50VARITY models performed well when evaluated via an independent set of de novo missense variants collected from neurodevelopmental case-control studies. Here, we only considered the subset of de novo missense variants that were in genes previously reported as being enriched for de novo missense and/or likely-gene-disruptive variants in neurodevelopmental disease cases. Because this certainly excluded some disease-relevant genes/variants, it would be interesting in future studies to expand the gene/variant set based on more permissive identification of relevant genes/variants. VARITY models might also be used for identification of candidate causal variants in other genes, e.g., in patients for whom a predicted-damaging variant has not been observed in an already-known disease gene.

51More broadly, we used nested cross-validation to estimate the generalization performance of VARITY models on the core set of ClinVar variants. We note that the set of ClinVar variants is influenced by historical ascertainment biases, in that they are enriched in disease genes discovered earlier, and some annotations were most likely influenced by early computational methods, such as BLOSUM and SIFT., Therefore, although ClinVar variants enable useful measures of relative performance, it is difficult to estimate performance on an absolute scale for variants that will be observed in the future. Although the de novo variant evaluation set was free from these biases, it was necessarily limited in scope so that VARITY’s generalization performance should in the future also be evaluated with additional independent test sets.

52We provided information about which features and feature combinations were most useful in training VARITY models and feature contributions for the prediction made for each variant. We also note that the quantitative VARITY scores themselves can provide some intuition, e.g., there is a clear trend relating VARITY_R scores to the probability of pathogenicity (Figure S16).

53The gradient boosted tree method has the advantage that training and testing examples for which some feature values are missing can be used directly without requiring imputation. Nevertheless, predictions may be less accurate for variants that are missing important features such as conservation-based scores and structure associated features. For example, if we remove all structural features from the core set, VARITY_R nested cross-validation performance (via AUBPRC) dropped from 92% to 91.5%.

54For VARITY training, we eliminated protein-identifying features to reduce protein-specific biases in our training set. However, we might want to augment VARITY scores in the future by explicitly deriving a prior probability of missense variant pathogenicity (i.e., the probability of pathogenicity before considering the identity of the missense variant) for each protein or protein region. However, because these priors may vary with the application setting, we suggest that the user should be explicitly informed about these priors and given the option to adjust them.

55There are several avenues for future improvement of VARITY models. More training examples could be added, e.g., adding putatively benign variants from primates as PrimateAI has done. One might also develop specialized models for different classes of proteins or variants, e.g., those in predicted membrane-spanning domains. More quality-informative properties could be added to improve add-on set weighting. More features could also be added, e.g., based on knowledge of genetic or protein-protein interactions. We are also currently lacking features related to the typical mode of inheritance for a given gene and associated disease (e.g., dominant or recessive), the typical mechanistic class (e.g., gain or loss of function), or typical penetrance or variable expressivity of variants. Although adding mode-of-inheritance features would be complicated by the fact that this information is often unavailable and can differ for different variants within the same gene, it could improve performance in general.

56To our knowledge, no previous computational variant effect prediction method has used the strategy of weighting variants used in training. However, it is not uncommon in machine learning to weight training examples. For example, the AdaBoost ensemble algorithm assigns higher weight to examples misclassified by previous classifiers in each learning iteration. Where multiple annotators have provided possibly noisy labels, the weight and the true label of each example can be learned together with the model via a maximum-likelihood approach. We also note that the modeling approach we describe is just one way to address the problem of model learning with noisy labels, which is the subject of intensive research especially within the deep learning field.,,,, Given that variant effect interpretation is greatly limited by the availability of high-confidence unbiased training data, this field will be worth watching closely for future inspiration.

57Other models tuned for other types of variation could employ this framework by simply changing the core and add-on training examples. For example, a computational model specialized for predicting the pathogenicity of variation in membrane-spanning domains could be trained by limiting the core set to variants from regions predicted to be in the membrane while still using all other variants as add-on training sets subject to differential weighting. Models for predicting deleteriousness of nucleotide-level variation in non-coding regions or spliced regions would require more extensive separate assembly of relevant features and both core and add-on variant sets but could then benefit from the algorithmic framework we describe here.

58As the performance of computational variant effect predictors improves, an important unsolved problem is how to translate these gains to have a commensurate impact on clinical variant interpretation. The increasing quantity of systematically collected large-scale functional evidence about missense variation has led to new recommendations for its use in the context of clinical interpretation. We expect that a similarly quantitative evaluation of computational predictions will support a more nuanced use of computational evidence in clinical variant interpretation.

59Users may search, browse, and download both VARITY_R and VARITY_ER scores (with associated feature contributions) for any of ∼18,000 human proteins via the VARITY web portal (see web resources).