Use of 3-Way Voting of Machine Learning Algorithms Improves Prediction Performance of the Efficacy of Antisense-Mediated Exon Skipping and Reduces the Computational Burden

1. Introduction

Antisense oligonucleotides (ASOs) are small single-stranded nucleotides that target specific mRNAs by binding to their sense strand through Watson-Crick base pairing, which can be employed to modulate gene expression through various mechanisms [1]. The therapeutic potential of ASOs was recognized in the 1970s [2]. However, unmodified ASOs have limited plasma persistence [3]. ASOs have gone through three generations, with improved stability and binding affinity due to modified sugar moieties, bases, and phosphodiester linkages [4]. For example, 2’-O-methyl nucleotides (2OMe) and phosphorodiamidate morpholino oligomers (PMOs) are 2nd and 3rd generation ASOs, respectively [4].

ASOs modify target mRNA expression through two main mechanisms: RNase H-dependent cleavage and steric block [5]. RNase H-dependent ASOs, designed as gapmers, bind to the target RNA and trigger cleavage by the endogenous RNase H enzyme, leading to target gene silencing [6,7,8].. Steric blocking ASOs, on the other hand, are often employed to specifically exclude (exon skipping) or retain (exon inclusion) a specific exon(s), leading to alternations in splicing decisions [2,9].

Exon skipping, where an ASO causes the exclusion of a specific exon in splicing, has emerged as a promising treatment for genetic diseases, especially muscular dystrophies. US Food and Drug Administration has approved multiple exon-skipping ASO treatments for Duchenne muscular dystrophy (DMD), including eteplirsen, golodirsen, viltolarsen, and casimersen [10,11,12,13]. Exon skipping has shown promising potential as a treatment option for many genetic diseases beyond DMD. Splicing defects are a common cause of many genetic diseases, and exon skipping can be used to restore proper splicing by skipping over faulty exons. Milasen, a patient-customized n-of-1 ASO drug targeted for a pseudoexon in the CLN7 gene, was recently approved by the FDA for the treatment of Batten’s Disease, demonstrating the potential of exon skipping for personalized medicine [14,15]. Exon skipping therapies are also being explored for other genetic diseases such as cystic fibrosis, retinitis pigmentosa, sarcoglycanopathy, dysferlinopathy, fibrodysplasia ossificans progressiva, epidermolysis bullosa, frontotemporal dementia and parkinsonism linked to chromosome 17 (FTDP-17), and cancer, among others [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34].

