Single-cell RNA combined with bulk RNA analysis to explore oxidative stress and energy metabolism factors and found a new prostate cancer oncogene MXRA8

Single-cell expression heterogeneity

Single-cell transcriptome analysis can provide deep insights into the cellular functions within tumor tissues, enabling broad applications in studying tumor cell heterogeneity, immunology, bone marrow transplantation, and disease classification.

In this study, we generated a Seurat object by aggregating expression data from six samples and obtained a total of 22,794 cells from prostate cancer tissue after quality control and filtering. After initial data reduction analysis (PCA), we selected the top 10 PCs for further dimensionality reduction and analysis. We then generated clustering results with resolutions ranging from 0.2–1.5, and created a clustering tree as shown in Figure 1A. When the resolution exceeded 1.1, frequent crossovers between clusters appeared. Based on previous analyses and consensus, we used a resolution of 1.1 for cell subtyping. We used SingleR to annotate cell subtypes (Supplementary Table 1: Meta_info) and corrected the annotations using marker genes from databases to obtain nine cell clusters, including B cells, CD8+ T cells, Mast cells, macrophages, monocytes, epithelial cells, endothelial cells, fibroblasts, and smooth muscle cells, as shown in the UMAP plot in Figure 1B. The distribution of these nine cell types across each patient is shown in Figure 1C, which reveals significant differences in the cell type composition among individuals, indicating tumor cell heterogeneity and complexity. We identified oxidative stress markers and energy metabolism markers for each cell cluster by identifying marker genes (Supplementary Table 1: Marker_Table). The number of oxidative stress markers and energy metabolism markers varied among different cell types, with macrophages having the most markers and immune T cells and B cells having the fewest, as shown in Figure 1D, 1E. Finally, we plotted the expression of oxidative stress and energy metabolism genes that intersected among cell types (OMG genes) as a violin plot, which showed strong heterogeneity in their expression across different cell clusters. The expression of OMG genes was highest in endothelial cells, macrophages, and fibroblasts, as shown in Figure 1F.

Figure 1. Heterogeneity of tumor single cells. (A) Cluster dendrogram of single-cell dataset; (B) UMAP plots of cell clustering and annotation results for single-cell dataset, with outer circle annotations indicating cell populations and inner circles representing patient origins; (C) Cumulative distribution histogram of cell clusters among patients; (D) Table showing the number of oxidative stress marker genes in each cell population; (E) Table showing the number of energy metabolism marker genes in each cell population; (F) Violin plots of OMG marker gene expression in different cell populations.

Inter-cellular differences in oxidative stress

We calculated the ssGSEA (single-sample gene set enrichment analysis) score for each cell’s oxidative stress gene set (Supplementary Table 1: OS_ssGSEA_Score) and compared the differences in oxidative stress function across cell populations using a violin plot. We calculated the ssGSEA score for each cell’s oxidative stress gene set (Supplementary Table 1: OS_ssGSEA_Score) and compared the differences in oxidative stress function between cell populations by plotting a violin plot, as shown in Figure 2A–2C. Macrophages showed strong responses to oxidative stress, oxidative detoxification, and oxidative stress-induced cell death, while CD8+ T cells and B cells had low enrichment scores for all three gene sets. Therefore, we focused on the gene functional features of macrophages. The KEGG GSEA enrichment results of macrophage marker genes (Supplementary Table 1: KEGG_GSEA_Macro) showed that these genes were mainly enriched in pathways related to immunity, such as systemic lupus erythematosus, autoimmune thyroid disease, antigen presentation and processing, and immunoglobulin production, as shown in Figure 2D.

Figure 2. Differences in oxidative stress and communication between subpopulations of cells. (A–C) Violin plots showing the differences in ssGSEA scores for oxidative stress response, oxidative stress-induced cell death, and cellular antioxidant defense among different cell populations; (D) Line plot of KEGG pathway GSEA enrichment results for all marker genes of macrophages.

Expression characteristics of the oxidative stress active population

