Var3PPred: variant prediction based on 3-D structure and sequence analyses of protein-protein interactions on autoinflammatory diseases

Dataset construction

A total of 35 genes pertinent to autoinflammatory diseases were listed along with their respective Online Mendelian Inheritance in Man (OMIM) accession number of associated diseases. A total of 702 missense variations found in these genes were collected from the Infevers database (Van Gijn et al., 2018). Interaction partners of these genes were obtained from the STRING and Intact databases using Cytoscape (Shannon et al., 2003; Szklarczyk et al., 2019; Hermjakob et al., 2004). A total of 191 distinct interacting partners of 35 autoinflammatory disease-related proteins were identified. The protein sequences of 226 genes, which include the 35 autoinflammatory disease-related proteins and their 191 interacting partners, were retrieved from UniProt (Consortium, 2019). 3D structures of these 226 proteins were modeled by trRosetta and AF2 (Yang et al., 2020; Jumper et al., 2021). Pathogenic and benign variations in the 35 autoinflammatory genes were analyzed for their ZDOCK and SPRINT interaction scores, weighted with 191 interaction partners sourced from the STRING and IntAct databases (Szklarczyk et al., 2019; Hermjakob et al., 2004). 3D structures of the proteins were used to obtain the protein-protein complexes by ZDOCK. Sequence-based PPI scores were obtained by SPRINT tools. An enhanced ΔΔG score has been computed using the PyRosetta Python package’s weighted score function (Chaudhury, Lyskov & Gray, 2010), in conjunction with the importance scores of HGPEC genes across 35 autoinflammatory diseases. For the 702 missense variations identified, scores for ASA, PSIC, pLDDT, and volume of variation were compiled. To overcome the imbalance between pathogenic and benign variations, the SMOTE algorithm was used to produce models with higher prediction rates by balancing the ratio between benign and pathogenic variants. The workflow of Var3PPred was presented in Fig. 1.

Protein structure modelling

The three-dimensional structures of SAID-associated proteins and interaction partners were modeled using trRosetta (Yang et al., 2020). The initial step of 3D structure modeling involved performing multiple sequence alignments (MSAs) for each protein sequence utilizing the HHblits tool with the Pfam database (Remmert et al., 2012; El-Gebali et al., 2019). These MSAs serve as crucial input for the subsequent trRosetta-based structure prediction, providing evolutionary information that aids in defining the structural constraints.

Analysis of variations’ protein structure and motifs: ΔΔG, pLDDT, ASA, and PSIC scores

The process of predicting the protein structure entails a deep residual network-based prediction of inter-residue orientations and distances, which guide the folding of the protein into its final 3D form. Furthermore, to assess the structural implications of 702 variations on the 3D structure of corresponding proteins, they were modified using PyRosetta (Chaudhury, Lyskov & Gray, 2010). Also, the change in the Gibbs free energy (ΔΔG) due to the variation was calculated with PyRosetta’s ΔΔG module.

From the AlphaFold Protein Structure Database, we extracted the AF2 computed structures of 35 proteins. Using these structures, we obtained/calculated the pLDDT scores, which quantify the local protein structure prediction confidence and the relative surface accessible surface area by the DSSP (Dictionary of Secondary Structure of Proteins) module available in the Biopython package (Cock et al., 2009). The PSIC scores for both the wild-type and variant amino acids were obtained from the pre-computed profiles of PolyPhen database. These scores provide insights into the possible impacts of amino acid substitutions on protein function and stability.

Handling missing data involves a systematic imputation strategy where each feature exhibiting missing values is modeled as a dependent function of other available features, employing a sequential, round-robin methodology. This process is operationalized using the fancyimputer package in Python, which facilitates advanced statistical techniques for effective and accurate data imputation (Rubinsteyn & Feldman, 2016). The imputation strategy employed is IterativeImputation, a multivariate imputer that estimates each feature from all the others. This approach imputes missing values by modeling each feature with missing values as a function of other features in a round-robin fashion.

