This study included two cohorts comprising 151 dMMR CRC patients who received neoadjuvant anti-PD-1/PD-L1 immunotherapy. The training cohort (SYSU-SAH, n=121) was derived from the Sixth Affiliated Hospital of Sun Yat-sen University (SYSU-SAH), with serum samples collected between January 2019 and January 2024 for model development. The external validation cohort (SYSUCC, n=30) was obtained from the Cancer Center of Sun Yat-sen University (SYSUCC), with samples collected from January 2023 to January 2024 for model performance evaluation. This study was approved by the ethics committees of both participating institutions: the Sixth Affiliated Hospital of Sun Yat-sen University (Approval No. 2025ZSLYEC-397) and Sun Yat-sen University Cancer Center (Approval No. B2023-188-04). Immunotherapy efficacy was evaluated by postoperative pathological assessment. Patients were classified as responders or non-responders based on major pathological response (MPR, defined as ≤10% residual viable tumor cells) or pathological complete response (pCR, defined as no residual viable tumor cells). Clinical data, including sex, age, tumor stage, smoking history, tumor location, tumor differentiation, and mismatch repair protein expression, were obtained from medical records. Missing clinical data were addressed using multiple imputation methods.

Serum samples from 15 immunotherapy responders and 15 non-responders in the training cohort were processed for targeted metabolomics analysis to identify potential metabolic biomarkers. Differential metabolites were screened based on the following criteria: Fold change (FC) > 1.5 or < 0.67 and False discovery rate (FDR) < 0.05, which identified nine candidate metabolites associated with immunotherapy response. Subsequently, the concentrations of these nine metabolites were quantified in serum samples. To avoid data leakage and ensure rigorous validation, feature selection was performed exclusively using the training cohort (SYSU-SAH, n=121). Five of these differential metabolites demonstrated significantly statistical difference in the training cohort and were selected as features for model development.

Two sets of predictive models were developed. The first set utilized the concentrations of five differentially expressed metabolites (PGE₂, tryptophan, arginine, citrulline, and histidine) as input features. The second set combined these five metabolites with clinical variables, including sex, age, tumor stage, smoking history, tumor location, tumor differentiation, and mismatch repair protein expression, as input features. To assess for multicollinearity among the five metabolite features, we calculated the Variance Inflation Factor (VIF). The VIF values for PGE2 (1.13), Tryptophan (1.08), Histidine (1.08), Arginine (1.16), and Citrulline (1.04) were all well below the common threshold of 5, indicating that multicollinearity was not a significant concern for our model development. Both sets of models were constructed using six machine learning algorithms: logistic regression (LR), k-nearest neighbors (KNN), random forest (RF), support vector machine (SVM), linear discriminant analysis (LDA), and extreme gradient boosting (XGBoost). The models were trained in the training cohort (n = 121) with 5-fold cross-validation (repeated five times) to optimize performance and mitigate overfitting. Hyperparameters for each algorithm were optimized using grid search to maximize the area under the receiver operating characteristic curve (ROC-AUC). Feature selection was conducted using recursive feature elimination (RFE) combined with the random forest algorithm to identify the most predictive features. The SHAP (SHapley Additive exPlanations) method was employed to evaluate feature importance, with SHAP values calculated for each feature based on the RF model. All models were implemented in the R environment using the “tidymodels” and “SHAP” (version 0.42.0) packages. Data preprocessing involved standard deviation normalization (StandardScaler) and multiple imputation for missing values (missing rate < 5%).

Model performance was assessed using multiple metrics, including ROC-AUC, accuracy, F1 score, Matthews correlation coefficient (MCC), negative predictive value (NPV), positive predictive value (PPV), recall, sensitivity, and specificity. These metrics were calculated separately in the training cohort (n = 121) using 5-fold cross-validation and in the external validation cohort (n = 30) to evaluate robustness and generalizability of the model. The predictive performance of each model was compared to identify the optimal algorithm for predicting neoadjuvant immunotherapy efficacy.

Statistical analyses were conducted using R (version 4.5.0). Continuous variables were compared using Student’s t-tests or Wilcoxon rank-sum tests, depending on data distribution. Categorical variables were analyzed with chi-square tests. The statistical comparison of the Receiver Operating Characteristic (ROC) curves between the 5-MPM model (containing only five metabolites) and the combined model (integrating the five metabolites with M stage) was performed using DeLong’s test. A p-value < 0.05 was considered statistically significant. Machine learning analyses were performed using the tidymodels package, with standardized data preprocessing steps, including normalization via StandardScaler and multiple imputation for missing values (missing rate < 5%).

Serum samples from 151 dMMR/MSI-H CRC patients with neoadjuvant anti-PD1/PD-L1 immunotherapy were analyzed, their clinical characteristics were provided in Table 1. In SYSUSAH cohort (n = 121), 95 patients were classified as responders and 26 as non-responders, whereas SYSUCC cohort (n = 30) consisted of 22 responders and 8 non-responders. Serum samples from 15 immunotherapy responders and 15 immunotherapy non-responders in SYSUSAH cohort were processed for targeted metabolomic, which revealed a panel of candidate metabolites distinguishing responders from non-responders. These metabolites were subsequently validated in the entire study population, and multiple predictive models were developed based on the validated biomarkers. The optimal model demonstrated robust performance and generalizability. The overall workflow of patient enrollment, metabolite screening, model construction, validation, and interpretation was presented in Figure 1.

