RNA-seq FPKM (Fragments Per Kilobase per Million) profiles and clinical information about the status of lymph nodes of the TCGA prostate adenocarcinoma (PRAD) cohort were acquired from the TCGA data portal (https://portal.gdc.cancer.gov). The corresponding RNA-seq splicing events profiles of PRAD were obtained from the TCGA SpliceSeq (16). In order to get reliable AS event data, we adopted a rigorous screening filter criteria with a sample proportion of PSI values of no less than 75% (17). The cases without the status of lymph node or lack of matched RNA-seq splicing events profiles were excluded. Finally, there were 413 cases included in our analysis cohort.

The overview of the workflow is shown in Figure 1 . The PRAD cohort was divided into two groups by the presence or absence of lymph node metastasis. We used the limma package (18) for differential analysis. An adjusted p value < 0.05 was applied as the threshold to determine the DEAS.

For the purpose of evaluating and analyzing the whole dataset structure. We performed unsupervised clustering on the DEAS feature pool by k-means clustering. Cluster analysis was performed in order to determine LNM patients clustering patterns without knowing the results in advance and evaluated by comparing the clustering outcome with the underlying facts.

The dataset was split randomly into a training set and a validation set, containing 351 cases in the training set and 62 cases in the validation set ( Supplementary Table 1 ). To lessen feature redundancy, two feature selection methods were used in the training set to identify features that might be essential for our model. In the first step, a correlation filter premised on the absolute values of pairwise Spearman’s correlation coefficient was applied. The threshold was set at 0.8 for ρ. In a nutshell, if two features have ρ > 0.8, the function examines the mean absolute correlation of each feature, and the feature with the higher mean absolute correlation will be eliminated. In the second step, the least absolute shrinkage and selection operator (LASSO) logistic regression algorithm was applied to choose the most optimized predictive features from the features selected from correlation filter algorithm.

In this study, we built eight machine learning classifiers, namely random forest (RF), multi-layer perceptron (MLP), logistic regression (LR), gaussian naive bayes (GNB), linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), support vector machine (SVM), and light gradient boosting machine (LGBM). The eight machine learning classifiers were trained on the training set, respectively. To evaluate the performance and overall error estimation of the eight classifiers, we applied a stratified five-fold cross-validation method with 100 iterations. Oversampling method was not applied in this study, as the ratio of positive samples to negative samples is approximately equal to 1:5. In the training set, each fold in turn was used as a validation set, and the other four folds were used as a training set. The validation outcomes from 100 five-fold cross-validation were then integrated to present a measure of global performance.

The statistical analysis was performed with the R (version 4.1.2) and the Python (3.9.7) programming languages and environments with the scikit-learn repository. The LGBM classifier was built on the lightgbm module. The performance of the classifiers was assessed by area under the receiver operating characteristic curve (AUC), accuracy, recall/sensitivity, specificity,and F1-score. P value of less than 0.05 indicates statistical significance.

SHAP (SHapley Additive exPlanations) was performed to explain the model. SHAP is a game theoretic method to explain the output of any machine learning model (19). SHAP values evaluate the significance of the output resulting from the inclusion of one feature for all combinations of other features.

Finally, there were 413 cases included in our analysis cohort, with 336 lymph node negative cases and 77 lymph node positive cases. 44070 AS events were preliminary identified from 10381 genes, including 3524 AA events in 2488 genes, 3101 AD events in 2185 genes, 9035 AP events in 3621 genes, 8663 AT events in 3781 genes, 16772 ES events in 6578 genes, 228 ME events in 221 genes, 2747 RI events in 1849 genes ( Supplementary Figure 1 ).

In this work, 333 important DEAS events in 255 genes were identified by the limma package as being correlated with lymph node metastasis by comparing LNM and non-LNM cases ( Supplementary Table 2 ). There were 157 ascending DEAS events and 177 descending events. Compared with the proportion in all AS types, the proportion of AT type in DEAS increased more apparently than any other type, which indicated that AT type may perform a crucial role in tumor metastasis ( Supplementary Figure 1 ).

Principal component analysis as well as k-means clustering was used to analyze the whole DEAS feature stream and two clusters were detected. A 76.6% compactness, the degree to which group members share similarities, was detected inside the clusters. The established cluster was validated by applying the silhouette coefficient (silhouette width), which is an algorithm to evaluate the cluster results. The silhouette plot showed that the clustering using two groups was perfect, with no negative silhouette width and the majority of cluster values greater than 0.03, as shown in Figure 2A . The lymph node positivity rate was 0.074 in cluster 1 and 0.347 in cluster 2. 76.6% of LNM were collocated in cluster 2 according to the unsupervised clustering algorithm ( Figure 2B ). Most of LNM patients appeared to be clustered together, potentially suggesting that DEAS is a reliable data source to distinguish between LNM and non-LNM.

First, 268 key features from 333 DEAS were acquired by using pairwise Spearman’s correlation filter with a 0.8 correlation coefficient threshold in the training set. Second, applying the LASSO-based feature selection method ( Figures 3A, B ), we further selected 96 crucial features from 268 key DEAS ( Supplementary Table 3 ).

The final selected 96 features were then treated as input layers in the subsequent classifiers. With stratified five-fold cross-validation on the training set, we ran 100 iterations on each of the eight distinct classifiers to evaluate their performance. Performance of the eight classifiers was shown in Supplementary Table 4 . Due to unbalanced data, we selected the optimal classifier based on F1-score and finally the LDA model was chosen with the best F1-score performance. The performance outcome showed that the LDA algorithm achieved an average AUC of 0.962± 0.026 and an accuracy of 0.929± 0.028, a specificity of 0.958± 0.024, a sensitivity of 0.812± 0.111, a F1-socre of 0.809± 0.079 on the stratified five-fold cross-validation set

Subsequently, the independent validation set was tested with the same LDA classifier. An AUC of 0.953 ( Figure 4A ) was reached by the classifier. The F1-score, specificity, sensitivity, and accuracy were all observed to be 0.815, 0.979, 0.917, and 0.919 respectively. In the validation set, eleven LNM and forty-six non-LNM cases were accurately detected. In the remaining cases, one LNM and four non-LNM cases were incorrectly categorized as non-LNM and LNM, respectively ( Figure 4B ).

In the validation set, there was a strong agreement between the observed LNM rate and the model prediction ( Figure 4C ) revealing good discrimination of the classifier. Hence, our model performed well in the internal testing set. The decision curve analysis ( Figure 4D ) demonstrated that the application of the LDA model to predict LNM in the validation set indicated a greater net benefit increase than the “treat everyone” or “no treatment” strategic scheme over a wide range of threshold probabilities, showing the model’s utility in clinical settings.

We used SHAP to explain the significance of each feature to the model output. Summary plot was drawn to display the top 20 features, which had the most impact on the model output. Figures 5A, B show how high and low the feature values were relative to the SHAP values in the training set. The features were listed from top to bottom in descending order by magnitude of impact on model output, with the first feature having the greatest influence. The parameter values of each feature variable are represented in color on the right side of the variable name, with red representing the high parameter value and blue representing the low parameter value. The higher the SHAP values, the higher the probability of lymph node metastasis. The reverse applies when the SHAP values decrease. The feature with the highest value is CALCOCO1|22108|RI. The lower feature value of CALCOCO1|22108|RI, the higher probability of lymph node metastasis, indicating protecting role of CALCOCO1|22108|RI in LNM.