We retrieved three publicly available metagenomic sequencing datasets of stool samples from CM patients prior to immunotherapy from the National Center for Biotechnology Information Sequence Read Archive (NCBI-SRA) using the following accession numbers: PRJEB43119 (17), PRJNA399742 (12) and PRJNA915098 (18). The PRJEB43119 dataset includes 165 samples, consisting of 22 complete responses, 42 partial responses, 30 stable diseases, and 71 progressive diseases. In subsequent analyses, we classified complete responses, partial responses, and stable diseases as responders, while progressive diseases were classified as non-responders. The PRJNA399742 dataset includes 14 responders, 23 non-responders, and 1 unknown sample. All samples in the PRJNA915098 dataset are derived from responders. The three cohorts were all of European descent, which helped minimize the potential impact of ethnic heterogeneity on the results. PRJEB43119 will serve as the discovery cohort and training set for machine learning, while PRJNA399742 and PRJNA915098 will be used as the replication cohort and external test sets for machine learning.

To remove low-quality base sequences and adapters, the fastq software (v0.21.0) was used for quality control of the raw sequencing data. Subsequently, the clumpify software (v38.90) will be used for duplicate removal in the data. Based on the human reference genome (GRCh38.p13), Bowtie2 (v2.4.2) software was used to filter out host-derived reads to obtain clean sequence data. Next, MEGAHIT and QUAST were used to assemble the cleaned data and perform gene prediction. CD-HIT software was used to perform redundancy removal on the gene prediction results, resulting in a non-redundant initial gene catalogue. Clustering was performed using default parameters: identity set at 95% and coverage at 90%, with the longest sequence selected as the representative sequence. Finally, Salmon software was used to calculate the gene abundance information for each sample based on the Unigenes gene sequences and the clean data from each sample.

For alpha diversity, the indices assessed include Ace, Chao, Richness, and Shannon. The significance of inter-group differences between the two groups was assessed using the Wilcoxon rank-sum test. For beta diversity, principal components analysis (PCA) was performed based on the OTU table, followed by the creation of a scatter plot. The statistical significance was then assessed using analysis of similarity (ANOSIM). The analyses were performed using the “Vegan” and “ade4” packages in R software.

Biomarker grouping was determined based on the effect sizes from LEfSe analysis (19). The differential analysis was performed using the Wilcoxon rank-sum test, and the results were adjusted for multiple hypothesis correction using the Bonferroni method (Bacterial comparisons: P < 0.0000327 is considered significant, and 0.0000327 < P < 0.05 is considered suggestive evidence. Fungal comparisons: P < 0.000568 is considered significant, and 0.000568 < P < 0.05 is considered suggestive evidence). The selection criteria were an LDA > 2 and a P < 0.05 in this exploratory study. The analysis was performed using the “LEfSe” package in R software.

The diversity and complexity of gut microbiota species level data limit the reliability and robustness of the models. To simplify the data and enhance stability, we opted to construct the RF model at the genus level (20). The biomarkers identified by LEfSe analysis will be used as potential features for the RF model. The “rfcv” function in the “randomForest” package was used to compute the minimum error across different feature subsets using 10-fold cross-validation. Subsequently, to enhance the model’s reliability, we introduced an external test set and used the “pROC” package to plot the receiver operating characteristic (ROC) curve. To evaluate the model’s performance, we employed the area under the curve (AUC), calibration curve, and clinical decision curve analysis (DCA). The accuracy, sensitivity, specificity, and other performance metrics were calculated and reported for test sets.

SHAP values explain the contribution of each feature to the model’s prediction, whether it be positive or negative (21). The feature importance plot is used to display the features that have the greatest impact on the model’s predictions, with feature importance ranked based on the mean absolute SHAP values. The hive plot illustrates the direction of the impact of changes in feature relative abundance on the model’s predictions, which helps us understand the decision-making process within the complex model and further validate the results of LEfSe. This process utilizes the “shapviz” package to interpret the predictions of the RF model.

After identifying the key bacterial and fungal taxa, Spearman’s correlation coefficient was used to assess the key bacterial-fungal interactions in the three cohorts (22). The “stats” package was used to compute the correlations.

All statistical analysis and data visualization were conducted using R version 4.3.0.

The four most abundant bacterial genera in responders and non-responders were Bacteroides, Bifidobacterium, Collinsella, and Fecalibacterium ( Figure 1A ). In the alpha diversity analysis, both the Ace and Richness indices did not show a significant statistical difference (P > 0.05). The Chao index of responders was significantly lower than that of non-responders, suggesting that the diversity of responders may be lower. However, the Shannon index of responders was significantly higher than that of non-responders, indicating a higher evenness in the gut bacterial composition of responders ( Figure 1B ). PCA analysis showed that the first two principal components explained approximately 52.93% of the diversity. Although the confidence intervals of the two groups partially overlap, ANOSIM analysis revealed a significant difference between responders and non-responders (R = 0.0367, P < 0.05) ( Figure 1C ). LEfSe analysis identified 45 potential biomarkers, including 40 bacterial taxa associated with responders and 5 bacterial taxa associated with non-responders. In responders, the top three bacterial taxa identified by LDA analysis were Akkermansia, Collinsella, and Dorea, while in non-responders, the top three were Jonquetella, Xanthovirga, and Roseiflexus ( Figure 1D ; Supplementary Table S1 ).

