A total of 193,714 missense variants (164,574 benign, 29,140 disease-causing) were collected from ClinVar (Landrum et al., 2018), gnomAD (Karczewski et al., 2020) and Dutch genome diagnostic laboratories (VKGL, 2019), which could be linked to protein structures, either directly or through homology with a sequence identity cut-off of 40%. This cutoff was selected based on previous research that suggests that a 40% identity corresponds to a good likelihood of functional equivalence (Pearson, 2013). Missense variants were mapped onto protein structures using 3DM systems as a guide (Kuipers et al., 2010). Independent test datasets were obtained from studies based on BRCA1 (Findlay et al., 2018), Gunning et al. (2021) and the InSIGHT database (Thompson et al., 2014). This resulted in a total of 217,679 missense variants that could reliably be mapped onto 57,551 structures; 25,856 structures were mapped to 40,369 pathogenic variants, and 31,695 structures were mapped to 177,310 benign variants. It should be noted that, at this stage, the structures are mapped regardless of the experimental method used for their determination. Missense variants from ClinVar were incorporated if they had a review status of at least one star, excluding those with conflicting interpretations. “Benign” and “Likely benign” ClinVar variants were included and categorized as benign, while “Pathogenic” and “Likely pathogenic” variants were incorporated and classified as pathogenic. The gnomAD variants with a minor allele frequency higher than 0.1% were selected and labeled as benign.

Our in-house database, HSSP (Touw et al., 2015) was consulted to obtain structure-based sequence alignments. Position-specific scoring matrices (PSSMs) were constructed for the alignments using PSI-BLAST (Altschul et al., 1997) with single iteration. Each of the PSSMs were then mapped back onto their respective structures using the PSSMGen package (https://github.com/DeepRank/PSSMGen).

The network used in our study includes a sequential organization of three 3D convolutional layers, alternating with one 3D max pooling layer followed by two fully connected layers (Figure 1). We include batch normalization layers, in addition to dropout layers between the fully connected layers to regularize the model. Details of the architecture are provided in Table 2 and the complete schema is provided in Supplementary Figure S2. Each 3D convolution layer comprises a set of learnable filters that traverse the input space (depth, height and width) with a stride of 1, capturing local spatial patterns in the variant environment. The output from convolution operations, i.e., the computed feature maps are transformed by a rectified linear activation function (ReLU), which allows the network to identify and extract meaningful spatial features. This is followed by dimension reduction using max pooling operation and a final 3D convolutional layer with ReLU. The transformed output is then flattened to a one-dimensional vector that serves as an input to two fully connected layers. The two final layers integrate the features and apply a set of weights that are optimized during the training step to map extracted features to target classes. The output is then passed through the softmax function which provides the final classification score, a probability estimate between 0 and 1, each for benign and pathogenic classes.

We performed 10-fold cross validation experiments while ensuring that the missense variants in the training and validation sets are from different proteins, to avoid type 1 circularity in predictions (Heijl et al., 2020). The test dataset included missense variants independent from the 10-fold training and validation sets. Most genetic variation is neutral, and it is therefore rather common to observe a higher number of benign variants than pathogenic variants in the training data, which has the potential to bias training and performance. We thus constructed balanced subsets of randomly sampled benign and pathogenic missense variants for each of the 10-fold runs. For efficient memory handling, we employed training in mini-batches of 256 variant instances which amounted to ∼1,200 mini-batches per epoch. An epoch refers to a single pass through the complete training data during which the model weights are adjusted to minimize the error between predicted and true label for each input. With the input dataset, one epoch in our approach referred to one pass through more than 280,000 variant instances. We used the AdamW optimizer (Loshchilov and Hutter, 2019) with a learning rate of 0.001 and weight decay of 0.005 to train our model for 10 epochs. We used cross entropy loss during training, which attempts to minimize the differences in probability distributions between predicted and ground truth labels by adjusting weights. A dropout rate of 0.5 was used to regularize the model. The hyperparameters including number of convolutional layers, number of max pooling layers, grid size, were optimized based on performance on validation set across 10 folds, starting from default parameters of the parent DeepRank.

Two metrics, Matthews Correlation Coefficient (MCC) and accuracy, were used to evaluate the performance of DeepRank-Mut. The primary metric used was MCC, as it offers a reliable statistical measure by taking all four categories-true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN) into account, proportional to the size of the binary classes (Eq. 2). The usefulness of MCC over accuracy or F1 scores for binary classification has been demonstrated previously (Chicco and Jurman, 2020). M C C = T P × T N − F P × F N T P + F P × T P + F N × T N + F P × T N + F N (2)

