Integration of Graph Neural Networks and multi-omics analysis identify the predictive factor and key gene for immunotherapy response and prognosis of bladder cancer

Abstract

Objective

The evaluation of the efficacy of immunotherapy is of great value for the clinical treatment of bladder cancer. Graph Neural Networks (GNNs), pathway analysis and multi-omics analysis have shown great potential in the field of cancer diagnosis and treatment.

Methods

A GNNs model was constructed to predict the immunotherapy response and identify key pathways. Based on the genes of key pathways, bioinformatic methods were used to generate a simple linear scoring model, namely responseScore. The intrinsic mechanism of responseScore was explored from the perspectives of multi-omics analysis. The relationship between each gene involved in responseScore and prognosis was also explored. Transfection experiments with human bladder cancer cells were used to investigate the biological effects of PSMB9 gene.

Results

The final GNNs model had an AUC of 0.785 on the training set and an AUC of 0.839 on the validation set. R-HSA-69620 and others were identified as key pathways. ResponseScore had a good performance in predicted the immunotherapy response and prognosis. Analysis results from genetic variation, pathways and tumor microenvironment, showed that responseScore was significantly associated with immune cell infiltration and anti-tumor immunity. The results of single-cell analysis showed that responseScore was closely related to the functional state of natural killer cells. Compared with the PCDH-NC group, cell migration and proliferation were significantly inhibited while cell apoptosis increased in the PCDH-PSMB9 group.

Conclusion

The GNNs predictive model and responseScore constructed in this study can reflect the immunotherapy response and prognosis of bladder cancer patients. ResponseScore can also reflect features such as tumor microenvironment, antitumor immunity, and natural killer cell function status in bladder cancer. PSMB9 was identified as a significant gene for prognosis. High expression of PSMB9 can inhibit bladder cancer cell migration and proliferation while increasing cell apoptosis.

Introduction

Bladder cancer is the most common tumor in the urinary system and the 13th leading cause of cancer death worldwide [1]. Bladder tumors are the sixth most common type of tumor [1]. Treatment options vary depending on disease stage and individual differences, such as surgical resection [2], chemotherapy based on different stages [3], neoadjuvant chemotherapy [4], and immunotherapy [5]. The therapeutic effect of bladder cancer depends to some extent on the clinical and pathological factors of patients, such as TNM staging and tumor grading, which are prognostic indicators. However, these clinical and pathological characteristics are not sufficient to predict patient prognosis and produce considerable differences in similar grading or staging [4]. This is essentially related to the heterogeneity of bladder tumors [6]. Therefore, there is an urgent need to discover biomarkers that can be used for better diagnosis, risk stratification, and clinical management of bladder cancer patient. Currently, the emergence of immunotherapy provides a new approach for cancer treatment. Currently, immunotherapies used for bladder cancer mainly include PD-L1 inhibitors and other immune checkpoint inhibitors. They provide first-line treatment options for patients with metastatic bladder cancer who cannot tolerate platinum-based chemotherapy. They are also used for maintenance therapy in patients who have achieved objective remission after platinum-based chemotherapy [7]. According to reports, avelumab (a PD-L1 inhibitor) can significantly prolong overall survival time for patients who have achieved objective remission. The median survival time of patients treated with avelumab is 21.4 months, while that of patients who have not received treatment is only 14.3 months [8]. A recent meta-analysis has also confirmed the enormous value of immunotherapy in the combined treatment of bladder cancer [9]. However, the effectiveness of immunotherapy is not significant for some patients, which poses challenges for the clinical management of bladder cancer patients. Finding biomarkers that can evaluate the efficacy of immunotherapy has become one of the hotspots in bladder cancer and even tumor research.

Currently, pathway analysis, especially pathway-based gene expression profiling analysis, has shown good stability and excellent performance. As a topological feature, pathways reflect to some extent the intrinsic connections between different genes. Therefore, compared with simple gene expression analysis, pathway-based gene expression profiling analysis can more effectively identify biomarker genes related to disease phenotypes. Among them, Graph Neural Networks (GNNs) have received widespread attention due to their ability to capture hidden topological information in pathways [10]. GNNs as a class of deep learning methods have gradually been adopted in recent bioinformatics research. For example, Peng et al. [11] proposed a GNNs model for cancer gene prediction, while Rhee S et al. used GNNs for breast cancer subtype classification. Simultaneously, the emergence of several new algorithms, such as the Prairie Dog Optimization Algorithm [12], Dwarf Mongoose Optimization Algorithm [13], and Metaheuristic Optimization Algorithms [14], has provided more powerful tools for optimizing data during the modeling process. In this study, we used GNNs to extract transcriptomic features and pathway-related topological features to model and predict the immunotherapy response of bladder cancer while identifying key pathways associated with immunotherapy response.

Multi-omics analysis is an important research method in biomedicine. Based on the key pathways identified by GNNs, we further constructed a more concise immunotherapy response predictive factor and explore its underlying mechanism using multi-omics analysis. At the same time, we also identified key genes for immunotherapy response and prognosis of bladder cancer. This study aimed to discover the predictive factor and key genes reflecting immunotherapy response and prognosis of bladder cancer, provide new targets for targeted therapy, and provide new insights for bladder cancer treatment.

Materials and methods

Data acquisition