3D and sequence-based protein–protein interaction analysis

Structural protein files in pdb format obtained from the trRosetta tool were preliminary for the ZDOCK tool (Pierce, Hourai & Weng, 2011). Structural protein files in pdb format received from the trRosetta tool were preliminary for the ZDOCK tool. Surface amino acids of proteins are marked, and electrostatic charges of atoms are calculated according to UNICHARMM data. The GNU-parallel tool is used to run processes in parallel on high-performance computers (Tange, 2018). Interaction partners of SAID-associated genes obtained from STRING and Inact databases were examined with ZDOCK, and the highest score was obtained. Proteins’ sequences in fasta format are accepted in the prepared pdb files. Similar sub-sequences are obtained using the PAM120 scoring matrix. Then, interaction scores of mutant interactions are obtained by trained interactions (Li & Ilie, 2017). Pairwise interaction scores of Protein interactions ZDOCK and SPRINT are shown in Table S2.

Weighting interaction scores with the HGPEC gene rank score in SAIDs

Random walk with restart algorithm on a heterogeneous network (RWRH) is used to get the genes’ SAIDs disease relatedness score. (1)${\displaystyle P^{t+1}=\left(1-\gamma \right)W^{\prime }{P}^t+\gamma P^0}$

The RWRH formula is represented at Eq. (1) as Le & Pham (2017) described, where P^t is the probability vector at iteration tt, γ is the restart probability, and W is the transition probability matrix of the network. The matrix W is constructed by combining the gene-disease and disease-gene interaction probabilities, normalized to ensure that each row sums to one. The initial probability vector P0 is set based on the diseases of interest, with non-zero probabilities assigned to the corresponding nodes in the network. The RWRH algorithm iteratively updates the probability vector until convergence, with the final probabilities indicating the relatedness of each gene to the specified diseases. By tuning the parameter γ, researchers can control the balance between exploration of the network and focus on the initial disease nodes, allowing for the identification of both closely and distantly related genes.

The pre-existing ‘Disease Similarity Network 15’, integrated within the HGPEC toolset, is a pivotal resource for discerning disease similarity in SAIDs. The prioritization process ranks genes based on their relevance in the recommended human Protein-Protein Interaction (PPI) network and similar selected diseases incorporated in the HGPEC (Le & Pham, 2017). These rank scores will later be used in weighting ZDOCK and SPRINT data, as shown in Table S3.

SPRINT and ZDOCK interaction scores were then calculated for all interaction partners of SAID-related genes corresponding to a single mutation, yielding a comprehensive interaction score.

The total score for a mutation was calculated by averaging all weighted interaction scores associated with the mutated protein, as depicted by Eq. (2): (2)${\displaystyle W_v=\frac{\sum _{}^{}\left[\mathrm{G}-{\mathrm{V}}_i\right]{\mathrm{V\; r}}_i}{{\mathrm{V\; r}}_{\mathrm{sum}}}.}$

G is the interaction score of the mutated protein with its interaction partner in the network. V _i is the interaction score of a wild-type protein with its interaction partner in the network. V r_i is the HGPEC gene rank score of the interaction partner of SAID-related gene. V r_sum is the sum of all V r_i values for all interactions of the mutated protein.

The overall fraction summarizes all the individual weighted interaction scores for the mutated protein. It is divided by the sum of all HGPEC gene rank scores of interaction partners of SAID-related genes. This fraction provides the average weighted interaction score for the mutated protein across all its interactions.

Machine learning models training and evaluation