We obtained 58 intersection genes as oxidative stress-specific markers by taking the intersection of oxidative stress genes and cell subtype-specific marker genes. We then identified the active oxidative stress-related cell population based on the expression of these 58 specific marker genes to study the expression characteristics of oxidative stress-related factors at the single-cell level. The results are shown in Figure 3A: using the optimal threshold (0.26), 2959 cells were identified as the active oxidative stress-related population (Supplementary Table 2: AUC_res). They were mainly distributed among monocytes, macrophages, endothelial cells, and some fibroblasts in the single-cell UMAP plot, as shown in Figure 3B. We then plotted the distribution of active and inactive oxidative stress cells in each cell cluster as a column chart, as shown in Figure 3C. To explore the functional differences between the active and inactive cell populations, we identified 712 differentially expressed genes between the active and inactive oxidative stress populations through differential expression analysis. We then performed KEGG and GO_BP gene set enrichment analyses based on these differentially expressed genes, as shown in Figure 3D, 3E. The enriched pathways identified by KEGG included rheumatoid arthritis, COVID-19, phagosomes, antigen processing and presentation, viral myocarditis, transplant rejection, Staphylococcus aureus infection, human T-cell leukemia virus 1 infection, and type 1 diabetes. The enriched functions identified by GO_BP included leukocyte-cell adhesion and regulation, positive regulation of leukocyte activation, leukocyte migration, positive regulation of lymphocyte activation, and biological processes involving MHC class II antigen processing and presentation of peptides or polysaccharide antigens.

Figure 3. Identification and functional characterization of oxidative stress-activated populations. (A) AUC scores for specific oxidative stress marker genes; (B) UMAP color-coded plot based on cell activity scores, where brighter colors indicate higher response scores and activity levels to the gene set; (C) Histogram showing the distribution of active and inactive oxidative stress cell populations in different cell groups; (D) Bubble plot of KEGG enrichment results for differentially expressed genes, with bubble size representing the number of enriched genes and color indicating the significance of the enrichment result; (E) Bubble plot of GOBP enrichment results for differentially expressed genes.

Expression characteristics of the active energy metabolism population

Similarly, we identified the active energy metabolism population based on the expression of 19 specific energy metabolism marker genes. The results are shown in Figure 4A: using the optimal threshold (0.18), 2969 cells were identified as the active energy metabolism population (Supplementary Table 2: AUC_res). They were mainly distributed among fibroblasts and some smooth muscle cells in the single-cell UMAP plot, as shown in Figure 4B. We then plotted the distribution of active and inactive oxidative stress cells in each cell cluster as a column chart, as shown in Figure 4C. Similarly, differential expression analysis identified 656 differentially expressed genes between the active and inactive energy metabolism populations, and KEGG and GO_BP gene set enrichment analyses were performed based on these differentially expressed genes, as shown in Figure 4D, 4E. The enriched pathways identified by KEGG included COVID-19, focal adhesion, ECM-receptor interaction, leukocyte transendothelial migration, protein digestion and absorption, malaria, ribosome, complement and coagulation cascades, and rheumatoid arthritis, as well as the hematopoietic cell lineage. The enriched gene sets identified by GO_BP included cell structure-related processes, wound healing, collagen metabolic processes, regulation of peptidase activity, collagen fibril organization, negative regulation of peptidase activity, cell-matrix adhesion, and negative regulation of endopeptidase activity.

Figure 4. Identification and characterization of energy metabolism-active populations. (A) AUC scores for specific energy metabolism marker genes; (B) UMAP color-coded plot based on cell activity scores, where brighter colors indicate higher response scores and activity levels to the gene set; (C) Histogram showing the distribution of active and inactive energy metabolism cell populations in different cell groups; (D) Bubble plot of KEGG enrichment results for differentially expressed genes; (E) Bubble plot of GOBP enrichment results for differentially expressed genes.

Differential expression between tumor and normal samples in the TCGA cohort

