For all human proteins in the Swiss-Prot database (UniProtKB/Swiss-Prot, Release 2018, 04) [6], we downloaded and curated their available experimental structures from Protein Data Bank (PDB) [16]. For those proteins without available experimental structures, their homology models from ModBase [25] were adopted. The processing details are described in Supplementary Fig. 1 and Supplementary Methods.

The humsavar.txt (2018, 04 of 25 Apr 2018) file provided by the UniProt FTP server contained manually annotated AAS, and is the main source for preparing the AAS datasets [6]. According to its documentation, AASs labeled with “Disease” serve as positive samples (i.e., daAAS), while those labeled with “Polymorphism” (similar to “neutral” or “benign”) are negative samples (i.e., nAAS). Another two data sources of AASs are the 1000 Genomes Project (1000G, Release 2 May 2013) [26] and the VariSNP (Release 16 Feb 2017) [27]. After removing those whose mutant allele frequency is less than 1%, the remaining AASs from 1000G and VariSNP were regarded as nAASs. After mapping to experimental structures and keeping only one instance for duplicates, the AASs from these three sources constituted the TotalDataset.

The TotalDataset was then split into TrainDataset and TestDataset 1 according to the following procedure: After an all-to-all BLASTP (E-value $\le$0.01) run among the protein sequences in TotalDataset [28], the proteins with sequence identity to any other one of less than 30% were selected, and their AASs were regarded as the TestDataset 1. The AASs in the remaining proteins served as the TrainDataset. On the other hand, the AASs in humsavar.txt located at homology structure models constituted the TestDataset 2. The detailed AAS mapping process is also illustrated in Supplementary Fig. 2.

A series of subsets were further prepared from these three datasets. In detail, the subsets were organized according to the proportion of daAASs on each protein: For proteins containing only daAASs or nAASs, their AASs constitute the “Pure” subsets; other subsets, namely, “Mixed” subsets at different mixing levels, were constructed by selecting proteins within specific ranges of daAASs proportion, including open interval between 0.0 and 1.0 (denoted as ]0.0, 1.0[) and close intervals of [0.1, 0.9], [0.2, 0.8], [0.3, 0.7], and [0.4, 0.6].

Based on protein structures, sequences, and sequence alignments, a total of 212 candidate features were calculated to characterize each AAS from various perspectives (gene/protein, AAS site, substitution itself, etc.). The full list of features is described in Supplementary Methods.

Each dataset was organized as a feature matrix with each row representing an AAS and each column corresponding to a feature. The missing data were filled with the mean values derived from the TrainDataset. Furthermore, each feature column was transformed into Z-score, i.e., all values in the same feature column were standardized with the mean and standard deviation derived from the TrainDataset.

Random forest (RF) is a machine learning framework consisting of an ensemble of decision trees, and has been widely applied in classification problems [24]. In the training stage, each tree is trained with only a part of the training samples to decrease over-fitting. In the stage of prediction, the consensus output from the majority of the trees was taken as the final result. Considering its successful application in many bioinformatics studies [29], we chose the RF technique in this work for automatic feature selection and model training.

We adopted a wrapper strategy in the feature selection procedure. That is, we iteratively generated feature subsets, and evaluated them by using the 10-fold cross-validation performance of RF classifiers trained on them accordingly. The brief logic was that we added two features with the most contribution and removed one feature with the least feature importance in each iteration, and this operation would be terminated to deduce the feature redundancy if the adding of new features would not improve the performance of the classifier any more. This procedure has been demonstrated effective in one previous study [30]. Notably, we conducted the cross-validation at the protein level, where AASs from the same protein were required to reside in the same fold. That is, AASs were grouped at the protein level, namely, Group10Fold cross-validation (G10F-CV). The purpose of this operation was to reduce the level of type 2 circularity, which has been discussed in-depth previously [31]: If protein/gene-level features were used, there would exist over-fitting in predicting the disease-association of those AASs whose proteins contain AASs used in the training procedure.