In selecting these six machine learning algorithms for our study, we aimed to encompass a diverse range of approaches to evaluate their performance on complex biological data. Random forest and AdaBoost, both tree-based ensemble methods, were chosen for their proven effectiveness in handling high-dimensional data and their ability to improve predictive accuracy through ensemble learning. Linear SVM and logistic regression represent linear models, which, despite their simplicity, are powerful tools for classification tasks and provide a baseline for comparison with more complex models. KNN, a non-parametric method, was included for its simplicity and intuitive approach to classification based on proximity in the feature space. Lastly, QDA, a parametric discriminant analysis technique, was selected for its ability to model complex decision boundaries when the data is assumed to follow a Gaussian distribution. Together, these algorithms offer a comprehensive overview of the performance of different machine learning approaches on biological data, with the Python Scikit-learn package providing a consistent and efficient framework for their implementation and evaluation (Pedregosa et al., 2011).

Before training models, the imbalance between benign (130) and pathogenic (572) variants number in the dataset was balanced by employing the Synthetic Minority Over-sampling Technique (SMOTE) (Chawla et al., 2002).

These models were trained on a dataset of 80% missense variations dataset, with respective ΔΔG values, weighted SPRINT scores, ZDOCK scores, pLDDT scores, ASA, mutant amino-acid volumes, as well as wild-type and mutant PSIC scores. Subsequently, 20% of the dataset was partitioned off for model testing. Each model’s hyper-parameters were optimized according to their performance on the training set. The explored hyper-parameters for each model included:

• ‘Random forest’: Number of estimators (700, 800, 900,1000), maximum depth of the trees (None, 5, 10, 20, 30), minimum samples required to split an internal node (2, 5, 10), and minimum samples required at each leaf node (1, 2, 4).

• ‘Linear SVM’ and ‘Logistic Regression’: Regularization parameter (0.1, 1, 10, 100).

• ‘K-nearest neighbor’: Number of neighbors to use (3, 5, 10, 15).

• ‘AdaBoost’: Number of estimators (10, 50, 100, 200).

Following training, models were evaluated on the test data using metrics, including Receiver Operating Characteristic Area Under Curve (ROC AUC), F1 score, balanced accuracy, accuracy, and Matthews correlation coefficient (MCC). Moreover, the models were subjected to rigorous validation via 20-fold cross-validation on the training data, with performance assessed using the same metrics as applied to the test predictions.

We compared various features and their interaction pairs within the context of the previously mentioned machine learning algorithms. To optimize performance, hyperparameter tuning was conducted for each of these models.

Feature importance, decision-making progress, and features paired interaction impact on model decision of selected models were evaluated with SHapley Additive exPlanations (SHAP) algorithm (Lundberg & Lee, 2017).

An in-depth analysis of the features and their interactions was conducted by generating scatter plots for paired features, providing a comprehensive visual representation of the data distribution. To delineate the decision boundaries established by the predictive models, 2D contour plots were produced. Thereby allowing an intuitive understanding of the models’ decision-making patterns. These models are also trained in 80% of the balanced dataset. The performance of each model was evaluated using the remaining 20% of the data, which was reserved for testing.

Comparison of variants features

Examined ZDOCK, SPRINT, ΔΔG, ASA, and volume change of variations; particularly, benign variations show higher means and standard deviations for ZDOCK and SPRINT values than pathogenic variations. Interestingly, the scenario is reversed for ΔΔG values, where pathogenic variations exhibit higher distribution and mean values than benign ones. This differential pattern can be visualized effectively via scatter plots, as shown in Fig. 3. Specifically, pathogenic variations present a broader ΔΔG distribution on the y-axis, whereas benign variations show a more concentrated distribution. Conversely, the ZDOCK data offer a wider range for benign variations. Moreover, the pLDDT score shows comparable separability to SPRINT in distinguishing between benign and pathogenic variants. Additionally, the wild-type PSIC and mutant PSIC scores, indicative of conservation values, are negatively correlated, further emphasizing their relevance in variant segregation.