To identify dysregulated metabolites associated with immunotherapy response in dMMR/MSI-H CRC patients, serum samples from 15 non-responders (NR) and 15 responders (R) in SYSUSAH cohort were processed for targeted metabolomic analysis. Principal component analysis (PCA) and orthogonal projections to latent structures–discriminant analysis (OPLS-DA) revealed distinct serum metabolic profiles between responders and non-responders (Figures 2A, B). Compared with responders, nine differential metabolites were identified in serum from the non-responders (FDR < 0.05; Fold change > 1.5 or < 0.67) (Figure 2C). These dysregulated metabolites included seven amino acids (threonine, isoleucine, histidine, tryptophan, citrulline, phenylalanine, arginine), an inflammatory mediator prostaglandin E2 (PGE2), and a methylation-related metabolite sarcosine (Figure 2C). Unsupervised clustering analysis showed that these nine metabolites could clearly classified the 30 samples into NR and R groups (Figure 2D). These differential metabolites levels were subsequently validated with enzyme-linked immunosorbent assay (ELISA) in training Cohort (SYSU-SAH, n = 121). This validation results revealed that five metabolites (PGE2, arginine, citrulline, histidine, and tryptophan) were consistently significantly different between R and NR groups (P < 0.05, Figure 2E).

Based on above five differential metabolites (PGE2, arginine, citrulline, histidine, and tryptophan), predictive models for immunotherapy efficacy were developed with six machine learning algorithms applied to the SYSUSAH (training) cohort (n = 121). Upon evaluation with five-fold cross-validation repeated five times, the RF model demonstrated the best overall performance, achieving an AUC of 0.85 (95% CI: 0.81–0.88) and an MCC of 0.50 (95% CI: 0.42–0.58), thereby demonstrating robust stability in handling class-imbalanced data. In an independent testing cohort (n = 30), the RF model exhibited strong generalization capability, as evidenced by 1000 bootstrap resampling validations: AUC increased to 0.88 (95% CI: 0.87–0.89) and MCC reached 0.76 (95% CI: 0.75–0.77) (Figures 3A, B). Although the Support Vector Machine model achieved a slightly higher AUC (reported as 0.89; 95% CI: 0.89–0.90), its performance on other metrics—including precision (0.77, 95% CI: 0.76–0.77), MCC (0.48, 95% CI: 0.47–0.49), sensitivity (0.75, 95% CI: 0.74–0.76), specificity (0.77, 95% CI: 0.77–0.78), and F1-score (0.62, 95% CI: 0.61–0.63)—was inferior to the RF model (Figure 3C), indicating that RF provided more balanced and comprehensive predictive performance. In summary, the metabolite-based 5-MPM model developed with Random Forest demonstrated excellent and well-balanced performance in both cross-validation and independent testing, underscoring its potential application in predicting immunotherapy efficacy.

To further enhance predictive performance, we combined above five metabolite biomarkers (PGE2, arginine, citrulline, histidine, and tryptophan) with M stage, which was associated with immunotherapy efficacy (p < 0.05) in our cohort. Using the training set from SYSUSAH cohort (n = 121), a combined predictive model was constructed and evaluated with five-fold cross-validation repeated five times. The RF integrated model demonstrated superior overall performance during cross-validation, with an accuracy of 0.849 (95% CI: 0.830–0.869) and AUC of 0.88 (95% CI: 0.85–0.91) (Figures 4A, C). For external validation, we tested this model in another independent dMMR/MSI-H cohort from external validation cohort (SYSUCC, n = 30) using 1000 bootstrap resamples. The RF integrated model achieved an accuracy of 0.87 (95% CI: 0.86–0.87) and AUC of 0.88 (95% CI: 0.88–0.88) (Figures 4B, C). Overall, although the RF-integrated model (five metabolites and M stage) achieved the best performance, its improvement over the metabolites-based 5-MPM model was not statistically significant (P > 0.05).

To interpret the metabolites-based 5-MPM model, SHAP method was applied to quantify the contribution of each metabolite to its predictive performance. As shown in the SHAP summary bar plot (Figure 5A), mean absolute SHAP values ranked metabolites by their importance in 5-MPM model in descending order: PGE₂, histidine, tryptophan, arginine, and citrulline. Furthermore, the SHAP summary bee swarm plot depicts the direction and magnitude of each metabolite’s impact on the model output. Elevated levels of PGE₂ and citrulline were associated with higher SHAP values, whereas lower levels of histidine, tryptophan, and arginine contributed to better prediction scores (Figure 5B). Moreover, the SHAP dependence plots further elucidated how individual metabolites influenced the model output. As shown in Figure 5C, a SHAP value greater than zero indicated a positive contribution toward predicting non-response (NR), corresponding to an increased risk of poor immunotherapy outcome. Thus, higher levels of PGE₂ and citrulline, together with lower levels of tryptophan, arginine, and histidine, were associated with increased SHAP values, suggesting a greater likelihood of inferior immunotherapy efficacy.

To verify the necessity of the selected metabolite panel, we performed a “leave-one-out” ablation study. Specifically, we trained five derivative models, each excluding exactly one metabolite, while maintaining identical hyperparameter settings and training/validation splits as the original 5-MPM model. Our analysis revealed that the full 5-MPM model (Accuracy = 0.90, AUC = 0.88) consistently outperformed all reduced models. As shown in Supplementary Table S1, the exclusion of any single metabolite resulted in a decrease of predictive performance (Accuracy range: 0.63–0.87; AUC range: 0.68–0.88). These findings confirm that each metabolite provides unique and essential information to the 5-MPM model.

To evaluate the clinical utility of the 5-MPM model, the decision curve analysis (DCA) was performed across the training and external validation cohorts. As illustrated in Figure 6, the DCA curves demonstrate that the 5-MPM model yielded a higher net benefit for predicting immunotherapy response across a broad range of threshold probabilities compared with the ‘treat-all’ or ‘treat-none’ strategies. These findings suggest that the 5-MPM model could serve as a valuable tool for guiding clinical decision-making in neoadjuvant immunotherapy for dMMR/MSI-H CRC patients.