After obtaining the final feature subset, we determined the optimal RF hyper-parameters by finding the maximum area under the ROC curve (AUC) of G10F-CV in a random search of 1000 trials (Supplementary Table 1). Then, by specifying the hyper-parameters with the optimal values, the final classifier was re-trained on the entire TrainDataset. All these processes were conducted by using the Scikit-learn toolkit v0.19.0 [32].

The two independent testing datasets, i.e., TestDataset 1 and TestDataset 2, were utilized to evaluate performance and to compare with other predictors by adopting standard performance measures, including accuracy (ACC), AUC, Matthews correlation coefficient (MCC), sensitivity (Sen), specificity (Spe), positive predictive value (PPV), and negative predictive value (NPV) (definitions in Supplementary Methods).

In addition to the evaluations on the full datasets of TestDataset 1 and TestDataset 2, more in-depth evaluations were performed on the ‘Pure’ and ‘Mixed’ testing subsets, which are designed to provide an examination of the extent of type 2 circularity [31].

Our predictor was compared with seven popular tools, including SIFT (version 5.2.2) [8], HumDiv-trained PolyPhen-2 (PPH2_HD) (version 2.2.2r405c) and HumVar-trained PolyPhen-2 (PPH2_HV) (version 2.2.2r405c) [11], PROVEAN (version 1.1.5) [13], FATHMM’s weighted method (FATHMM-W) and unweighted method (FATHMM-U) (version 2.3) [14], and PANTHER-PSEP (version 1.01) [15]. The settings of these tools are described in detail in Supplementary Methods.

In summary, 6430 experimental structures of 5278 proteins and 4238 homology models of 3682 proteins were obtained for preparing the AAS datasets (Supplementary Fig. 1).

After initial data cleaning, the humsavar dataset retained 29,328 daAASs and 39,679 nAASs. Among them, 11,157 daAASs and 7072 nAASs were mapped to experimental structures, and 2471 daAASs and 2778 nAASs were mapped to ModBase [25] homology models. Overall, the mapped percentages were 46.5% and 24.8% for daAASs and nAASs, respectively. For the VariSNP dataset, the initial cleaning, keeping those with MAF at no less than 1%, and mapping to UniProt canonical sequences, resulted in 18,233 nAASs. Among them, 1281 (7.0%) were mapped to experimental structures. For the 1000G dataset, 3296 of 34,791 (9.5%) nAASs were mapped to experimental structures. After integration and proper data partition described in section 3.1, these mapped AASs were separated into TrainDataset, TestDataset 1, and TestDataset 2 (Supplementary Fig. 2).

The TrainDataset consists of 14,117 AASs mapped on the experimental structures of 1979 proteins, with 7385 daAASs on 611 proteins and 6732 nAASs on 1750 proteins (Table 1). The ratio of daAAS to nAAS is highly balanced (7385:6732), which will ensure that the trained classifier would not suffer from data bias [33, 34].

The TestDataset 1 contains 5485 AASs mapped on the experimental structures from 834 proteins, with 3772 daAASs on 321 proteins and 1713 nAASs on 733 proteins (Table 1). The TestDataset 2 contains 5249 AASs mapped on homology models of 1418 proteins, with 2471 daAASs on 383 proteins and 2778 nAASs on 1208 proteins (Table 1). The TestDataset 1 and TestDataset 2 were utilized to evaluate the performance of the trained predictor on AASs with features extracted from experimental structures and reliable homology models, respectively. Notably, both the TestDataset 1 and TestDataset 2 have no overlap with TrainDataset at either the substitution level or the protein level, which would presumably avoid the two types of circularity that may cause overly optimistic estimation of performance [31]. Moreover, the sequence identity of any pairwise comparison between TestDataset 1 and TrainDataset is less than 30%, but 947 of 1418 proteins in TestDataset 2 have high scoring pairs with sequence identity $\ge$30% with those in TrainDataset, indicating that TestDataset 1 is more rigorous than TestDataset 2. The details of TrainDataset, TestDataset 1, and TestDataset 2 are provided at http://www.wdspdb.com/AAS3D-RF/.