At the minimum error, the top eight biomarkers ranked by feature importance using the RF algorithm were included in the construction of the diagnostic model ( Figure 2A ). The SHAP summary plot of the model shows the effect of features on the prediction model ( Figure 2B ). The selected features were ranked from highest to lowest based on their average absolute SHAP values. From high to low are: Romboutsia, Endomicrobium, Aggregatilinea, Colwellia, Akkermansia, Candidatus Moduliflexus, Mucispirillum, and Microbacter. The hive plot shows that as the feature value (Relative abundance) of all features increases, their SHAP values tend to predict responders ( Figure 2C ). It is worth noting that these features were identified as responder associated bacterial taxa in the LEfSe analysis, further highlighting their potential value in predicting immunotherapy responses in CM patients.

The test set, composed of two cohorts, was used to evaluate the performance of the RF model. Using the test set, the RF model demonstrated poor performance (AUC = 0.637, Accuracy = 0.578, Sensitivity = 0.536, Specificity = 0.647, F1 Score = 0.584) ( Figure 3A ; Supplementary Table S2 ). We evaluated the accuracy of the RF model in predicting the response of CM patients to immunotherapy by analyzing the calibration curve and DCA. The calibration curve of the test set showed some degree of bias, reflecting the model’s generalization ability to unseen data ( Figure 3B ). The DCA analysis indicated that the model provided limited net benefit for clinical decision-making across most threshold probabilities ( Figure 3C ).

The four most abundant fungal genera in responders and non-responders were Beauveria, Pyricularia, Ceratobasidium and Valsa ( Figure 4A ). None of the four alpha diversity indices showed significant differences between groups ( Figure 4B ). PCA analysis revealed that the first two principal coordinates captured approximately 79.92% of the diversity. The confidence intervals of the two groups almost overlapped. However, the ANOSIM analysis revealed a significant difference between the groups (R = 0.0257, P < 0.05) ( Figure 4C ). The subsequent LEfSe analysis identified four potential gut fungal biomarkers, among which Rasamsonia, Rutstroemia, and Ganoderma were associated with responders, while Zancudomyces was exclusively associated with non-responders ( Figure 4D ; Supplementary Table S3 ).

When the error was minimized, the top three potential fungal biomarkers, ranked by importance based on the RF algorithm, were incorporated into the model construction ( Figure 5A ). The SHAP summary plot of the model illustrates the influence of the features on the prediction ( Figure 5B ). The selected features were ranked from highest to lowest based on their average absolute SHAP values. From high to low are: Ganoderma, Rutstroemia and Zancudomyces. The SHAP hive plot reveals that an increase in the feature value (Relative abundance) of Ganoderma and Rutstroemia is associated with predictions of responders, while an increase in the feature value of Zancudomyces is associated with predictions of non-responders, which is consistent with the results from the LEfSe analysis ( Figure 5C ).

Similarly, a test set composed of data from two cohorts was used to evaluate the model’s predictive performance. Compared to the bacterial based model, the fungal RF model demonstrated better performance. However, its performance was generally suboptimal (AUC = 0.654, Accuracy = 0.489, Sensitivity = 0.393, Specificity = 0.647, F1 Score = 0.489) ( Figure 6A ; Supplementary Table S4 ). The calibration curve of the test set showed some degree of bias, reflecting the model’s generalization ability to unseen data ( Figure 6B ). The DCA analysis indicated that the model provided limited net benefit for clinical decision-making across most threshold probabilities ( Figure 6C ).

Models constructed separately based on bacteria or fungi both exhibited poor predictive performance. We further combined bacteria and fungi to explore the performance of the merged model. Under the condition of minimized error, the top nine potential biomarkers ranked by importance according to the RF algorithm were included in the model construction ( Figure 7A ). Among the biomarkers included in the model, seven bacterial biomarkers were identified, including Candidatus Moduliflexus, Colwellia, Romboutsia, Mucispirillum, Endomicrobium, Aggregatilinea, and Akkermansia. Only two fungal biomarkers were included in the model, namely Rutstroemia and Zancudomyces. The SHAP summary plot of the model illustrates the influence of the features on the prediction ( Figure 7B ). The selected features were ranked from highest to lowest based on their average absolute SHAP values. From high to low are: Romboutsia, Endomicrobium, Aggregatilinea, Zancudomyces, Candidatus Moduliflexus, Colwellia, Akkermansia, Mucispirillum and Rutstroemia. The SHAP hive plot reveals that an increase in the feature value (Relative abundance) of Zancudomyces is associated with the prediction of non-responders, while an increase in other feature values is associated with the prediction of responders, which is consistent with the results from the LEfSe analysis ( Figure 7C ).

After combining the bacterial and fungal data, the performance of the RF model showed a significant improvement compared to the models constructed separately using bacteria or fungi, with the following metrics: AUC = 0.707, Accuracy = 0.644, Sensitivity = 0.679, Specificity = 0.588, and F1 Score = 0.648 ( Figure 8A ; Supplementary Table S5 ). Compared to the previous models, the calibration curve of the merged model shows a better fit, indicating a higher consistency between the model predictions and the actual incidence ( Figure 8B ). DCA indicates that the merged model provides higher net benefit for clinical decisions across most threshold probabilities ( Figure 8C ). Thus, merging the data from gut bacteria and fungi provides a more comprehensive understanding of the gut microbiome and further supports the potential interactions between bacteria and fungi.

Bacteria and fungi involved in the construction of the merged model were identified as important biomarkers. To investigate the interactions between bacteria and fungi, we constructed a microbiome association network by calculating the Spearman correlation coefficient to identified key interactions. In the PRJEB43119, PRJNA399742, and PRJNA915098 datasets, three, two, and three positive correlations were identified, respectively ( Figures 9A–C ). The positive correlation interaction between Akkermansia and Rutstroemia was identified in all three cohorts, which may represent a key interaction associated with immunotherapy response in CM patients. Moreover, the positive correlation interaction between Candidatus Moduliflexus and Rutstroemia was identified in two cohorts, providing relatively strong evidence for it being a key interaction ( Figure 9D ).