For comparative evaluation with popular state-of-the-art variant effect predictors, we used precomputed pathogenicity prediction scores of 8 algorithms from dbNSFP v4.3 database (Liu et al., 2020; 2011), including SIFT4G (Vaser et al., 2016), PolyPhen2 (Adzhubei et al., 2010), MutationTaster (Schwarz et al., 2014), MutationAssessor (Reva et al., 2011), FATHMM (Shihab et al., 2013), VEST4 (Carter et al., 2013), PROVEAN (Choi and Chan, 2015) and MutPred (Li et al., 2009), as well as prediction scores from Helix (Vroling and Heijl, 2021). Where available, we used “converted rankscores” from dbNSFP to ensure that a higher score always indicated higher likelihood of pathogenicity. We excluded meta predictors from this comparison, as well as those that combine annotations from other tools, to account for methods that rely on first principles to predict functional effects of missense variants.

Training, validation, and test sets are often generated using a simple random split. However, this can result in over fitting and misleading results due to data leakage between the training and evaluation sets (Heijl et al., 2020). Data splitting at the level of proteins or genes, where training sets never include any data samples from proteins that occur in the validation or test set is used to mitigate this. We split our dataset into 10 pairs of training and test sets, each containing 90% and 10% of the full dataset, respectively, allowing 10-fold cross validation on the full dataset. Independent test sets were gathered from three studies, as described in methods, to aid in the final assessment of the tool’s performance. These test sets have been selected as they not only cover genes in-depth (Thompson et al., 2014; Findlay et al., 2018) but are also aimed at benchmarking pathogenicity predictors specifically (Gunning et al., 2021).

After splitting the data, the balanced subsets of randomly sampled benign and pathogenic variants, each mapped to at most three structures, comprised a total of ∼50,000 instances in the training set, ∼4,700 instances in validation and 6,571 in the test set, per fold. The test set was kept identical across all cross-validation folds for an unbiased evaluation of the model.

DeepRank-Mut retains its modularity in implementing data pre-processing steps and training the deep neural network, similar to its parent DeepRank (Renaud et al., 2021). It allows for flexibility in tasks including feature calculations, setting the grid size and grid resolution, data augmentation, as well as optimizing hyperparameters of the neural network. The base requirements of DeepRank-Mut include a dataset of variants with labels (benign or pathogenic), a dataset of variant-structure maps where each variant is linked to a 3D structure (either experimentally determined or evolutionarily related), a dataset of 3D structures and an optional dataset of PSSM profiles derived for each structure. As detailed in the methods, the framework computes physicochemical properties of the amino acid at the variant site as well as its environment within a radius of 10Å, followed by voxelization to encode the atomic neighborhood of residues (Figures 1A, B). Our approach relies on leveraging local properties of sites characteristic of benign or pathogenic variants, as pathogenic variants generally tend to occur in regions important for structural/functional integrity of the protein (Iqbal et al., 2020), like its hydrophobic core. We thus compute a total of 31 features (Table 1), encompassing structural and sequence-based properties, for the residue at the variant site and residues spatially proximal to it. The computed features are mapped to a 3D grid where each voxel is parameterized with the feature channels (Figure 1C), which is then followed by data augmentation. As a given variant environment can differ in orientation within or across proteins, the data augmentation step accounts for rotational invariance, thereby improving the model’s robustness to variations in input data (Supplementary Figure S3). From our dataset of structures and missense variants, we generated ∼300,000 augmented grids per fold dataset, which were used as input to 3D-CNN (Figure 1D). Each augmented 3D grid is treated as a separate variant instance, thus our model outputs 6 predictions per missense variant (origin grid +5 augmented grids) which are averaged to give one final classification score.

Our approach achieved a mean accuracy of 0.77 and an average MCC score of 0.52 across the test datasets, with an average sensitivity (true positive rate) of 0.75 and an average specificity (true negative rate) of 0.78 (Figure 2; Table 3).

We examined our model’s predictive ability by analyzing missense variants in the test-set that were consistently predicted incorrectly across all 10 folds. We explored the aspects that promoted incorrect classification. A total of 2,883 missense variants were found to be incorrectly classified across the cross-validation experiments, of which more than half (1,732) consisted of misclassified pathogenic variants. We computed relative solvent accessibilities (RSA) for each variant residue, by dividing their absolute solvent accessibilities in Ų by their maximum allowed solvent accessibilities obtained from Rost and Sander (Rost and Sander, 1994). Residues were categorized as solvent-exposed if the RSA values were >20% and buried if below 20%. Using these a substantial proportion of the misclassified pathogenic variants was found to be solvent-exposed (Supplementary Table S2).

We constructed 2 × 2 contingency tables based on the correct and incorrect classifications with respect to solvent accessibility of the associated variants. Figure 5 illustrates the role of solvent accessibility in the predicted outcomes. The misclassified variants pertained to solvent-exposed pathogenic variants and buried benign variants (Figure 5A, odds ratio = 0.27). That we are relatively successful in predicting pathogenicity in buried variants is consistent with the notion of buried enrichment of pathogenic variants (Iqbal et al., 2020; Savojardo et al., 2020). The distribution of raw atom-level solvent accessibility values across benign and pathogenic classes calculated in our approach is illustrated in Supplementary Figure S9. Two reasons for the quality of the predictions could be postulated: a) considering the contribution of SASA in the model’s performance, it is likely that the model is unable to generalize on missense variants that fall outside the purview of typical SASA distribution observed in benign and pathogenic variants, or b) the 3D input grids for solvent-exposed missense variants are sparsely populated which leads to a lack of discernible patterns/features for the model to learn from.