We inspected the datasets by examining the “Pure” and “Mixed” subsets (Fig. 1A,B). A substantial portion of AASs came from the “Pure” subsets (32.6% in TestDataset 1 and 65.1% in TestDataset 2), i.e., many proteins contributed only daAASs or only nAASs. These “Pure” or “Mixed” testing subsets can provide more detailed comparisons of the predictors’ performance, as some predictors confounded by type 2 circularity could not be evaluated properly on the entire testing dataset [31].

Fig. 1.

The composition of different AAS subsets and the performance of eight predictors on them. Each subset was prepared by incorporating the AASs whose proteins harbor a specific range of daAASs proportion (interval labels under the X axis) for TestDataset 1 (A) and TestDataset 2 (B). Among them, “All” represents the whole testing dataset, while “Pure” contains the AASs whose proteins harbor either daAASs or nAASs only. The numbers of daAASs and nAASs of each subset are given on top of each bar. The percentage of the AAS number in each subset is also indicated. The AUC scores on different AAS subsets are plotted for TestDataset 1 (C) and TestDataset 2 (D).

A similar pattern can be found in the TrainDataset, where a large part of AASs (44.3%) are from the “Pure” subset (Supplementary Fig. 3). In addition, considering that some features adopted in the classifier describe the whole gene or protein but not the AAS itself, the trained classifier may be skewed to classify protein but not the AAS (type 2 circularity). As we had adopted protein level cross-validation (G10F-CV in this work) in the training procedure [31], this skewness could presumably be removed.

A total of 21 features (Supplementary Table 2) were finally selected after the iterative feature selection procedure. After randomly searching 1000 hyper-parameter combinations of RF, we obtained the optimal one (n estimators: 338, max depth: 10, max features: 0.5233959) corresponding to the highest AUC (0.873) in the G10F-CV. By specifying the optimal hyper-parameter combination, we re-trained the final RF classifier, namely AAS3D-RF, on TrainDataset. We also implemented an automatic pipeline that can calculate the required features for a given AAS, and can then load AAS3D-RF to make predictions. The codes were implemented in Python 2.7 and are available at http://www.wdspdb.com/AAS3D-RF/ and https://github.com/PKU-XiongYao/AAS3D-RF.

The measures of the performance on TestDataset 1 and TestDataset 2 are listed in Table 2 and Supplementary Table 3. The ACC and MCC values are 0.811 and 0.591 for TestDataset 1, and 0.839 and 0.684 for TestDataset 2, demonstrating its superior overall performance. The better performance on TestDataset 2 may partially stem from the similar proteins between TestDataset 2 and TrainDataset, since we did not require that the proteins having HSP with sequence identity $\ge$30% be excluded from TestDataset 2. When removing this part of data from TestDataset 2, we can observe that the MCC drops from 0.684 to 0.627, supporting our speculation very well (Supplementary Table 3). In the real-world application, one may often need to predict the disease-association of an AAS whose protein is similar to some protein in the training set, so the expected performance would presumably be better than reported here.

We compared the performance of AAS3D-RF with seven other popular tools. Except FATHMM-W (discussed below), AAS3D-RF outperforms all the other predictors on TestDataset 1 and TestDataset 2 in terms of ACC, AUC, MCC, and PPV (Fig. 2A,B, Supplementary Table 4). In particular, the MCC values are 8.7 and 11.3 percentage points higher than the second-best predictors (except FATHMM-W) on TestDataset 1 and TestDataset 2, respectively (Supplementary Table 4). Notably, many other tools, including SIFT, PPH2_HD, PPH2_HV, PROVEAN, and PANTHER_PSEP, have achieved much higher Sen scores, but their Spe scores are very low, indicating that they are biased to positive predictions and prone to higher false positive rates (1-Spe). As MCC is a performance metric considering both positive and negative predictions in balance, we can conclude that AAS3D-RF achieved superior balanced performance without bias to a specific class. The detailed ROC curves demonstrate a similar conclusion: The curve of AAS3D-RF covers the largest area under the curve except FATHMM-W (Fig. 2C,D).

