An interpretable machine learning model for predicting prognosis of medulloblastoma integrating genetic and clinical features

Machine Learning


Baseline clinical information

This retrospective study involved 729 patients with MB in the Chinese cohort for the identification of the prediction model. These 729 patients were allocated into separate training and testing sets (Fig. 1a). The international cohort consisting of 201 patients was used as the external validation set (Fig. 1b). Supplementary Tables 1 and 2 summarize the comparison of clinicopathological and molecular characteristics of patients who received postoperative radiotherapy and/or chemotherapy among the training, testing, and external validation sets. The details of study design are displayed in Fig. 2. The design consists of five parts: data preparation, model development, optimal model selection, model interpretation, and an online calculator construction.

Fig. 2: Overview of the study methodology.
Fig. 2: Overview of the study methodology.

This figure illustrates the five core phases of the research, ranging from data preparation to the construction of clinical application tools. a Data preparation: flowcharts detailing the screening process for the Chinese discovery cohort (n = 1043) and the international consortium external validation cohort (n = 471). Based on exclusion criteria, 729 Chinese patients and 201 international patients with complete clinical data were ultimately included. The Chinese cohort was further partitioned into training and testing sets at a 7:3 ratio. b Model development: integration of clinical features, molecular features, and radiotherapy parameters. Six algorithms were employed to construct predictive models across four distinct application scenarios: CMR, CM, CR, and CO. c Selecting the best model: evaluation of model performance for 5- and 10-year survival predictions via ROC curves, AUC, calibration curves, and DCA. d Model interpretation: utilization of the SHAP method to rank feature contributions and provide individualized explanations for the final selected model. e Online calculator: translation of the optimized models into interactive web-based Shiny applications that support the input of clinical and molecular features to dynamically assess patient survival probabilities at various time points.

Among the 729 patients in the derivation cohort, molecular subgroups were identified in 424 (58.2%) patients, including 40 WNT-MB, 84 SHH-MB, 84 Gr.3-MB, and 216 Gr.4-MB. Patients were followed up for a median of 6.9 years (95% CI, 6.6 to 7.4). The 5-year and 10-year cumulative OS rates were 82.7% and 67.5%, respectively. The median age at diagnosis was 8 years (IQR, 6 to 11). Based on the presence of dissemination on cytology or gadolinium-enhanced craniospinal MRI, 110 patients (15.1%) had metastases (M + ). Histopathological classification was performed for 525 (72.0%) patients. Of those, the proportion of CMB (41.8%) was the highest, followed by DNMB (20.6%) and LC/AMB (5.6%), and the lowest percentage was found in MBEN (4.0%). GTR/NTR was achieved in 665 (91.2%) patients and STR in 64 (8.8%). A total of 533 cases (73.1%) received postoperative radiotherapy combined with chemotherapy, while the other 196 cases (26.9%) had only radiotherapy after the operation. The median CSI dose was 30.6 Gy (IQR, 28.8 to 36.0) with PFTB boosted to a median dose of 55.8 Gy (IQR, 54.0 to 55.8). According to the classification of the observed outcomes, 729 cases were divided into alive and deceased groups, of which 542 cases were survivors. The clinicopathological and molecular characteristics of the training and testing sets are listed in Supplementary Table 1 and 2. There was no statistically significant difference between the two sets for all the analyzed characteristics (all p > 0.05).

Supplementary Tables 1 and 2 provide the clinicopathological and molecular characteristics of 201 patients in the international cohort, which had a median follow-up of 5.8 years (Q1-Q3: 5.3–6.7). The median age was 8 (5.0, 12.6) years, and 135 (67.2%) were male. Molecular subgroups were identified in 200 patients (99.5%), including WNT-MB (n = 4), SHH-MB (n = 61), Gr.3-MB (n = 45), and Gr.4-MB (n = 90). The M+ stage was observed in 60 (29.9%) patients and the M0 stage in 141 (70.2%). Histological subgroups were identified in 179 patients (89.1%), including CMB (n = 121), DNMB (n = 27), MBEN (n = 6) and LC/AMB (n = 25). GTR/NTR was achieved in 181 (90.0%) patients and STR in 20 (10.0%). A total of 154 cases (76.6%) received postoperative radiotherapy in combination with chemotherapy, while the other 47 cases (23.4%) received postoperative radiotherapy alone. The median CSI dose was 24.0 Gy (IQR, 23.4 to 36.0) with PFTB boosted to a median dose of 54.0 Gy (IQR, 54.0 to 55.8).

Feature selection

In the CMR and CM scenarios, candidate predictors were selected using univariate Cox analysis (Fig. 3a), based on NCCN guidelines, previous studies, clinical expertise, and data accessibility. The expression levels of molecular markers MYC, MYCN, OTX2, and GFI1 (p = 0.026) were ultimately included in the analysis. Specific data for each variable from the training set were utilized to develop six algorithms—CoxPH, RSF, XGBoost, ENET, DeepSurv, and GBM—to predict 5- and 10-year prognoses in the four scenarios.