It is also worth mentioning that pathogenic variants show broader separation in the ASA values. Although the mean Volume of variations is similar for benign and pathogenic variants, this feature demonstrates better separation of variants in combination with the pLDDT value. This suggests the potential utility of incorporating these diverse features in comprehensive variant classification frameworks. ZDOCK and SPRINT have a positive correlation of 0.089. This might indicate that variations have a more substantial impact on protein-protein interactions and also tend to have a stronger impact on domain-domain interactions. ΔΔG and the Volume of the mutation have a positive correlation of 0.14, suggesting that variations causing more extensive changes in residue volume might also cause more significant changes in protein stability. SPRINT and pLDDT have a strong negative correlation of −0.54, suggesting that variations affecting domain-domain interactions tend to occur in regions of lower structural confidence. ASA and pLDDT have a strong negative correlation of −0.71, indicating that residues with higher structural confidence tend to be less exposed. ASA and SPRINT have a positive correlation of 0.462608, suggesting that variations in more exposed residues might have a stronger impact on domain-domain interactions. Wild-type PSIC and ASA have a strong negative correlation of −0.49, suggesting that more conserved positions tend to be less exposed. According to Fig. 3, the most separable pair-wised features are pLDDT and volume of mutation, and SPRINT features separate benign and pathogenic variants when paired with all other features.

Machine learning models

Comparison results with other models are given in Table 1. The best model is the random forest model algorithm according to ROC AUC and balanced accuracy metrics. The comparative analysis of machine learning algorithms for predicting SAIDs associated variations reveals significant variations in the performance of all the tested models. Tree-based methods such as random forest and AdaBoost displayed superior predictive power across multiple evaluation metrics, particularly regarding the ROC AUC, accuracy, and F1 score, achieving scores of approximately 0.995, 0.978, and 0.95 for both models, respectively. Random Forest, an ensemble learning method that operates by constructing multiple decision trees, demonstrated the highest MCC of 0.90, indicating its robust reliability and strong binary classification abilities. Similarly, the AdaBoost algorithm, known for focusing on classification problems where instances are weighted, also presented a high performance with an MCC of around 0.89.

Meanwhile, other models, including linear support vector machine (SVM), logistic regression, K-nearest neighbour, and quadratic discriminant analysis, yielded relatively lower performance metrics, with the ROC AUC ranging from 0.90 to 0.92 and MCC from 0.51 to 0.77. These differences underline the effectiveness of ensemble tree-based methods in handling the complexity and high dimensionality of biological data, yielding promising results for predicting variants. The random forest algorithm has been selected for subsequent investigations based on the performance metrics across multiple machine learning models.

The superior performance of tree-based methods such as random forest and AdaBoost in our study can be attributed to their inherent ability to handle the complexity and high dimensionality of biological data. These methods effectively capture non-linear relationships and interactions between features, which are common in biological datasets. The ensemble approach of random forest, which combines multiple decision trees, helps in reducing overfitting and improving the generalization of the model. Each tree in the ensemble is built on a random subset of the data and features, leading to diverse trees that capture different aspects of the data. This diversity helps in achieving high accuracy and robustness in the predictions.

AdaBoost, on the other hand, focuses on instances that are difficult to classify, iteratively adjusting the weights of these instances to improve the model’s performance. This adaptive weighting scheme allows AdaBoost to concentrate on the areas of the data where the model’s performance is weak, leading to a more refined and accurate classification (Tang, Henderson & Gardner, 2021).

In contrast, other models like linear SVM, logistic regression, K-nearest neighbor, and quadratic discriminant analysis may not be as effective in capturing the complex patterns present in biological data. These models have limitations in modeling non-linear relationships and may struggle with the high dimensionality of the data, leading to lower performance metrics compared to tree-based methods. The linear nature of SVM and logistic regression, the distance-based approach of KNN, and the assumption of a Gaussian distribution in QDA may not be well-suited for the intricate structure of biological data, highlighting the advantage of using ensemble tree-based methods for such tasks (Singh, 2019; Dreiseitl et al., 2001; Kim et al., 2011).