Fig. 2.

Comprehensive performance comparison of eight predictors on the two independent testing datasets. (A) and (B) show the performance in terms of seven metrics for TestDataset 1 and TestDataset 2, respectively. (C) and (D) display the ROC curves for TestDataset 1 and TestDataset 2, respectively.

Several previous studies have reported that the prediction performance of FATHMM-W was significantly confounded by type 2 circularity, i.e., it intended to predict all the variants from the same protein as pathogenic or neutral as a whole [31, 35, 36]. In other words, FATHMM-W would perform worse on a dataset with proteins containing a nearly equal number of daAASs and nAASs on each of them (e.g., the subsets with daAAS proportion in the range of [0.4, 0.6]) than on a dataset with proteins containing only daAASs or nAASs (e.g., the “Pure” subsets). Prediction methods sensitive to type 2 circularity would perform poorly in discriminating daAASs from nAASs within a given protein. In our work, the evaluations on “Pure” and different levels of “Mixed” subsets have demonstrated this: FATHMM-W performs the worst on the subsets with daAAS proportion in [0.4, 0.6], but obtains an impressively high AUC on the “Pure” subsets (Fig. 1C,D). As shown, AAS3D-RF consistently ranks at the top tier of all methods among “Mixed” subsets, indicating that the type 2 circularity in AAS3D-RF has been removed at the highest level. As their values are in fact the same for all variants from the same protein, the gene/protein-level features should be the source of type 2 circularity. In our cross-validation and independent testing, the variants from the same protein were either all in the training set or all in the testing sets, hence this grouped cross-validation at the protein level here may have served as an effective strategy to decrease type 2 circularity.

As published tools were trained on different datasets, it is often difficult to obtain a proper and fair testing dataset that contains enough data and has no overlap with any of their training datasets. Here, TestDataset 1 and TestDataset 2 are sufficiently rigorous to our predictor, since there are no overlapped AASs or proteins between them and TrainDataset. However, this scenario does not apply to other tools. Hence, the performance comparisons based on TestDataset 1 and TestDataset 2 will presumably provide a realistic performance estimate for AAS3D-RF but may supply overly optimistic estimates for other tools. In a truly fair comparison, the improvement of AAS3D-RF over other tools would presumably be larger.

To interrogate the interpretability of features, which may offer clues for understanding the molecular basis of the daAASs, we plotted the distributions of some features in daAAS and nAAS separately (Fig. 3). As shown, these features can be grouped within four categories: physicochemical properties of residues, conservation-related properties, gene/protein-level features, and structure-based features. While many of them have been analyzed in previous studies, we here mainly focus on several structure-based features.

Fig. 3.

Z-score distributions of a part of selected features for the daAASs (red) and nAASs (blue) in TrainDataset. The p-values of two-tailed Mann-Whitney test for each feature comparison are shown.

Five features in the conservation-related group show the most evident contrast between daAASs and nAASs (Fig. 3). Conservation often indicates structural stability or functional importance. As a comprehensive result of many underlying mechanisms, the conservation-related features have superior ability in separating daAASs from nAASs. Many previous studies have also repeatedly revealed the power of this type of features [34, 37, 38, 39, 40, 41, 42, 43]. However, these features cannot offer further molecular mechanistic clues of the disease-associated AASs, as mentioned in the Introduction.