Despite these promising developments, there are still significant challenges in developing effective exon-skipping therapies. A major hurdle is a difficulty in selecting an optimal sequence for exon skipping, as the efficacy of ASOs is often unpredictable due to numerous factors involved in the exon-skipping process [35]. Designing effective ASO sequences requires consideration of various criteria [36], particularly for exon skipping [37]. Software tools such as eSkip-Finder can aid in this process [38]. eSkip-Finder (https://eskip-finder.org) is a web-based tool developed by Chiba et al. that provides a solution for identifying optimal ASO sequences for exon skipping by using machine learning models built from a curated database of publications and patents [38].

The selection of important features is a crucial step in the tool’s approach, and the eSkip-Finder uses an exhaustive search of subsets of features to identify these critical components. However, due to the high computational cost, the subset size was limited to seven features. To optimize the performance of the models, hyperparameters in the support vector regressor are optimized through a grid search. This optimization process is computationally intensive, requiring a significant amount of computing power, and can take several days to complete.

This paper seeks an alternative solution to reduce the computational cost associated with the eSkip-Finder. Some machine learning algorithms such as decision-tree or random forest have built-in feature ranking capabilities [39]. Ensemble methods are also proven to have good performance with reasonable computation cost [40,41]. We explored their utility in ASO efficacy prediction and demonstrated that a combination of three algorithms, namely random forest, gradient boosting, and XGBoost, through a 3-way voting mechanism can significantly reduce computation time while maintaining or slightly improving the prediction performance. This approach offers a promising solution for reducing computational cost in the ASO efficacy prediction process.

2. Materials and Methods

The datasets used in this study were the same as those used in Chiba et al. [38]. That is, for PMO, 369 and 57 measurements were used for training and testing and there were 98 and 11 unique ASO sequences in each split without overlapping; for 2OMe, 197 and 31 measurements were used for training and testing and there were 111 and 13 unique ASO sequences in each split without overlapping. As PMO and 2OMe have different chemistry thus different binding affinity, the datasets were handled separately.

For each measurement, there were 32 numerical features calculated via bioinformatics tools as discussed in Chiba et al. (such as dose). The categorical feature, Malueka’s category, was excluded from modeling. As reported in [38]. this feature is not important in determining the ASO efficacy. The feature was specifically linked to dystrophin exons [42]. Models developed with this feature included will be difficult to generalize to other genes.

The efficacy was measured as a percent in the range 0 to 100. The efficacy is the value to predict, making this a regression problem. All 32 features were inputted into the machine learning models and feature selection was left to the models.

The machine learning libraries included scikit-learn (0.42.2) [43] and XGBoost (1.6.1) [44]. The following regressors were used: support vector, random forest, gradient boosting, and XGBoost. The last three were also used to vote by the simple average of the individual predictions. The support vector regressor was included for comparison purpose, as it was used in Chiba et al. All those regressors were built without hyperparameter tuning, i.e., default parameters were used in each regressor (except random seeds). The computation code was developed using Python (3.9.7) on Mac (Quadcore i5, 2 GHz CPU, 16 GB RAM).

Two metrics were used to assess model performances: R² and mean absolute error (MAE) between true efficacy values and predictions. The models were first assessed on the training data via 10-fold cross-validation. The best model was then selected and applied to the reserved test data. The R² and MAE on each fold were collected and their mean and standard deviation were further computed to aid the best model selection.

While the random forest, gradient boosting, and XGBoost models were trained, they also collected data to compute the feature importance score. The voting regressor had no feature importance score, however. We therefore used the model-agnostic method, permutation feature importance provided by scikit-learn, to rank the feature importance.

3. Results

The performance metrics for various models using 10-fold cross-validation on the training data are shown in Table 1. 5-fold and 20-fold cross-validations were also attempted and the results were similar to what was reported here. The data splitting was based on ASOs, i.e., there were no overlapping ASOs in training and validation splits. As can be seen from Table 1, for both PMO and 2OMe ASOs, the 3-way voting approach gives the largest R² and smallest mean absolute error (MAE). We thus chose this approach and applied it to the test datasets. The support vector regressor performed noticeably poorly as there was no hyperparameter optimization in the current study. It shall also be noted that the whole computing took about 10 seconds on a laptop computer.

When the 3-way voting models, trained on the training data with all features, were applied to the test data, the predictions were similarly assessed. For PMO, we have R² = 0.706 and MAE = 12.25, and for 2OMe, R² = 0.795 and MAE = 9.237. The R² values are higher than those reported [38], which were 0.6 and 0.7 respectively. The true efficacy and predicted one have a good linear correlation, as depicted in Figure 1. It shall be noted that, unlike the support vector regressor which can generate unrealistic, negative efficacy values, the 3-way voting approach will not possibly predict a negative efficacy as long as the input data has no negative efficacy.

The feature importance ranking using the training data as reported by the 3-way voting is shown in Figure 2. The rankings using the test data are similar on top-ranked features, suggesting that overfitting is not a concern. Among top 5 and 10 features using training or test dataset, 3 and 8 are common for PMO and 4 and 9 are common for 2OMe. The 4 PMO features used in Chiba et al. here were ranked at 1, 24, 11, and 15. The 6 2OMe features used in Chiba et al. here were ranked at 2, 25, 4, 3, 17, and 11. In both cases, some correlation can be observed. We also noted that some features were strongly correlated, e.g., niscore and niscore_per_base. Niscore_per_base was ranked 17th, but niscore was ranked 5th in our 2OMe model. Therefore, at least some discrepancies can be attributed to the feature correlations. Due to the randomness in the algorithms, the rank order can be slightly different in each run.

To check if the voting approach works for different genes and exons, we applied the trained PMO model to the exon 73 skipping of collagen type VII alpha 1 chain [9]. The results are summarized in Table 2. The predictions by the voting approach preserve the ranking order of ASO efficacy experimentally measured. Cautions must be taken when one extends the model to a different application domain, however. As more data is accumulated in databases such as eSkip-Finder, we expect predictive models will be validated rigorously and extended as needed.

4. Discussion

We applied machine learning algorithms with built-in feature selection capabilities to train on and predict exon-skipping PMO and 2OMe ASO efficacy. The model build process requires much less time. Among various algorithms assessed, the voting strategy yielded the best-performing predictors in terms of R² and mean absolute error (MAE) between the true and the predicted efficacy. R² were 0.706 (PMO) and 0.795 (2OMe), which were slightly higher than were reported [38]. The MAE was also reported as a reference. This observation on the voting approach was consistent with the general consensus in the machine learning community. Due to the model itself, no negative efficacies are predicted in our approach, whilst the support vector regressor does not have this guarantee. Important features used in our approach were similar to what eSkip-Finder discovered. Features used by Chiba and colleagues were overall ranked high in our voting approach and some differences can be explained by feature correlations. Thus, our modeling approach has similar interpretability.

As mentioned above, the voting approach predicts non-negative efficacies as long as there are no samples with negative efficacies in the training data. However, this can be a drawback, i.e., the approach will not predict any efficacies larger than the highest efficacy in the training data, since decision trees are used essentially in the individual algorithms. This potential limitation can be easily remedied by collecting training samples with large efficacies.

The proposed voting approach has a very short training time. We believe that the same approach might be applicable to developing predictive models for other diseases where ASO efficacy data is available. The voting scheme still relies on engineered features scientists hand-picked. As a possible future extension, one could consider machine learning algorithms in combination with natural language processing techniques, which has been successfully applied to biological sequence analysis [45].