We created separate training subsets of buried and solvent-exposed variants to understand 3D-CNN’s generalizability to either subset. We observed that the predictions on pathogenic variants improved with the model trained on buried missense variants alone, however, this model misclassified much of the benign variants, whereas the model trained on solvent-exposed variants alone showed a performance comparable to that of the full model trained on all variants (Figure 5B; Supplementary Figure S5). It is possible that the presence of a large proportion of solvent-exposed variants in our training data may have impacted the performance (Supplementary Figure S9). Furthermore, to assess whether sparsity of 3D grids of solvent-exposed variants affected the model’s performance, we calculated the ratio of solvent (void) voxels to atom-contained (non-void) voxels in the 3D grids in test dataset and compared the distribution of these ratios against the corresponding pathogenic and benign prediction scores. We find no correlation for pathogenic variants (Pearson’s r = −0.09), while we find that the presence of void voxels is weakly indicative of correct classifications for benign variants (Pearson’s r = −0.24) (Supplementary Figure S10). This overall suggests that grid sparsity has weak effect on the correct classification of benign variants, whereas the incorrect classifications of solvent-exposed pathogenic variants is possibly due to other reasons, such as lack of function-specific features, and/or incomplete knowledge of their interaction partners.

Since data augmentation and feature normalization strategies, typically used to circumvent lack of generalizability and potential biases, are already incorporated in our approach we experimented with inclusion of other structural features: secondary structural content and normalized B-factors. The premise behind use of secondary structural content was based on the report by Abrusán and Marsh (Abrusán and Marsh, 2016), who showed differences in the ability of alpha helices and beta strands to tolerate mutations. Secondary structural assignments for protein structures were obtained from our in-house database (DSSP v.3.1.4) (Kabsch and Sander, 1983), and were stored as one-hot encoded features in 3D grids. B-factors or temperature factors are obtained from X-ray crystallography experiments that indicate atomic flexibility in the protein’s crystalline state, and are known to correlate with flexible regions of the protein. Based on the earlier reports of active/functional sites associated with lower B-factors as compared to non-functional residues (Sun et al., 2019), we used normalized B-factors as a feature to potentially capture such differences. However, the two additional features did not serve as strong determinants of pathogenicity (Supplementary Figure S6). The relatively low quality of predictions for solvent-exposed pathogenic variants and buried benign variants could be due to lack of function-specific features.

We further investigated DeepRank-Mut’s generalizability with respect to mutation mechanisms. Most available pathogenicity predictors do not make a distinction between different types of mutation mechanisms such as loss-of-function (LoF) or gain-of-function (GoF), that are often linked to mode of inheritance. LoFs are function-disrupting mutations that usually cause damage to protein structures and are straightforward to comprehend and identify, as they are generally not tolerated at sites of high structural and/or functional importance, and lead to degradation of the protein. In contrast, GoFs exhibit milder effects on protein stability while giving rise to altered protein functions that lead to diseases (Gerasimavicius et al., 2022). In terms of mode of inheritance, autosomal recessive (AR) diseases are predominantly linked to LoFs, while autosomal dominant (AD) diseases manifest through mechanisms such as GoFs, dominant-negative mutations (DN) as well as through LoFs, i.e., haploin sufficiency (Veitia et al., 2018).

To understand how DeepRank-Mut generalizes on distinct modes of inheritance of pathogenic variants, we split our test datasets into variants with AD inheritance (n = 1,363; 550 benign, 813 pathogenic) and variants with AR inheritance (n = 563; 244 benign, 319 pathogenic), based on information obtained from ClinVar (Landrum et al., 2018). Only a smaller subset could be mapped to crystal structures: 585 structures mapped to 515 AD variants, and 77 structures mapped to 132 variants. We did not filter the AD dataset further to segregate mutations into haploinsufficient genes (LoFs) and non-LoFs (GoFs, DNs), due to lack of detailed annotations of non-LoFs in ClinVar. However, it is worth noting that mutations in the AD dataset could consist of higher proportion of LoFs than non-LoFs due to smaller mutational target for non-LoFs, i.e., fewer mutations alter protein function than disrupt it. Figure 6 illustrates a marked difference in the model’s performance between the two datasets, suggesting dependence on underlying effects of the variant on the protein. It is apparent from the figure that our model is able to generalize AR mutations (LoFs) better than AD mutations (LoFs and non-LoFs). Details on the pathogenicity predictions obtained for AD and AR datasets, are provided in Supplementary Table S3. To further examine our relative success in correctly classifying buried pathogenic variants and AR variants we analyzed the distribution of solvent accessibility in the AD and AR datasets. Interestingly, the typical distribution of solvent-exposed benign variants and buried pathogenic variants was found to be more pronounced in AR datasets than in the AD datasets (Figure 7), explaining the relative success of our model in distinguishing LoF mutations from benign (Supplementary Figure S7, MCC = 0.67).