Relative solvent accessibility (RSA) and folding free energy change ($\mathrm{\Delta }$$\mathrm{\Delta }$G) are two widely used features derived from 3D structures in discriminating daAASs from nAASs [11, 41, 44, 45, 46, 47, 48]. Fig. 3 shows that daAASs are more likely to be located at buried sites (small RSA values). The AASs occurring in buried sites may destabilize the protein by distorting the hydrophobic core, thus resulting in pathogenic effects [49, 50]. More directly than RSA, $\mathrm{\Delta }$$\mathrm{\Delta }$G measures the stability change caused by the substitution itself. In our dataset, the $\mathrm{\Delta }$$\mathrm{\Delta }$G values of daAASs are much larger than those of nAASs (Fig. 3), directly demonstrating that daAASs are more prone to large depletion of stability. This also suggests that stability loss resulting from AAS serves as a major mechanism of their disease-association.

The context-dependent substitution score (CDSS) is a substitution score between different amino acids under specific structural environments or contexts. Three features, including Miyata, Abs_dF1, and Abs_dpI, describe general physicochemical differences independent of residues’ structural environments. In contrast, CDSS provides a more specific residue substitution score in given RSA levels and secondary structure states [51]. In Fig. 3, the CDSS values of daAASs are smaller than those of nAASs, indicating that unfavorable substitutions with respect to CDSS tend to be associated with diseases.

Environment-dependent residue contact energy (ERCE) is another selected structural environment-related feature, accounting for both the secondary structure and residue contacts in the 3D structure [52]. The selected ERCE feature for an AAS describes the summation of contacting energy values between the wild-type residue and all its neighboring residues within a distance of 6.5 Å. The lower the contact energy, the more stability it contributes. In Fig. 3, the ERCE values of daAASs are smaller than those of nAASs, suggesting that the residues with higher stability contribution, if substituted, tend to become pathogenic.

To our knowledge, it is the first time that CDSS and ERCE have been utilized in developing methods for predicting the disease-association of AAS. The effectiveness of CDSS and ERCE indicates an understanding that substitutions introducing incompatible residues into a specific structural context will tend to be deleterious to the protein and thus be pathogenic.

In addition to better interpretability of the structural features, to what extent they contribute to the prediction performance is another aspect that must be interrogated. Toward this end, we re-trained an additional predictor with all seven structural features removed according to the same procedure adopted in training AAS3D-RF, namely, the structure-removed predictor. Then, we compared its performance with AAS3D-RF on TestDataset 1 and TestDataset 2 (Table 2 and Supplementary Table 3).

First, the ACC and MCC of structure-removed predictor were 2.8 and 4.9 percentage points less than the AAS3D-RF on the TestDataset 1, respectively (Table 2 and Supplementary Table 3), demonstrating that the structure features evidently contributed to the performance.

Second, when evaluated on TestDataset 2, the observed improvement of ACC and MCC after adding structure-related features was not as large as that observed on TestDataset 1 ($\mathrm{\Delta }$ACC = 0.004 from 0.835 to 0.839, $\mathrm{\Delta }$MCC = 0.008 from 0.676 to 0.684) (Table 2 and Supplementary Table 3).

What factors have affected the extent of performance improvement of structural features? As several previous studies have proposed that the contribution of structural features is more evident when reliable conservation-related features are unavailable [53, 54], we checked whether this applies in our case. On TestDataset 2, by using the number of aligned sequences (the feature of nal_1e-45) as a proxy of reliability of conservation-related features, we observed that the improvement of MCC was indeed more evident in AAS data with nal_1e-45 200 than those $\ge$200 ($\mathrm{\Delta }$MCC = 0.147 from 0.645 to 0.792 vs. $\mathrm{\Delta }$MCC = 0.007 from 0.676 to 0.683) (Supplementary Table 5). Similar results were obtained on TestDataset 1 (0.244 vs. 0.041) (Supplementary Table 5). The importance of structural features in predicting pathogenic AASs was also frequently emphasized in previous studies, since they can improve the prediction of nAASs in conserved/constrained regions and daAASs in regions with loose constraints, in addition to offering hints of pathogenic mechanisms [55, 56, 57].

In summary, incorporation of structural features further improves prediction performance, especially in scenarios wherein conservation-related features are of low reliability.