To begin, we conducted an analysis of differential gene expression between normal and tumor tissue samples within the TCGA PRAD cohort (Supplementary Table 3: deg_res) and obtained 1974 differentially expressed genes, including 762 upregulated genes and 1212 downregulated genes in tumors. Among all the differentially expressed genes, there were 88 oxidative stress-related genes, and their differential expression heatmap is shown in Figure 5A: among them, 28 oxidative stress-related genes were upregulated in tumors, while 60 were downregulated in tumor tissue. Similarly, among all the differentially expressed genes, there were 81 energy metabolism-related genes, and their differential expression heatmap is shown in Figure 5B: among them, 24 energy metabolism-related genes were upregulated in tumors, while 57 were downregulated in tumor tissue. In addition, there were 7 OMG differentially expressed genes, and their expression patterns are shown in Figure 5C: they were all downregulated in tumor tissue. Finally, we used ssGSEA scores to evaluate the enrichment differences of oxidative stress-related pathways between tumor and normal tissues, as shown in Figure 5D–5F: the oxidative stress response, oxidative detoxification, and oxidative stress-induced cell death in prostate cancer tissue were significantly weaker than those in normal tissue.

Figure 5. Dysregulation of oxidative stress and energy metabolism-related genes. (A) Heatmap of significantly differentially expressed oxidative stress genes, with red indicating high expression and blue indicating low expression, and orange and green indicating tumor and normal tissue samples, respectively; (B) Heatmap of significantly differentially expressed energy metabolism-related genes; (C) Box plot showing the expression of commonly differentially expressed oxidative stress and energy metabolism-related genes, with orange and green representing tumor and normal samples, respectively; (D–F) Violin plots showing ssGSEA score differences for oxidative stress-related gene sets in tumor and normal tissues, including oxidative stress response, cellular antioxidant defense, and oxidative stress-induced cell death.

Identification of potential patients based on differential OM genes

Next, we performed correlation analysis between the differentially expressed oxidative stress-related genes and energy metabolism-related genes, and selected genes that were highly correlated with at least one gene from the other functional group (|correlation|>0.4 and p < 0.001). Finally, we obtained a set of 63 genes with differential associations between oxidative stress and energy metabolism, and used them for clustering analysis to identify tumor subtypes. We then plotted Kaplan-Meier survival curves for each cluster subtype. The best clustering result was achieved when using the KM algorithm and Euclidean distance, with a best K value of 2 (Figure 6A–6D), and there was a significant difference in survival curves between the two subtypes (Figure 6E–6H). The sample clustering results are shown in Supplementary Table 4. The consistency cluster cumulative distribution function (CDF) is shown in Figure 6I, which displays the CDF for different values of k. Figure 6J shows the change in the area under the CDF curve ask increases relative to k-1. We then used PCA to reduce the dimensionality of the data and plotted a two-dimensional scatter plot to validate the clustering results, as shown in Figure 6K: the two subtypes were clearly separated, indicating good clustering results. Ultimately, we identified two significant oxidative stress and energy metabolism molecular subtypes with significantly different prognoses, where cluster 1 (C1) had a significantly better prognosis than cluster 2 (C2).

Figure 6. Identification of prostate cancer subtypes. (A–D) Clustering results for k = 2, k = 3, k = 4, and k = 5; (E–H) Survival curves for k = 2, k = 3, k = 4, and k = 5, with different colors representing different clusters; (I) CDF curve distribution for consensus clustering; (J) Distribution of area under the CDF curve for consensus clustering; (K) PCA scatter plot showing the classification results, with red representing C1 and blue representing C2.

Differential expression between tumor subtypes