Gene expression data and related clinical information of TCGA-BLCA were downloaded from The Cancer Genome Atlas Program (TCGA, https://portal.gdc.cancer.gov/) database. The IMvigor210 dataset was obtained from the R package IMvigor210CoreBiologies. RNA sequencing data was normalized to transcripts per thousand base million (TPM). Gene mutation data was obtained from the Genomic Data Commons (GDC; https://portal.gdc.cancer.gov/) and analyzed using the “maftools” package. Copy number variation data was obtained from the UCSC Xena database (http://xena.ucsc.edu/). Gene Expression Omnibus (GEO) database was used to obtain GSE130001 and GSE146137 datasets, and only human samples were extracted to obtain single-cell RNA sequencing (scRNA-seq) data.

GNNs model construction and identification of key pathways

Bilin Liang et al. [15] proposed PathGNN for predicting long-term survival of cancer patients. The GNNs model in this study refers to its architecture, as shown in Supplementary Fig. 1. First, pathway and gene network information were obtained from Reactome database using graphite R package. Pathways with less than 15 genes or more than 400 genes and duplicate pathways were removed, and 927 pathways were retained for GNNs model construction. In total, 2439 pathways were downloaded from the Reactome database, while 927 pathways were retained for analysis. Here we used “GSVA” R package to calculate pathway enrichment scores for each sample. Then we utilized Pytorch Geometric module in Python to construct graph dataset as the input of GNNs model. In this way, the input for GNNs models includes pathway enrichment scores for all included pathways in each sample, as well as the gene expression profiles associated with each included pathway. Mathematical principles and specific mechanisms of GNNs model can be found in Supplementary Table 1. Model construction was based on Pytorch (version 1.8) and Pytorch Geometric (version 2.0.3). The IMvigor210 dataset was randomly divided into training and testing sets in a 4:1 ratio. The training set comprised 238 samples, while the test set contained 60 samples. Model training was based on IMvigor210 training set with immunotherapy response as prediction target for 150 epochs, using Adam algorithm for model optimization. Learning rate was set to 0.001 for first 100 epochs and 0.0005 for remaining 50 epochs. The determination of the parameters of the model was mainly based on the research of predecessors and the accumulated tuning experience of the technical team. Integrated Gradients (IG) algorithm in Captum module was used to explain importance of each pathway to prediction results, and IG scores were obtained to reflect importance of each pathway. The principles and processing procedures of the Integrated Gradients algorithm can be found in Supplementary Table 1. After z-transformation, pathways with the absolute value of z-score greater than 4 were selected to construct new simplified version of GNNs model. Here we chose 4 as the cutoff of z-score mainly based on the balance between model performance and the number of pathways included. This model was trained for 250 epochs, with learning rate of 0.001 for first 150 epochs and 0.0005 for remaining 100 epochs. Model and code for dataset preprocessing, model construction, and training were uploaded to GitHub (https://github.com/dongbianrichu/Graph-Neural-Networks-for-immunotherapy-response-in-bladder-cancer). Model performance was evaluated using area under the receiver operating characteristic curve (AUC) metric.

Screening and clustering analysis of relevant genes.

The genes involved in key pathways were extracted. Single-factor logistic regression was used to screen out significant genes (p < 0.05) related to immunotherapy response, and a total of 269 genes were obtained for subsequent analysis. Consensus clustering was performed using the R package “ConsensusClusterPlus” to divide patients into different subgroups based on the expression matrix of the included genes. Principal component analysis (PCA) was performed to evaluate the clustering effect.

Development and validation of immunotherapy response scoring model

Based on the IMvigor210 dataset, we used LASSO regression algorithm in R package “glmnet” to further determine the key genes related to immunotherapy response. Based on the linear combination of the expression values of key genes and LASSO regression coefficients as weights, we developed an immunotherapy response scoring model named responseScore. TCGA-BLCA was used as an external validation set to explore the performance of responseScore.

Survival analysis and pathway enrichment

The cohort was divided into high and low groups based on responseScore using the median as the boundary. Kaplan–Meier curves were plotted using R package “survminer”. Multivariate Cox analysis was performed with responseScore or responseScore grouping and various clinical factors. Nomograms for predicting overall survival of 1 year, 3 years, and 5 years were constructed based on the Cox analysis results using R package “rms”. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed on differentially expressed genes between high and low responseScore groups using the “clusterProfiler” [16] package.

Antitumor immune-related analysis

The R package “estimate” was used to calculate the immune score and stromal score of patients. The R package “xCell” was used to evaluate the infiltration of various immune cells in patients. Bagaev et al. [17] proposed 29 functional gene expression signatures (Fges) to represent the cellular and functional properties of the tumor microenvironment. Microsatellite instability (MSI) was calculated for each sample using R package “PreMSIm” [18] and divided into two categories: MSI-high (MSI-H) and MSI-low/microsatellite stability (MSS). Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/login/) was used to evaluate the effect of immune checkpoint inhibition (ICI) treatment on individuals and obtain individual Exclusion score, myeloid-derived suppressor cells (MDSC) score, tumor-associated macrophages (TAM) M2 subtype score, cancer-associated fibroblasts (CAF) score, etc.

Single-cell analysis

Low-quality cells were excluded based on the criteria of nFeature_RNA > 500, nFeature_RNA < 7000, and percent.mt < 15. All four samples were standardized using the SCTransform method, excluding mitochondrial interference factors. Additionally, the glmGamPoi package was employed to accelerate the standardization process. All 3000 genes with the most significant mutations identified by SCTransform were used to find anchor points through the FindIntegrationAnchors function. Based on these anchor points, we integrated the four samples using the IntegrateData function for subsequent analysis. The FindVariableFeatures function in R package “Seurat” was used to identify 2000 highly variable genes for subsequent principal component analysis (PCA)-based dimensionality reduction. The t-distributed stochastic neighbor embedding (t-SNE) method was used for visualization of single-cell clustering based on principal components (PCs) 1–30. Using R package “SingleR” and calling HumanPrimaryCellAtlasData as the reference database, we inferred the cell type of each single cell. The R package “inferCNV” can identify evidence of large-scale chromosomal copy number changes in cells. Using endothelial cells, NK cells, and tissue stem cells as reference cells and setting 0.1 as the cutoff value, we predicted copy number variations (CNVs) of epithelial cells. Epithelial cells with CNA_value > 0.7 and cor.estimate > 0.7 were defined as malignant epithelial cells. R package “Monocle 2” was applied to generate pseudotime trajectories of single cells. The SCENIC package can evaluate the transcription factor activity of each cell based on scRNA-seq data. In addition, we used R package “CellChat” to describe the interaction between different types of cells.

Experiment on functional verification of PSMB9

One human bladder epithelial immortalized cell line (SV-HUC-1) and four human bladder cancer cell lines (EJ, 5637, HT-1376, UMUC3) were cultured, and the expression of PSMB9 was detected by qPCR and Western Blot (WB). Empty vector and PSMB9 overexpression vector were constructed and transfected into human bladder cancer cells (5637), resulting in PCDH-NC group and PCDH-PSMB9 group. After 24 h of transfection, the expression of PSMB9 was detected by qPCR. After 48 h of transfection, the expression of PSMB9 was detected by WB. At the same time, scratch experiments were performed to detect cell migration of the two groups at 0/6/12/24 h after transfection. After 48 h of transfection, transwell was used to detect cell migration, flow cytometry was used to detect apoptosis and cell cycle of the two groups. After 12/24/48/72 h of transfection, CCK-8 was used to detect cell proliferation of the two groups. For details of the experiments, please refer to Supplementary Table 5.

Statistical analysis

Wilcoxon test was used for comparison of continuous variables. Chi-square test was used for intergroup comparison of categorical variables. Log-Rank test was used for intergroup analysis of survival rate. Delong’s test was used for comparison of AUC. Pearson or Spearman correlation coefficient was used for correlation analysis as appropriate. A p value less than 0.05 was considered statistically significant (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).

Results

Construction of GNNs model and development of immunotherapy response score

The overall flowchart of this study is shown in Fig. 1. In this study, we first trained a GNNs model with 927 pathways as input. After 150 epochs of training, the model had an AUC of 0.961 (95% CI 0.844–1.000) on the training set. At the same time, we obtained the IG score of 927 pathways, which are detailed in Supplementary Table 3. Among them, there were 11 pathways with the absolute value of z-scores greater than 4 (recorded in Supplementary Table 3), ranked by IG score from high to low: R-HSA-69620, R-HSA-159231, R-HSA-194138, R-HSA-399721, R-HSA-381070, R-HSA-6811434, R-HSA-74158, R-HSA-193807, R-HSA-975634, R-HSA-399719, and R-HSA-174143. The information of these 11 pathways was recorded in Supplementary Table 4. Then we constructed a GNNs model with 11 pathways as input. After 250 epochs of training, the AUC of the model on the training set was 0.785 (95%CI: 0.719–0.851, as shown in Fig. 2A), and the AUC on the validation set was 0.839 (95% CI 0.718–0.959, as shown in Fig. 2B). TIDE value as a common predictor for the effect of immunotherapy performed far worse than the GNNs model on IMvigor210 training (P < 1.0 × 10^–3) and validation sets (P = 0.023). In addition, we also explored the performance of the GNNs model on TCGA-BLCA dataset. Due to the lack of immunotherapy data in TCGA-BLCA dataset, we could only explore the correlation between model output and TIDE value. As shown in Fig. 2C, the TIDE value of response group predicted by the model was much lower than that of non-response group (P < 1.0 × 10^–3). The distribution of TIDE values between the two groups was significantly different (as shown in Fig. 2D).

Fig. 1

The overall flowchart of the study

Fig. 2

The performance of GNNs model for predicting immunotherapy response. A ROC curves of GNNs model and TIDE score for predicting immunotherapy response on the training set of IMvigor210 dataset, and comparison of their AUCs. B ROC curves of GNNs model and TIDE score for predicting immunotherapy response on the validation set, and comparison of their AUCs. C Boxplot showing the difference in TIDE values between immunotherapy response group and non-response group predicted by GNNs model on TCGA-BLCA dataset. D Distribution of TIDE values in response group and non-response group predicted by GNNs model on TCGA-BLCA dataset

Subsequently, we extracted the genes involved in the key pathways, totaling 643. Based on the IMvigor210 dataset, single-factor logistic regression analysis of immunotherapy response was performed, and the results are shown in Supplementary Table 5. Significant genes related to immunotherapy response were screened out, totaling 269 genes for subsequent clustering analysis. According to the relative change in area under consensus clustering cumulative distribution function curve (Fig. 3B), the k value of clustering was finally set to 2. The clustering situation is shown in Fig. 3A. The PCA distribution maps of two different clusters are shown in Fig. 3C. As shown in Fig. 3D, there was a significant difference on immunotherapy response between different clusters (P < 1.0 × 10^–3). The distribution of clinical variables such as gender between different clusters is shown in Fig. 3E. Therefore, we can conclude that the screened genes can largely distinguish patients’ immunotherapy response.

Fig. 3

The results of consensus clustering and LASSO regression, with the performance of responseScore for predicting immunotherapy response. A Consensus clustering results of IMvigor210 dataset. B Relative change in area under consensus clustering cumulative distribution function curve for k. C Two-dimensional principal component analysis plot showing the distribution of different clusters in IMvigor210 dataset. D Chi-square test results of immunotherapy response between different clusters in IMvigor210 dataset. E Sankey diagram showing the basic statistical information of clinical variables in each cluster. F The relationship between lambda value of LASSO regression and binomial deviance. Finally, 7 genes were included to construct immunotherapy response score, namely responseScore. G ROC curve of responseScore predicting immunotherapy response on IMvigor210 dataset. For comparison, currently prevalent predictive models for immunotherapy response, including models from the Wei Meng team, the Tianlei Xie team, IL4I1, and MTHFD2, were also presented. H Correlation coefficient map between responseScore and expression levels of 12 immune checkpoint genes on TCGA-BLCA dataset

Subsequently, LASSO regression was performed on the 269 included genes. The relationship between the lambda value of LASSO regression and binomial deviance is shown in Fig. 3F. Based on the balance between binomial deviance and the number of genes, a scoring system for immunotherapy response was finally constructed with 7 genes, named responseScore. With 10 times the amplification value of LASSO regression coefficient as weight, the calculation formula of responseScore is -12.940213*NFIA + 3.745066*COPA + 1.939527*RFC4 + 0.508840*MCM10 + 0.393219*EXO1 + 4.341316*PSMB9 + 0.832620*CHEK2. We found that the AUC of responseScore in predicting immunotherapy response on IMvigor210 dataset was 0.739 (95% CI 0.662–0.804, as shown in Fig. 3G). To validate the performance of the responseScore, we compared it with currently prevalent predictive models and biomarkers for immunotherapy response, including models from the Wei Meng team [19], the Tianlei Xie team [20], IL4I1 [21], and MTHFD2 [22]. The results (as shown in Fig. 3G) highlight the outstanding performance of our model. Moreover, to verify the scientificity of the construction method of responseScore, we randomly extracted 6 genes involved in 11 key pathways, constructed another score based on logistic regression, and compared the AUC of the two. As shown in Supplementary Fig. 2, the performance of responseScore is significantly better than that of the other score. In addition, on TCGA-BLCA dataset, responseScore was positively correlated with the expression of immune checkpoint genes such as PDCD1, CD274, and CTLA4 (as shown in Fig. 3H). Therefore, responseScore has demonstrated its good ability to distinguish immunotherapy response.

The relationship between responseScore, each gene and prognosis

Patients were divided into high and low responseScore groups using the median as the boundary. The distribution of the expression of related genes and clinical variables such as gender between high and low groups on IMvigor210 dataset are shown in Fig. 4A. In K-M analysis of the entire IMvigor210 dataset, there was a significant difference in survival rate between high and low groups (P < 1.0 × 10^–3, as shown in Fig. 4E), which further stratified analysis also showed its extensive ability to distinguish prognosis (Fig. 4B). Multivariate Cox analysis showed that responseScore grouping was an independent prognostic factor regardless of clinical factors such as gender, race, history of tobacco use, tumor size, etc. (as shown in Fig. 4C). According to the results of multivariate Cox analysis (Fig. 4F), we constructed a nomogram (Fig. 4G) and found that its AUC for predicting 1-, 3-, and 5-year survival was 0.717 (95% CI 0.651–0.786), 0.699 (95% CI 0.635–0.774), and 0.736 (95% CI 0.631–0.819, as shown in Fig. 4H), respectively, which showed superior predictive performance. In this section, we also conducted a systematic review of prognostic models published by other research teams (as shown in Supplementary Table 6). By comparing them, we found that the overall performance of our prognostic model in this study surpasses that of most existing models.

Fig. 4

ResponseScore serves as an independent factor in survival analysis. A Heatmap shows the distribution of gene expression and clinical variables in high and low responseScore groups of IMvigor210 dataset. B Kaplan–Meier plot of high and low groups in IMvigor210 dataset and various stratified analysis results. C Forest plot shows the results of multivariate Cox regression with responseScore grouping and various clinical factors as variables in IMvigor210 dataset. D Forest plot shows the results of multivariate Cox regression with responseScore grouping and selected clinical factors as variables in TCGA-BLCA dataset. E Kaplan–Meier plot of high and low groups in TCGA-BLCA dataset. F Forest plot shows the results of multivariate Cox regression with responseScore and various clinical factors as variables in TCGA-BLCA dataset. G Nomogram constructed based on responseScore and various clinical factors in TCGA-BLCA dataset. H ROC curve of nomogram for predicting 1, 3, and 5-year survival

We further studied the relationship between each gene involved in responseScore and prognosis. The results of single-factor Cox analysis of each gene on IMvigor210 dataset are shown in Supplementary Fig. 3A, and all 7 genes are significant prognostic factors. Except for the high expression of NFIA, which is associated with poor prognosis, the expression of other genes is positively correlated with good prognosis. The K-M analysis of each gene (Supplementary Fig. 3B) further verified the ability of each gene to stratify prognosis using the median of gene expression as the cutoff point. The results of single-factor Cox analysis on TCGA-BLCA dataset are shown in Supplementary Fig. 3C. Using the best threshold determined by R package “survminer” as the cutoff point for grouping, the survival rate of high-expression group of PSMB9 and CHEK2 was significantly higher than that of low-expression group (P = 0.035 and 0.045, respectively, as shown in Supplementary Fig. 3D), while NFIA was the opposite (P < 1.0 × 10^–3, as shown in Supplementary Fig. 3D).

Immune landscapes of different responseScore groups

Tumor microenvironment is a complex environment composed of tumor cells, immune cells, and stromal cells, which affect the invasion and metabolism of tumor cells [23]. The result of ESTIMATE analysis found that the immune score of the high responseScore group was significantly higher than that of the low group (P = 0.003, Fig. 5A), indicating that the infiltration degree of immune cells in the high group was higher. The stromal score of the high group was not significantly different from that of the low group (P = 0.890, Fig. 5B). Figure 5C shows the infiltration of immune cells evaluated by XCELL algorithm. The abundance of activated dendritic cells, B cells, CD4 + memory T cells, CD8 + naive T cells, CD8 + cells, macrophages M1, NK cells, neutrophils, Tgd cells, Th1 cells, Th2 cells and regulatory T cells in patients with high responseScore was higher than that in patients with low score (Fig. 5C). In addition, we also performed correlation analysis between immune cell infiltration and responseScore, and found that the infiltration of 23 immune cells, such as activated dendritic cells, B cells, CD4 + memory T cells, CD8 + naive T cells, macrophages M1, NK cells, Th1 cells, Th2 cells and regulatory T cells, were positively correlated with responseScore (Fig. 5D). We compared the scores of 29 Fges between the two groups and found that the high responseScore group had significantly higher score of anti-tumor immune signature, such as MHC II, Antitumor cytokines, Co-activation molecules, NK cells, Checkpoint molecules and Effector cells (Fig. 5E). IMvigor210 dataset classified samples into three immune subtypes: desert, excluded and inflamed. As shown in Fig. 5F, the proportion of excluded and desert phenotypes in high group (40% and 22%) was lower than that in low group (52% and 34%), while inflamed phenotype (38%) was higher than that in low group (13%). There was a significant difference in immune phenotype between the two groups (P < 1.0 × 10^–4).

Fig. 5

Immune landscapes of different responseScore groups. A Comparison of Immune Score between high and low responseScore groups in the IMvigor210 dataset. B Comparison of Stromal Score between the two groups. C Boxplot showing the comparison of the infiltration of various immune cells evaluated by the XCELL algorithm between the two groups. The items marked with red stars indicate statistically significant differences. D Correlation graph between responseScore and various immune cells. E Heatmap showing the comparison of the enrichment scores of 29 Fregs between high and low groups. The items marked with red stars indicate statistically significant differences. F Comparison of the distribution of three immune phenotypes between high and low groups in the IMvigor210 dataset. G Comparison of MSI distribution between high and low groups in the IMvigor210 dataset. H Comparison of responseScore between MSI-high group and MSI-low/microsatellite stability group. I Correlation graph between responseScore and expression levels of 12 immune checkpoint genes in the IMvigor210 dataset. J Boxplot showing the comparison of myeloid-derived suppressor cells (MDSC) score between high and low groups in the IMvigor210 dataset. K Comparison of tumor-associated fibroblasts (CAF) score between high and low groups. L Comparison of tumor-associated macrophages (TAM) M2 subtype score between the two groups. M Comparison of Exclusion score between the two groups

Immunotherapy response is closely related to microsatellite instability [24]. Microsatellite instability is caused by defects in the mismatch repair system, leading to errors in DNA replication. In many solid tumors, a close correlation has been found between MSI-H status and immunotherapy response. Currently, MSI has been clinically approved as a predictive factor for immunotherapy response. We explored the association between responseScore and microsatellite instability. As shown in Fig. 5G, the MSI-H and MSS/MSI-L subtypes accounted for 91% and 9%, respectively, in the high responseScore group, while they were 87% and 13%, respectively, in the low group. There was no difference in responseScore between MSI-H subtype and MSS/MSI-L subtype (P = 0.220, Fig. 5H), indicating that responseScore is an immunotherapy effect marker different from microsatellite instability. The results of correlation analysis between responseScore and 12 immune checkpoint genes (Fig. 5I) showed that three common immune checkpoint genes PDCD1, CD274, and CTLA4 were positively correlated with the score, indicating that the higher the responseScore, the more effective the immune checkpoint inhibitor treatment may be. In addition, we observed that although there was no significant difference in MDSC score between the two groups (P = 0.170, Fig. 5J), the low group had higher CAF score (P = 0.010, Fig. 5K), TAM M2 subtype score (P < 1.0 × 10^–4, Fig. 5L), and therefore higher Exclusion score (P < 1.0 × 10^–3, Fig. 5M).

Genetic changes

We further studied the genetic variations of high and low responseScore groups in TCGA-BLCA dataset, including somatic mutations and CNVs. The waterfall plots of mutation genes in the two groups are shown in Supplementary Fig. 4A and 4B, respectively. Somatic mutations occurred in 87.88% of low-group samples and 93.10% of high-group samples. The top 6 genes with the highest mutation frequency were the same in both groups, namely TP53, TTN, KMT2D, MUC16, ARID1A, and KDM6A. Supplementary Fig. 4C shows the difference in somatic gene mutations between the two groups. TP53 (P < 1.0 × 10^–3), MAGI3 (P < 1.0 × 10^–3), EPG5 (P < 1.0 × 10^–2), HECW1 (P < 1.0 × 10^–2) and other mutations were significantly higher in the high responseScore group than in the low group, while FGFR3 (P < 1.0 × 10^–2), RNF213 (P < 1.0 × 10^–2), IFT122 (P < 1.0 × 10^–2) and SOX5 (P < 1.0 × 10^–2) mutations were higher in the low group than in the high group. Supplementary Fig. 4D and 4E respectively show the basic consistency of tumor-related signaling pathways enriched by mutation genes between high and low groups, including RTK-RAS, NOTCH, WNT, Hippo, PI3K, MYC, TGF-Bata, TP53 and NRF2 signaling pathways. As for the tumor mutation burden (Supplementary Fig. 4F), the burden of the high group was significantly higher than that of the low group (P < 1.0 × 10–3). The frequency distribution maps of copy number variations in the two groups are shown in Supplementary Fig. 4G and Supplementary Fig. 4H.

Functional analysis

To further study the underlying mechanism of responseScore, relevant pathways were analyzed through GO and KEGG databases. The main cellular locations/compartments involved in GO analysis were endoplasmic reticulum lumen, leading edge membrane, intermediate filament cytoskeleton, receptor complex and so on (Supplementary Fig. 5A). The main pathways involved in KEGG pathway enrichment analysis included Wnt signaling pathway, Focal adhesion, PI3K-Akt signaling pathway, MAPK signaling pathway, DNA replication, Mismatch repair, etc. (Supplementary Fig. 5B).

Single-cell analysis

After strict screening, 7725 cells were obtained for subsequent analysis. To determine the cell type, we set the resolution to 0.5 and identified 11 cell clusters. Using SingleR to infer cell types, we finally determined four types of cells: epithelial cells, endothelial cells, natural killer cells (NK cells), and tissue stem cells, as shown in Fig. 6A. Supplementary Fig. 6 displays the confidence of cell type identification performed on this dataset. The responseScore of NK cells was higher than that of other types (Fig. 6B). We then extracted NK cells for analysis and divided them into high-score and low-score subgroups using the median as the boundary. To better understand the functional state of NK cells, we used Monocle 2 to determine their developmental trajectory. Pseudo-time analysis showed that NK cells were divided into different stages on the time trajectory (Fig. 6C), that is, different functional states. ResponseScore grouping can represent different functional states. The results in Fig. 6D also show that responseScore reflects the functional state of NK cells to some extent. We extracted differential genes between the two subgroups for enrichment analysis (Fig. 6E), and the results showed that biological processes enriched in high-score NK cell subgroup included positive regulation of leukocyte cell_cell adhesion, positive regulation of cell activation, regulation of cell_cell adhesion and leukocyte cell_cell adhesion.

Fig. 6

The association of responseScore with NK cells. A t-SNE plot showing the clustering information and cell types of scRNA-seq data. B Comparison of responseScore between various cell types. C Results of pseudotime analysis. The distribution of responseScore group (left), state (middle), and pseudotime (right) reflects a coherent and unified pattern. Left: Natural killer cells (NK cells) are divided into two different responseScore groups and are color-coded differently. Middle: In the developmental trajectory, NK cells are divided into different states. Right: The trajectory of pseudotime is shown by the gradient from black to blue. D Distribution of responseScore of NK cells on the trajectory of pseudotime. E Results of GO analysis. Highly expressed genes in the high-scoring subpopulation of NK cells were extracted and subjected to GO analysis, which showed the pathways enriched in the high-scoring subpopulation

Analysis of the activity of transcription factors in two subgroups showed that the activity of XBP1, ETS1, FOS, STAT3, IRF8 and other transcription factors in the high-score subgroup of NK cells was higher than that in the low-score subgroup (Supplementary Fig. 7A). To study the interaction between NK cells and bladder cancer cells, we used R package “inferCNV” to infer malignant epithelial cells. With CNA_value > 0.7 and cor.estimate > 0.7 as parameters, we identified 2000 malignant epithelial cells, as shown in Supplementary Fig. 7B. The number and intensity of interactions between different types of cells are shown in Supplementary Fig. 7C and Supplementary Fig. 7D respectively. The signal output of two different subgroups of NK cells to other types of cells is shown in Supplementary Fig. 7E. It could be seen that the number and intensity of interactions between high-score subgroup of NK cells and bladder cancer cells were higher than those between low-score subgroup. The high-score subgroup had higher incoming signaling patterns such as ICAM, CXCL, MHC I (Supplementary Fig. 7F) and outgoing signaling patterns such as ITGB2 (Supplementary Fig. 7G) than the low-score subgroup.

Functional experimental results of PSMB9

To further explore the genes involved in responseScore and find potential bladder cancer treatment targets, we conducted biological experiments on PSMB9 and explored its biological effects. The results of Fig. 7A and B showed that PSMB9 gene was less expressed in human bladder cancer cells than in bladder epithelial immortalized cells (SV-HUC-1). Figure 7C and D confirmed the success of plasmid transfection experiment. PCDH-PSMB9 group showed significantly higher expression of PSMB9 compared to PCDH-NC group. Then, we conducted different functional experiments and compared the two groups. Scratch test (Fig. 7E) and transwell results (Fig. 7F) showed that overexpression of PSMB9 inhibited bladder cancer cell migration and the osmotic changes of bladder cancer cells after treatment of 5637 cells. CCK-8 experiment results (Fig. 7G) also confirmed the inhibitory effect of PSMB9 overexpression on the proliferation ability of bladder cancer cells. Shown in flow cytometry result (Fig. 7H), the apoptosis of 5637 cells was increased after overexpression of PSMB9. Flow cycle experiment results (Fig. 7I) showed that the proliferation of 5637 cells was slowed down after overexpression of PSMB9.

Fig. 7

Functional experimental results of PSMB9. A qPCR detection of PSMB9 mRNA expression in bladder epithelial immortalized cells (SV-HUC-1) and human bladder cancer cells (EJ, 5637, HT-1376, UMUC3). B The expression of PSMB9 protein in bladder epithelial immortalized cells (SV-HUC-1) and human bladder cancer cells (EJ, 5637, HT-1376, UMUC3) was detected by Western Blot. C qPCR detection of PSMB9 mRNA expression in PCDH-PSMB9 group and PCDH-NC group. D The expression of PSMB9 protein in two groups was detected by Western Blot. E The effect of overexpression of PSMB9 on cell proliferation was detected by scratch test. The area between the two horizontal lines represents the location of the scratches. F Transwell detected the effect of overexpression of PSMB9 on cell migration. G Effect of overexpression of PSMB9 on cell proliferation detected by CCK-8. H The effect of overexpression of PSMB9 on apoptosis was detected by flow cytometry. I Effect of overexpression of PSMB9 on cell cycle was detected by flow cytometry

Discussion

This study first determined the key pathways for predicting the immunotherapy response of bladder cancer through convolutional neural networks method. Pathways are interactions between genes and other substances that exhibit specific biological functions and are usually viewed as a type of graphical data. GNNs are a class of deep learning methods that perform better than convolutional neural networks or other machine learning methods on some graph data because they can capture the topological information of graph data. In this study, GNNs were used to model and predict the immunotherapy response of bladder cancer, achieving excellent results. The GNNs model with 927 pathways as input had an AUC of 0.961 (95% CI 0.844–1.000) on the training set. To improve the generalized predictive ability of the model, we selected 11 key pathways based on the IG score and constructed a GNNs model with 11 pathways as input. The model had an AUC of 0.785 (95% CI 0.719–0.851) on the training set and an AUC of 0.839 (95% CI 0.718–0.959) on the validation set, with excellent predictive performance. The model not only used gene set information but also effectively used topological information hidden in pathways, demonstrating the great application prospects of pathway analysis and GNNs in predicting immunotherapy response and promoting the application of pathway-based analysis in the biomedical field. In addition, we have identified 11 key pathways in immunotherapy response for bladder cancer. R-HSA-69620 is the cell cycle checkpoints pathway, which monitors the progress of the cell cycle to ensure genome integrity, chromosome separation fidelity, and mutual dependence between S phase and mitosis. In fact, the inability to monitor cell cycle progression and constrain genome integrity may be a prerequisite for molecular evolution required for tumor development. Therefore, cell cycle checkpoints pathway is crucial for tumor biological processes. This pathway is also closely related to immunotherapy response. Wei-Wei Shi et al. [25] constructed a scoring system based on this pathway that can distinguish immunotherapy response to some extent. R-HSA-159231 is the transport of mature mRNA derived from an transcript pathway, which involves the process of mRNA transfer from the nucleus to the cytoplasm and plays an important role in correcting RNA processing. R-HSA-194138 is Signaling by VEGF, where vascular endothelial growth factor (VEGF) is a key regulator of angiogenesis. In tumor progression, activation of VEGF pathway promotes tumor angiogenesis, tumor growth and metastasis. From the perspective of antitumor immune, VEGF can induce MDSC growth and enhance its immunomodulatory activity [26], thereby changing tumor microenvironment and affecting the effect of immune checkpoint therapy. It has been reported that anti-VEGF drugs can alleviate immunosuppression, allowing CD8 + T cells to enter tumor tissue and enhance the effect of immune checkpoint blockade therapy [27]. R-HSA-399721 is glutamate binding, activation of AMPA receptors and synaptic plasticity pathway. R-HSA-381070 is IRE1alpha activates chaperones pathway. This study identified key pathways that affect bladder cancer immunotherapy response, providing ideas for improving the efficacy of immune therapies such as immune checkpoint blockade.

This study derived a linear and computationally efficient predictive factor, namely responseScore, based on the GNNs model, and explored its value in tumor immunity and tumor treatment. The AUC of responseScore in predicting immunotherapy response in the IMvigor210 dataset was 0.739, and it was positively correlated with the expression of immune checkpoint genes such as PDCD1, CD274, and CTLA4 in the TCGA-BLCA dataset, demonstrating its good ability to distinguish immunotherapy response. Meanwhile, this study explored the value of responseScore in prognosis. K-M analysis and Cox analysis on the IMvigor210 dataset and TCGA-BLCA dataset showed that responseScore was an independent prognostic factor. We also constructed a nomogram to predict bladder cancer prognosis and found that its AUC for predicting 1-, 3-, and 5-year survival was 0.717, 0.699, and 0.736 (Fig. 4H), respectively, with good prediction performance. By comparing it with models and biomarkers previously published by other teams, we observed that responseScore has significant value in predicting both immunotherapy response and prognosis in bladder cancer.

This study then explored the intrinsic mechanism of responseScore from the perspectives of genetic variation, immune microenvironment, and single-cell analysis. The analysis results demonstrate a broad association between responseScore and the tumor microenvironment. Single-cell analysis further reveals a close connection between responseScore and NK cells within the tumor microenvironment, providing valuable insights into the relationship between the responseScore and immunotherapy response.

In the results of single-cell analysis, the responseScore of NK cells was higher than that of other cell types (Fig. 6B), further indicating that responseScore was positively correlated with the infiltration of NK cells. The results of Monocle 2 showed that responseScore grouping represented different functional states of NK cells to some extent (Fig. 6C), among which high-scoring NK cell subpopulations enriched biological processes such as positive regulation of leukocyte cell_cell adhesion, positive regulation of cell activation, leukocyte cell_cell adhesion (Fig. 6E), etc. Therefore, we personally speculate that the high- and low-scoring groups represent two functional states of NK cells: high-efficiency state and low-efficiency state (or resting state). The results of transcription factor activity analysis (Supplementary Fig. 7A) showed that the activity of XBP1, ETS1, STAT3, IRF8 and other transcription factors in the high-scoring subpopulation of NK cells was higher than that in the low-scoring subpopulation. According to reports, XBP1 can positively regulate the cytotoxicity of NK cells against tumor cells and help NK cells survive through anti-apoptotic mechanisms [28]. ETS1 can induce the expression and function of various NK cell receptors (NKR), which are necessary for NK cell-mediated cytotoxicity; ETS1-deficient NK cells have decreased degranulation ability, and cannot kill target cells because they cannot respond to NKR ligand activation [29]. STAT3 can increase the expression of NKG2D in NK cells and promote NK cell-mediated target cell killing in some tumor models [25]. IRF8 regulates TLR9 expression in NK cells by directly binding to promoters to activate NK cells; IRF8 leads to a nearly tenfold decrease in TLR9 expression in NK cells, leading to activation failure [30]. These transcription factors are all closely related to the cytotoxicity effect of NK cells on target cells, indirectly verifying that high responseScore represents the cytotoxicity state of NK cells. The results of cell communication showed that the number and strength of interactions between high-scoring subpopulations of NK cells and bladder cancer cells were higher than those between low-scoring subpopulations (Supplementary Fig. 7D), and incoming signaling patterns such as ICAM and CXCL (Supplementary Fig. 7F) were higher than those between low-scoring subpopulations. Outgoing signaling patterns such as ITGB2 (Supplementary Fig. 7G) were also higher than those between low-scoring subpopulations. ICAM-1 can form stable contacts between NK cells and endothelial cells to play an important role in homing NK cells to solid tumors [31], while promoting complex formation between NK cells and tumor cells, activating NK cells to kill tumor cells [32]. Chemokines CXCL-9, CXCL-10, and CXCL-11 enhance the migration of NK cells and activate their effects by acting on CXCR3 [33]. MHC class I molecules bind to killer cell immunoglobulin-like receptors (KIRs) on NK cells to produce downstream signals that inhibit cell killing. This is crucial for the function of NK cells. Tumor cells trigger NK cell killing because they lack MHC class I molecules [34]. These aspects all verify from different angles that the high-scoring subpopulation of NK cells is in a high-efficiency state of killing target cells. After receiving immune checkpoint blockade therapy, inhibitory pathways such as PD-1 are blocked, and patients with high-efficiency state NK cells in the tumor microenvironment have significant natural efficacy.

These genes involved in responseScore also showed the ability to distinguish prognosis (see Supplementary Fig. 3A, Supplementary Fig. 3C and Supplementary Fig. 3D), demonstrating their potential as tumor markers and bladder cancer treatment targets. The discriminative ability of these genes in prognosis also explains to some extent the sensitive reflection of responseScore to prognosis. K-M analysis of transcriptome data (IMvigor210 dataset and TCGA-BLCA dataset) showed that the survival rate of the high-expression group of PSMB9 was significantly higher than that of the low-expression group. To find potential bladder cancer treatment targets, we experimentally studied the biological effects of PSMB9 in bladder cancer. The results showed that overexpression of PSMB9 can inhibit bladder cancer cell migration, cell proliferation, and cell invasion changes while increasing bladder cancer cell apoptosis. Therefore, high expression of PSMB9 can reduce the malignant behavior of bladder cancer cells. PSMB9 gene encodes one of the catalytic subunits of proteasomes, which are involved in the metabolism of macromolecules related to proteasomes and may be involved in antigen processing and presentation, making PSMB9 an attractive target gene for cancer research. PSMB9 has been revealed to have anticancer effects in skin melanoma [35], triple-negative breast cancer [36, 37], and liver cancer [37]. It also has been confirmed to be highly associated with prognosis and immune responses to checkpoint inhibitors in skin melanoma [38]. Yutao Wang and colleagues reported that PSMB9 may enhance CD8 + T cell infiltration by upregulating the expression of immunoproteasome 20S core subunits in bladder cancer [39]. This positive regulation affects MHC-I surface expression and antigen presentation, ultimately promoting CD8 + T cell infiltration [40]. Notably, CD8 + T cell infiltration is highly correlated with the immunotherapy response observed in bladder cancer. Nasser A Elhawary et al. [41] found that the single nucleotide polymorphism of PSMB9 can increase the risk of bladder urothelial carcinoma, indirectly revealing the anti-cancer effect of PSMB9. This study directly proved the anti-cancer biological effect of high expression of PSMB9 through transfection experiments with cancer cell lines. The PSMB9 gene has potential as a bladder cancer suppressor and can be applied to clinical treatment.

In conclusion, this study constructed a GNNs prediction model for the immunotherapy response of bladder cancer, identified key pathways, and derived an easily calculated predictive factor based on the GNNs model. The application value and intrinsic mechanism of responseScore were explored from multi-omics aspects. Finally, one of the genes of responseScore, PSMB9, was identified as a significant gene and specifically studied for its biological effects on bladder cancer.

Nevertheless, this study has some limitations. First, the construction and validation of the prognosis model are based on TCGA and IMvigor210 datasets and have certain biases; second, the GNNs prediction model and responseScore of this study still need real-world data to verify their clinical utility; thirdly, this study is mainly based on modeling transcriptome data, primarily considering the widespread application of transcriptome technology and the ease of obtaining corresponding datasets. Then, the process from gene transcription to function realization is a very complex process with many intermediate steps, so it is necessary to consider these biological processes. Incorporating proteomics and metabolomics data into model analysis may enhance the predictive accuracy of models and provide further insights into the mechanisms underlying bladder cancer. According to reports [42, 43], Graph Neural Networks have been applied to model proteinomics, metabolomics, and other data, playing a role in predicting tissue phenotypes and developing cancer drugs. In the future, we will expand the input range of our GNNs model, striving to include data from methylation, proteomics, metabolomics, etc., to achieve higher prediction accuracy, and to identify more meaningful bladder cancer pathways and targets. Simultaneously, we are applying GNNs to the field of pan-cancer research to uncover broader pathways and mechanisms underlying cancer. Nevertheless, this study still has great value.