Fig. 3: Performance of six algorithms in the prediction of medulloblastoma prognosis in CMR scenario.
Fig. 3: Performance of six algorithms in the prediction of medulloblastoma prognosis in CMR scenario.

a Univariate COX regression feature selection for model construction. b Comparison of Area Under the Receiver Operating Characteristic curves (AUROCs) for predicting medulloblastoma prognosis using six algorithms within different time spans in testing set (left) and external validation set (right). c Receiver operating characteristic (ROC) curves with different algorithms for predicting 5- (left) and 10-year (right) survival predictions in medulloblastoma patients in testing set (n = 59). d ROC curves with Extreme Gradient Boosting algorithm (XGBoost) for predicting 5- (left) and 10-year (right) survival predictions in medulloblastoma patients in external validation set (n = 73). e Calibration plots of overall survival predictions for XGBoost algorithm in testing (left) (n = 59) and external validation set (right) (n = 73). f Decision curves for evaluating the clinical utility (net benefit) of XGBoost algorithm for 5- (left) and 10-year (right) survival predictions (n = 73).

Model performance comparison

Given the significant differences between WNT-MB and other molecular subgroups, we included interaction terms in the multivariate Cox analysis. As shown in Supplementary Tables 3 and 4, no significant interactions were found when PFTB boost dose and treatment strategies were analyzed separately with molecular subgroups. Therefore, we combined WNT with other subgroups to construct the predictive model.

In four scenarios (including CMR, CM, CR, and CO), we utilized six algorithms to predict patient prognosis. Their comparative performances are summarized in Supplementary Table 6, demonstrating that the XGBoost and GBM algorithms exhibited moderate yet consistent discrimination and favorable calibration across different scenarios. In the CMR scenario for Gr.3/4-MB in testing set, XGBoost algorithm (IBS = 0.122, C-index = 0.612) exhibited exceptional predictive performance for the prognosis of MB patients, achieving an AUC of 0.601 at 5 years and 0.734 at 10 years (Figs. 3b, 3c), followed by CoxPH algorithm (C-index = 0.578), RSF algorithm (C-index = 0.531), GBM algorithm (C-index = 0.510), ENET algorithm (C-index = 0.495), and DeepSurv algorithm (AUC = 0.465). The calibration plots in testing set presented in Fig. 3e illustrated that XGBoost algorithm maintained good consistency between its predictions and the observations for the 5- and 10-year OS rates (IBS = 0.122). Moving on to the CM scenario, XGBoost algorithm continued to show its superiority. It also boasted the best predictive performance (IBS = 0.131, C-index = 0.609), achieving an AUC of 0.618 at 5 years and 0.737 at 10 years, and the Time-dependent AUC values plot and ROC curves and for ML algorithms are presented in Supplementary Fig. 2a and 2b respectively. The calibration curves in Supplementary Fig. 2d show that the predicted probabilities and observed outcomes for the XGBoost algorithm were similar for the 5- and 10-year overall survival rates (IBS = 0.131). Moreover, in the CR scenario for the testing set, the GBM model demonstrated the best predictive performance in predicting the prognosis of MB patients (IBS = 0.114, C-index = 0.637), with 5-year AUCs of 0.662 and 10-year of 0.736 (Supplementary Fig. 3a, b). As shown in Supplementary Fig. 4a, b, of the abovementioned six algorithms, the GBM model fared best in terms of predicting the prognosis of MB patients in the CO scenario (IBS = 0.112, C-index = 0.635). GBM algorithm also showed good agreement between predicted and observed 5-year and 10-year OS rates in the CR and CO scenarios (Supplementary Figs. 3d, 4d). Overall, the superiority of XGBoost and GBM method was quantitatively validated through C-index, IBS metrics, and calibration plots across testing cohorts, demonstrating their potential for clinical practice.

Predictive performance of disease-free survival (DFS) models

To further address the clinical significance of tumor recurrence, we analyzed DFS, defined as the interval from treatment initiation to the first recurrence or death. In the testing set, the predictive performance across the four scenarios was modest, with C-indices of 0.657 for the CoxPH-based CMR model, 0.632 for the GBM-based CM model, 0.659 for the RSF-based CR model, and 0.643 for the RSF-based CO model. Furthermore, external validation was unfeasible due to the absence of recurrence data in the international cohort. Considering the limited predictive efficacy and the current lack of external generalizability, these DFS models do not yet meet the requirements for robust clinical application. Consequently, we did not develop an interactive web-based calculator for DFS at this stage (Supplementary Fig. 6).

External validation