Based on the clustering results of tumor samples, we performed differential expression analysis between the two subtypes and identified 1208 genes with differential expression (Supplementary Table 4: cc_deg). Among them, there were 64 differentially expressed oxidative stress-related genes, and their volcano plot is shown in Figure 7A: only one oxidative stress gene was upregulated in C2, while the rest were upregulated in C1. There were also 52 differentially expressed energy metabolism-related genes, and their volcano plot is shown in Figure 7B: most genes were upregulated in C1. To identify functional differences between prostate cancer molecular subtypes, we performed KEGG and GOBP pathway enrichment analysis based on the differentially expressed genes (Supplementary Table 4) and plotted a bubble chart of the top 10 significantly enriched pathways, as shown in Figure 7C–7D: the differentially expressed genes mainly enriched in functional pathways such as cell adhesion molecules, ECM-receptor interaction, Staphylococcus aureus infection, malaria, amoebiasis, protein digestion and absorption, hematopoietic cell lineage, PI3K-Akt signaling pathway, viral myocarditis, and biological processes related to cell structure, cell adhesion regulation, leukocyte migration, and wound healing. We then compared the enrichment differences of six immune-related biological factors using the ssGSEA algorithm, as shown in Figure 7E: C1 had higher enrichment scores than C2. Similarly, we calculated the enrichment of four immune function gene sets, as shown in Figure 7F: C1 had higher immune activity than C2.

Figure 7. Differential expression and characterization of molecular subtypes in prostate cancer. (A, B) Volcano plots of differentially expressed genes between subtypes, with red indicating upregulation, blue indicating downregulation, and gray indicating no significant difference in gene expression. The labeled genes in the two plots are oxidative stress genes and energy metabolism genes, respectively; (C, D) Bubble plots of KEGG pathway and GOBP functional enrichment analysis results for differentially expressed genes between subtypes; (E, F) Box plots showing differences in enrichment scores for immune-related factor gene sets and immune function gene sets between subtypes.

Differences in immune infiltration between prostate cancer molecular subtypes

The composition of the tumor microenvironment includes stromal and immune cells, which have been proposed as valuable components for diagnosing and prognosticating tumors. In this study, we estimated the proportions of immune infiltrating cells within each tumor sample using four common methods (Supplementary Table 5), and displayed differences in immune infiltration between two subtypes of prostate cancer via box plots. Supplementary Figure 1A depicts the differences in the infiltration proportions of 28 immune cell types calculated by ssGSEA, with 25 cell types exhibiting significant differences in infiltration proportions. Supplementary Figure 1B shows box plots indicating differences in the infiltration proportions of six immune cell types calculated by the TIMER algorithm, where all six immune cell types exhibited significantly higher infiltration proportions in C1 compared to C2. Box plots shown in Supplementary Figure 1C display differences in scores for four metrics calculated using ESTIMATE, with the stromal score, immune score, and ESTIMATE score of C1 being significantly higher than those of C2, while tumor purity was lower in C1 compared to C2. The results for immune infiltrating cell proportions calculated by CIBERSORT and xCell are presented in Supplementary Figures 2 and 3.

Identification of prostate cancer prognostic signature

In the process of cancer development, oxidative stress and energy metabolism interact with each other and influence the evolution of tumor cells together. In order to study the impact of both on tumor prognosis, we first identified the intersection of differentially expressed genes identified by three differential analyses (differential expression between TCGA molecular subtypes, differential expression of oxidative stress active cell populations, and differential expression of energy metabolism active cell populations) using a Venn diagram, as shown in Figure 8A: we obtained 185 common DEGs. Then, we analyzed these common differentially expressed genes using univariate Cox analysis and finally identified 20 genes associated with prostate cancer prognosis (Supplementary Table 6: cox_res), as shown in the forest plot of univariate Cox results in Figure 8B. Subsequently, we randomly sampled 2/3 of the TCGA_PRAD overall set (n = 478) as the training set (n = 320), used Lasso linear regression to remove redundant genes, and set seed = 1110. We ultimately selected 9 prognostic-related signatures (Supplementary Table 6) through this process, as shown in Figure 8C–8E. We utilized the median gene expression value to stratify high and low expression groups, and subsequently constructed Kaplan-Meier survival curves for the selected model genes within the overall TCGA cohort. Remarkably, all model genes exhibited significant differences in KM curves between high and low expression groups as illustrated in Figure 8F.