In addition to the models mentioned, we also evaluated a basic shallow fully connected neural network (multilayer perceptron), which achieved an accuracy of approximately 0.85. While this performance is commendable, it still falls short of the 0.98 ROC AUC achieved by the random forest model. Furthermore, a single decision tree model exhibited lower accuracy compared to the random forest, with its performance being sensitive to the random seed used during model creation, indicating variability in its predictive capability.

The most important figures are the wild-type and mutant PSIC scores. The weighted SPRINT interaction score is the third most important feature in the random forest model. Also SHAP impact score of the feature correlated with feature importance (Figs. 4D, 4E). However, in the decision plot, which was employed with 20% test data, some of the test variant’s prediction scores were determined with other features like ASA, pLDDT, and ZDOCK (Fig. 4B).

According to the results of this model, the ROC AUC value is 99%, according to the average of 20 split cross-validations. Other metrics are given below.

• F1 test result: 0.990280

• Precision: 0.993148

• Recall: 0.987684

• Accuracy: 0.956031

• ROC AUC: 0.9946581

• Balanced Accuracy: 0.990332

• MCC: 0.981044

In a comprehensive evaluation across six machine learning models, the pairing of SPRINT and wild-type PSIC stood out, achieving an accuracy of 72% in random forest, linear SVM, KNN, and AdaBoost. This combination also exhibited 71% and 68% accuracy in logistic regression and QDA, respectively. Other notable pairings include pLDDT with wild-type PSIC (71% accuracy) and pLDDT with SPRINT (70% accuracy). Interestingly, the ZDOCK and SPRINT combination was examined in two scenarios, registering accuracies of 73% and 71%. Additionally, combinations like wild-type PSIC with ASA, pLDDT with ASA, ΔΔG with SPRINT, and SPRINT with ASA achieved accuracies hovering around 67% to 70%. The least-performing pairing was ZDOCK with ASA, which reached 61%. The results emphasize the importance of feature pairing in enhancing predictive accuracy in machine learning models (Fig. 5). The learning curve analysis reveals convergence between test and training data, suggesting proper model fitting Fig. S2.

The train and prediction regions, according to the dual selections of the features of the created models, are shown in Fig. 5. Here, blue regions are determined to predict benign variations, and red areas are determined to predict pathogenic variations. The locations of the dual feature data in these regions are seen as a scatter plot.

Comparison with other prediction methods (SIFT, Polyphen, CADD)

A comparative analysis was conducted to benchmark the predictive capabilities of our model against three widely used variation prediction methods: SIFT, PolyPhen, and CADD. The test set for this comparison included 59 variations from exons 2 and 10 of the MEFV gene.

Our analysis revealed that our model outperforms these established methods in its predictive capacity. The SIFT model yielded a higher number of false positives (five) and false negatives (13), thereby suggesting a lower accuracy. Similarly, the PolyPhen tool, while showing no false positives, had 11 false negatives. The CADD method, which classified variations with a Phred score above ten as pathogenic, resulted in eight false positives and nine false negatives. In stark contrast, our model demonstrated superior performance with eight false positives and, notably, no false negatives. This result underscores the robustness of our model, highlighting its potential to offer improved accuracy for predicting the pathogenicity of genetic variants in the context of SAIDs Fig. 6. Moreover, Var3PPred achieves a higher ROC AUC of 0.99 (Fig. 4C) compared to 0.92 for AlphaMissense, 0.92 for BayesDel, 0.91 for MetaLR, 0.90 for MetaSVM, 0.85 for DEOGEN2, 0.85 for ClinPred, 0.80 for DANN, 0.79 for FATHMM, and 0.75 for Eigen-raw (Fig. 6E).

The prediction results for 420 missense VUS in genes associated with SAIDs and associated domains are available on our GitHub page for further reference and analysis. Our prediction model evaluated VUS, 23.61%) as benign and 76.39% as pathogenic.