To evaluate the generalizability of the model, we performed external validation using the XGBoost and GBM algorithms in the international MB cohort. When combining molecular information and radiotherapy strategy in the CMR scenario, as depicted in Fig. 3d, the externally validated ROC curve attained a 5-year AUC of 0.807 (95% CI: 0.685–0.930) and a 10-year AUC of 0.787 (95% CI: 0.610–0.963), which was comparable to that of the testing set (p = 0.12). For the external validation the CM scenario (Supplementary Fig. 2c), the XGBoost algorithm achieved a 5-year AUC of 0.692 (95% CI: 0.584–0.801) and a 10-year AUC of 0.729 (95% CI: 0.581–0.877), which was similar to the AUC obtained in the testing set (p = 0.29). Additionally, for the CR scenario, as illustrated in Supplementary Fig. 3c, the externally validated ROC curve demonstrated a 5-year AUC of 0.722 (95% CI: 0.608–0.836) and a 10-year AUC of 0.727 (95% CI: 0.547–0.906), comparable to that observed in the testing set (p = 0.96). For the CO scenario (Supplementary Fig. 4c), the GBM algorithm achieved a 5-year AUC of 0.730 (95% CI: 0.638–0.822) and a 10-year AUC of 0.712 (95% CI: 0.591–0.833), which was not statistically distinguishable from the AUC obtained in testing set (p = 0.92).

Model calibration is depicted in the external validation set (Fig. 3e, Supplementary Figs. 2d, 3d, and 4d), which show favorable consistency between the predictions and the observed outcomes of the four scenarios. Collectively, the XGBoost and GBM algorithms in the four scenarios showed favorable calibration and consistent performance in external validations.

The DCA further demonstrated the predictive and clinical application potential for XGBoost and GBM algorithms in the four scenarios. In the CMR and CO scenarios, the DCA presented in Fig. 3f and Supplementary Fig. 4e revealed the XGBoost and GBM algorithms each perform optimally across the wide threshold range for predicting the 5- and 10-year OS rates in their corresponding scenarios. Furthermore, regarding the clinical applicability of the CM and CR scenarios, the XGBoost and GBM algorithms achieved a robust net benefit only within a narrow range of threshold probabilities (Supplementary Figs. 2e and 3e).

Model explanation

Given that the SHAP method interprets the final model output by calculating each variable’s contribution to prediction, we employed this method to analyze the results of the XGBoost algorithm. We evaluated the feature-importance rankings based on SHAP values for the CMR scenario(Fig. 4a). In this plot, the contributions of each indicator to the prediction model were assessed using the average SHAP values and presented in descending order as the five most essential features: GFI1 expression level, M stage, Subgroup, MYCN expression level, and MYC expression level. We performed dot plot analysis to uncover the direction and strength of the influence of each feature on model prediction. Features, such as M+ stage and high MYC and GFI1 expression, significantly resulted in the poor prognosis, further underscoring the significance of molecular events in predictive modeling (Fig. 4b). The top 5 most influential variables in the summary plot of the CM scenario were roughly the same as the CMR scenario (Supplementary Fig. 5a). In CR and CO scenarios, histological and molecular subgroups ranked as the top two most important variables (Supplementary Fig. 5b, c). Notably, under these two clinical scenarios, the M stage significantly increases the risk of poor prognosis, while a higher PFTB dose (≥ 55.8 Gy) and RT + CT can dramatically reduce this risk (Supplementary Fig. 5b, c).

Fig. 4: Global model explanation by the SHapley Additive exPlanations (SHAP) method for each feature variable within the final model in testing set (n = 59).
Fig. 4: Global model explanation by the SHapley Additive exPlanations (SHAP) method for each feature variable within the final model in testing set (n = 59).

a Visualization of SHAP summary bar plot depicting the contribution ranking of the XGBoost algorithm’s features in CMR scenario. The bars represent the importance of the variables and their overall contribution to the model prediction. b SHAP summary dot plot. Each patient gets one dot per feature in the model, with the dot color (dark blue = high, light blue = low) showing the actual feature value. Dots stack vertically to show density.

Implementation of the web calculator

The XGBoost-based survival predictor was integrated into a web application for utilization in clinical scenarios. To improve clinical practicality, two interactive web-based Shiny apps by the CMR and CM scenarios were created: https://prognosticmodel.shinyapps.io/Scenario1_CMR/ for cases where the radiotherapeutic dose information was available and https://prognosticmodel.shinyapps.io/Scenario2_CM/ for cases where it was not. The web applications for the CR and CO scenarios are accessible online at the following links: https://prognosticmodel.shinyapps.io/Scenario3_CR/ and https://prognosticmodel.shinyapps.io/Scenario4_CO/. Practical demonstration using a representative case in CMR scenario. By inputting the actual values of the features required for the scenarios, the application can automatically predict the survival rates and clinical risk groups of individual patients with MB. In this case, users input complete data entry by responding to 11 queries and the calculator could predict survival rates at different time spans and the importance of variables (Fig. 5a, b). The results in Fig. 5c showed that LC/AMB, Gr.3, and MYC high level were associated with poorer prognosis, while CSI dose contributed positively to the favorable prognosis.

Fig. 5: Online web application for clinical utility.
Fig. 5: Online web application for clinical utility.

a Practical demonstration using a representative case in CMR scenario. Users input complete data entry by responding to 11 queries. b Survival line chart illustrating predicted survival probabilities at multiple key future time points (12, 36, 60, 80, 100, and 120 months) from the current moment. The y-axis indicates survival probability (%), while the x-axis represents time points. Specific survival values are displayed at the top of the graph. c Feature contribution plot based on SHAP values. Positive SHAP values indicate an increased risk of death, whereas negative values suggest a protective effect.



Source link