Figure 8. Lasso regression analysis results for TCGA dataset. (A) Venn diagram of differentially expressed genes between subtypes and active cell populations; (B) Forest plot of single-factor Cox analysis results for common differentially expressed genes; (C) Trajectory of changes in Lasso regression independent variables, with the x-axis representing the logarithm of the independent variable Lambda and the y-axis representing the coefficient of the independent variable; (D) Confidence intervals for each Lambda in Lasso regression; (E) Lasso regression coefficients for nine key prognostic factors; (F) KM curves for model genes, with red indicating high expression group, blue indicating low expression group, and p-value indicating the significance of survival curve differences.

Validation of model performance using TCGA dataset and external dataset

To further evaluate the impact of the 9 signatures on the recurrence of the training set, we first calculated the RiskScore for each sample in the dataset using the formula RiskScore = MYLK × −0.144 + FBLN1 × −0.059 + PTGDS × −0.047 + TGFB1 × 0.372 + MXRA8 × 0.151 + PTN × −0.056 + COL1A1 × 0.034 + THBS2 × 0.082 + MT1E × −0.059 (Supplementary Table 6: TCGA_Train), the median value of RiskScore (0.4329028) was used as a critical threshold to categorize samples into high-risk and low-risk groups. The gene expression heatmap of the model genes is presented in Figure 9A, showing their distribution in the high-risk and low-risk groups of the training set. To observe the relationship between survival and score, scatter plots of survival time, survival status, and sample risk scores were generated (Figure 9B, 9C). Kaplan-Meier curves were also plotted to compare the prognostic differences between the two groups, revealing that the high-risk group had significantly lower survival rates compared to the low-risk group (Figure 9D; p-value < 0.001). Moreover, ROC curves were constructed based on the risk model, demonstrating its good predictive capability with AUC values of 0.755/0.731/0.699 for 1/3/5 years (Figure 9E).

Figure 9. Validation of prognostic performance of the model using TCGA training set. (A–C) Risk triad plots of TCGA training set, including heatmap of model gene expression in risk score groups, scatter plot of risk scores, and scatter plot of survival time, with red representing high-risk group and blue representing low-risk group; (D) KM curve of TCGA training set; (E) ROC curve of TCGA training set.

To further validate the reliability of the model score in predicting overall survival of prostate cancer patients, we conducted similar analyses on the TCGA PRAD test and overall sets, as well as external datasets GSE70770 and GSE70769, and obtained similar conclusions. These findings suggest the robustness of the model score in predicting overall survival of prostate cancer patients across different datasets (Supplementary Figures 4–7).

Independence of risk score

To verify whether RiskScore can serve as an independent prognostic factor, this analysis combined prostate cancer Tumor_type, Age, T_Stage, N (N stage), M (M stage), MSI, and Race information (Supplementary Table 6: Clinical_stat) for single and multiple factor Cox regression analysis. The results are shown in Figure 10A: in the single-factor Cox regression, Riskscore, T_Stage, N, MSI, and Race was significantly different from the Reference. In the multiple-factor Cox analysis, Riskscore and T_Stage was significantly different from the Reference in the presence of other clinical features, proving that they are independent prognostic factors. In addition, a nomogram (Figure 10B) was drawn based on survival time, survival status, Age, and Stage, showing that Riskscore contributed the most as a clinical factor. Then, cumulative distribution plots of Age, T_Stage, and N grouping samples between high-risk and low-risk groups were drawn and Fisher’s exact test was performed (Figure 10C–10E): all p-values were less than 0.05, indicating a significant correlation between their distribution and the risk model grouping. We then studied the differences in Riskscore distribution between clinical feature groups, as shown in Figure 10F–10J: Riskscore was significantly different among Age, T_Stage, N stage, M stage, and prostate cancer subtype, indicating that multiple clinical features of cancer are related to the risk score.

Figure 10. Independence of risk score with tumor clinical characteristics. (A) Forest plot of single and multiple Cox analysis results for Riskscore and clinical factors in TCGA cohort; (B) Column chart of the prediction model, with the square plus line segment representing the contribution of each clinical factor to outcome events, Total Points representing the total score obtained by adding up the scores of all variables, and the three lines at the bottom representing the cumulative survival probability at 5/3/1 years for each value point. (C–E) Cumulative distribution histogram of different clinical feature groups between high and low-risk groups, with p-value indicating the significance of differences; (F–J) Box plot showing differences in risk score distribution among different clinical feature groups, with different colors representing different groups and p-value indicating the significance of differences.

Differences in immune microenvironment between model groups

To determine variations in the tumor immune microenvironment between high-risk and low-risk groups, we initially computed associations between Riskscore, the expression of nine model genes, and 23 immune checkpoint inhibitors (Supplementary Table 7). Next, we displayed these correlations as a heatmap (Figure 11A) for visualization purposes. Results demonstrated significant positive correlations between model gene expression and immune checkpoint inhibitor expression for CSF1R, CTLA4, HAVCR2, and TGFBR1. We then compared expression differences for 23 immune checkpoint inhibitors between high-risk and low-risk groups, presenting box plots depicting results in Figure 11B: 16 immune checkpoint inhibitors demonstrated differential expression, with higher expression levels observed in the high-risk group. Finally, based on analysis of immune infiltrating cell proportions, we compared differences in immune cell infiltration between high-risk and low-risk groups. Box plots in Figure 11C illustrate a comparison of four ESTIMATE scores between high-risk and low-risk groups. Our results indicated that the high-risk group exhibited higher stromal score, immune score, and ESTIMATE score, but lower tumor purity, compared to the low-risk group, with significant differences observed. Figure 11D depicts differences in the infiltration proportions of 28 immune cell types calculated by ssGSEA between high-risk and low-risk groups, with infiltration proportions of 22 cell types exhibiting significant differences. Additional algorithm results are presented in the supplementary material.

Figure 11. Differences in immune microenvironment among model groups. (A) Correlation coefficient heatmap between model and immune checkpoint expression, with point color representing the level of correlation and ^*indicating significance; (B) Box plot showing differential expression of 23 immune checkpoint inhibitors between high and low-risk groups; (C) Box plot showing differences in stromal score, immune score, ESTIMATE score, and tumor purity, with red and blue representing high and low-risk groups, respectively; (D) Box plot showing differences in immune cell infiltration proportions calculated by ssGSEA algorithm between high and low-risk groups, with red and blue representing high and low-risk groups, respectively.

Enrichment differences

Based on the HALLMARK pathway enrichment scores of prostate cancer samples (Supplementary Table 8), combined with model grouping information, we explored the pathway enrichment differences between high-risk and low-risk groups, which can help study the relationship between cancer characteristic pathways and prognosis. Supplementary Figure 8A shows that multiple pathways have significantly different enrichment scores between model groups, including fatty acid metabolism, interferon α/γ response, and inflammatory response. We then calculated the correlation between Riskscore and model gene expression and HALLMARK pathway enrichment scores, as shown in Supplementary Figure 8B: RiskScore was significantly negatively correlated with pathways such as fatty acid synthesis, bile acid metabolism, cholesterol homeostasis, and fatty acid metabolism, and significantly positively correlated with pathways such as angiogenesis and G2M checkpoint.

Correlation between risk model and prostate cancer genome mutations

Genomic mutations play an important role in cancer progression and therapy development. To investigate the distribution of somatic variations between high-risk and low-risk groups and gene mutations among samples with different clinical features, we selected the top 30 genes with the highest mutation frequency in both groups and created waterfall plots (Figure 12A, 12B). PRAD featured TP53, TTN, and SPOP as top three genes with the highest mutation frequency, but their mutation frequencies varied between high- and low-risk groups. Comparing mutation frequencies between the two groups revealed significant differences for TP53, FOXA1, RYR1, and PTEN, which are genes with high mutation rates, as shown in cumulative distribution bar charts (Figure 12C–12F). TP53 mutation was further associated with prostate cancer recurrence, resulting in worse prognosis (Figure 12G). Finally, we compared tumor mutation burden (TMB) between high- and low-risk groups, revealing higher TMB in the high-risk group (Figure 12H). By dividing samples into high-TMB and low-TMB groups based on median TMB value, we plotted KM survival curves for both groups (Figure 12I), demonstrating significantly worse prognosis and survival differences for patients with high TMB.

Figure 12. Genomic mutation differences among model groups. (A) Waterfall plot of the top 30 genes (mutation frequency) with SNV in high-risk group; (B) Waterfall plot of the top 30 genes (mutation frequency) with SNV in low-risk group; (C–F) Cumulative distribution histogram of frequently mutated genes with significant differences between high and low-risk groups; (G) Survival curve showing the difference between mutant and wild-type groups of gene TP53; (H) Box plot showing the difference in TMB between high and low-risk groups; (I) KM survival curve showing the difference between high and low TMB groups, with red representing high TMB group and blue representing low TMB group.

Evaluating chemotherapy drug resistance using the risk model

Using the expression data of TCGA PRAD, we predicted the IC₅₀ values of 138 drugs from the GDSC database (Supplementary Table 9). We then performed statistical tests to examine whether there were significant differences in IC₅₀ values between high-risk and low-risk groups based on the model grouping results and IC₅₀ values. We selected drugs with significantly different IC₅₀ values and calculated their correlation with Riskscore, finding that five drugs showed a highly positive correlation with Riskscore (Figure 13A–13E). Next, we plotted box plots displaying the distribution differences of IC₅₀ values for these five drugs between high-risk and low-risk groups (Figure 13F–13J). These results indicated that the low-risk group exhibited greater sensitivity to these five chemotherapy drugs, as evidenced by significantly lower IC₅₀ values compared to the high-risk group.

Figure 13. Correlation between risk score and chemotherapy drug resistance. (A–E) Scatter plots showing the correlation between risk score and drug IC₅₀ values, with R representing the correlation coefficient and p-value indicating the significance of the correlation; (F–J) Box plots showing differences in drug IC₅₀ values between high and low-risk groups.

Predicting the efficacy of immunotherapy using the risk model

We started by calculating the risk scores for each sample in the immunotherapy dataset and dividing them into high- and low-risk groups based on median (Supplementary Table 10: IMvigor210_res). Significant survival differences were observed between the two groups using KM curves (Figure 14A). Next, we compared immune response differences between high- and low-risk groups in the immunotherapy cohort, demonstrating that responsive samples were significantly more frequent in the low-risk group than in the high-risk group (Figure 14B). Further analysis revealed a significant difference in risk score between responders and non-responders depicted in Figure 14C, suggesting that the low-risk group is more likely to benefit from immunotherapy. Additionally, we used TIDE online analysis to predict immune response of samples in the TCGA_PRAD dataset and evaluated differences in immune response between oxidative stress and metabolic subtypes (Figure 14D), showing that the proportion of responders was significantly higher in C2 subtype compared to high-risk group and C1 subtype. We also assessed the IPS score to determine the tumor immunogenicity and predict response to immunotherapy. Analysis of four categories of IPS scores between the high-risk and low-risk groups (ips_ctla4_neg_pd1_neg, ips_ctla4_pos_pd1_neg, ips_ctla4_neg_pd1_pos, ips_ctla4_pos_pd1_pos) showed that two IPS scores were significantly higher in the low-risk group, while the other two did not exhibit significant differences between the groups, indicating that patients in the low-risk group are more likely to benefit from immunotherapy (Figure 14E, 14F).

Figure 14. Prediction of patient response to immunotherapy using the risk model. (A) KM curves of high and low-risk groups in the IMvigor210 cohort; (B) Cumulative distribution histogram of response and non-response groups between high and low-risk groups in the IMvigor210 cohort; (C) Violin plot showing differences in risk score between response and non-response groups in the IMvigor210 cohort; (D) Cumulative distribution histogram of TIDE-predicted immune response between response and non-response groups in the clustering subtype of the IMvigor210 cohort; (E, F) Violin plots showing differences in IPS score between high and low-risk groups, with red and blue representing high and low-risk groups, respectively.