Survival differences existed between patients in cluster 1 and 2

Firstly, the expression of m6A methylation regulatory-related and glucose metabolism-related genes was retrieved in the TCGA-BLCA dataset. A total of 326 m6A-related glycolytic genes were identified by Spearman correlation analysis with |Cor|> 0.5 and p < 0.05 as the threshold of significance in BLCA samples (Additional file 1: Table S1), of which 70 were prognosis-related m6A glucose metabolism genes, as identified by univariate Cox analysis. Among these, IGF2BP3, IGF2BP2, SLC16A1, DSC2, TGFBI, P4HA2, ANXA1, and CALU were negative prognostic factors, and CDK10 was the most positive prognostic factor (Additional file 4: Table S4). Based on the expression pattern of prognosis-related m6A glucose metabolism genes, BLCA were divided into two subtypes by consensus clustering analysis (Fig. 1A, B).

Fig. 1

Survival differences existed between patients in cluster 1 and 2. A, B Results of consistent cluster analysis used to divide BLCA patients into two subtypes. C KM curves of the two subtypes (p = 0.00738). D Clinical features of the two subtypes. **P < 0.01, ***P < 0.001

The KM curves indicated a significant difference in survival between the two subtypes of patients (p = 0.00738); cluster 2 patients had a better prognosis than that of cluster 1 patients (Fig. 1C). In addition, four clinical features, including gender, pathological T, M, and tumor stage, were significantly different between subtypes (Fig. 1D).

34 target genes were associated with metabolism-related pathways in BLCA

There were 7119 upregulated and 4,216 downregulated genes between 411 BLCA and 19 HC samples (Fig. 2A, B), and 22 upregulated and 66 downregulated genes between cluster 1 and cluster 2 samples (Fig. 2C, D). Subsequently, a total of 34 genes associated with BLCA occurrence and development (target genes) were screened by intersecting the differential expressed genes (DEG)1 between BLCA and HC samples and DEG2 between cluster 1and 2 (Fig. 2E, Additional file 5: Table S5).

Fig. 2.

34 target genes were associated with metabolism-related pathways in BLCA. A, B DEGs between 411 BLCA and 19 HC samples. C, D DEGs between cluster 1 and cluster 2. E Screening of m6A-related glycolytic genes by intersecting DEG1 and DEG2. F GO analysis of m6A-related glycolytic genes. G KEGG analysis of m6A-related glycolytic genes

GO functional analysis showed that the 34 target genes were related to several terms, including carbohydrate catabolism, glycolysis, ATP metabolism, and pyruvate metabolism (Fig. 2F, Additional file 6: Table S6). As for KEGG pathways, these genes were significantly associated with amino acid biosynthesis, carbon metabolism, glycolysis, and gluconeogenesis (Fig. 2G, Additional file 7: Table S7).

The survival risk model could be used to predict the prognosis of BLCA patients

Firstly, 34 intersection genes were submitted to univariate Cox regression analysis, and seven genes were obtained with p-value < 0.05 (Fig. 3A). Then, stepwise multivariate Cox regression analysis was used to identify two biomarkers (IP6K2 and PLA2G2F), both of which were positive factors (Hazard Ratio [HR] < 1) (Fig. 3B). The risk score was calculated for each patient, and the median risk score (0.961) was used as the cutoff to divide patients (n = 406, with complete survival data) into the high-risk group (n = 203) and low-risk group (n = 203). The results of KM and risk curves showed significant survival differences between the two groups (Fig. 3C, D). The AUC value of the ROC curve at 1, 3, and 5 years was all greater than 0.6, indicating that the survival risk model could be used as a prognostic model effectively (Fig. 3E). In addition, a significant negative correlation was found between disease risk and expression of IP6K2 and PLA2G2F.

Fig. 3

The survival risk model could be used to predict the prognosis of BLCA patients. A, B Results of Cox analysis identifying IP6K2 and PLA2G2F as positive factors (Hazard Ratio < 1). C, D Significant survival differences observed between the high- and low-risk groups. E ROC curve revealing the potential use of the survival risk model as a prognostic model. F–H Results of KM, risk and ROC curves in the validation cohort, consistent with the training cohort

GSE32894 was used as a validation dataset to verify applicability of the model. The results of the KM, risk, and ROC curves were consistent with those of the training dataset (Fig. 3F–H).

Risk score was an independent risk factor for BLCA patients

Univariate and multivariate Cox regression analyses were conducted to evaluate the clinical parameters and risk score to assess their prognostic value. The results revealed that gender, tumor stage, and pathological T stage were associated with risk score (Fig. 4A). Four clinical factors, including age, pathology T, N, and risk score, were significantly associated with patient survival (Fig. 4B, C). A nomogram based on these clinical factors was constructed to predict the patients’ 1-, 3-, and 5-year survival rates; the calibration curve showed that the slopes were closest to 1 (Fig. 4D, E). In addition, the nomogram model had a higher benefit rate than those of others, and the “number high risk” curves coincided with “number high risk with event” curves. These results indicated that the nomogram model has accurate predictive ability for BLCA (Fig. 4F, G).

Fig. 4

Risk score was an independent risk factor for BLCA patients. A Gender, tumor stage, and pathologic T were associated with risk score. B, C Age, pathology T, N, and risk score were negatively associated with patient survival. D, E The nomogram for predicting survival in BLCA patients (1-, 3- and 5-year). F, G Evaluation of the accuracy of prediction using the DCA curve (F) and Clinical Impact Curve (G). T (Tumor) refers to the extent of the primary tumor; N (Node) refers to the involvement of regional lymph nodes; M (Metastasis) indicates the presence of metastasis

Difference in function enrichment between the high- and low-risk groups

To investigate the potential molecular mechanisms of risk model, we performed GSEA enrichment analysis between the high- and low-risk groups. The GSEA results for the two groups are shown in Fig. 5A–D. T-cell activation, antigen receptor-mediated signaling pathway, B-cell differentiation, B-cell-mediated immunity, and B-cell receptor signaling pathway were enriched in the high-risk groups (Fig. 5A), whereas cellular glucuronidation, estrogen metabolism, miRNA-mediated gene silencing by inhibition of translation, sodium ion homeostasis, and translation repressor activity were enriched in the low-risk groups (Fig. 5B). The genes expressed in the high-risk groups were involved in antigen processing and presentation, cell adhesion molecules (CAMs), chemokine signaling, hematopoietic cell lineage, and JAK/STAT signaling (Fig. 5C). Meanwhile, genes expressed in the low-risk groups were involved in ascorbate and aldarate metabolism, drug metabolism, other enzymes, metabolism of xenobiotics by cytochrome p450, porphyrin and chlorophyll metabolism, and steroid hormone biosynthesis (Fig. 5D). The GSEA results revealed potential pathways or mechanisms that were activated during tumorigenesis and development, which can help us evaluate the prognosis of BLCA.

Fig. 5

The difference in function enrichment between the high- and low-risk groups A GO analysis showing enrichment of T-cell activation, antigen receptor-mediated signaling pathway, B-cell differentiation, B-cell mediated immunity, and B-cell receptor signaling pathway in the high-risk groups. B GO analysis showing enrichment of cellular glucuronidation, estrogen metabolism, miRNA-mediated gene silencing by inhibition of translation, sodium ion homeostasis, and translation repressor activity in the low-risk groups. C KEGG analysis showing enrichment of antigen processing and presentation and of CAMs in the high-risk groups. D KEGG analysis showing enrichment of ascorbate and aldarate metabolism, drug metabolism, and of other enzymes, in the low-risk groups

Significant differences were observed in immune microenvironment between the two risk groups

The CIBERSORT algorithm was used to investigate the differences in immune cell infiltration between the high- and low-risk groups. Seven immune cells, including naive B cells, plasma cells, activated memory CD4 + T cells, regulatory T cells, M1 macrophages, M2 macrophages, and activated dendritic cells, were found to have significantly different expression levels in the high- and low-risk groups. Among these, the immune infiltration of activated memory CD4 + T cells, M1 and M2 macrophages was significantly higher in the high-risk groups (Fig. 6A, B). IP6K2 expression correlated positively with M2 macrophages, while PLA2G2F expression correlated positively with M2 macrophages and negatively with activated dendritic cells (Fig. 6C). Meanwhile, risk score correlated positively with activated memory CD4 + T cells, M1 and M2 macrophages, and negatively with regulatory T cells and activated dendritic cells (Fig. 6D). In addition, 28 immune functions were significantly different between the high- and low-risk groups; among these, IP6K2 expression correlated positively with MHC class 1 and NK cell function (Fig. 6E, F). In summary, we speculated that immune microenvironment was associated with the occurrence and development of BLCA.

Fig. 6

Significant differences were observed in immune microenvironment between the two risk groups. A, B Seven differential immune cells were identified in the high- and low-risk groups. C Heat map of correlation between biomarkers and differential immune cells. D Correlation between risk scores and differential immune cells. E Box plot of differential immune functions between the high- and low-risk groups. F Heat map of correlation between biomarkers and differential immune function. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

The risk score can predict the effectiveness of immunotherapy

As immunotherapy became a promising strategy for the treatment of tumors, we evaluated differences in immunotherapy response between high- and low-risk groups. The expression levels of seven immune checkpoints (CD274, CTLA-4, LAG-3, HAVCR2, PDCD1, PDCD1LG2 and TIGIT), 14 immune checkpoint inhibitors (LAG3, CD40LG, CD86, CD48, CD274, CD244, CD70, CD27, CTLA4, ICOS, CD28, TIGIT, BTLA, TNFRSF4) and 32 immune factors were significantly higher in the high-risk groups, while those of two immune checkpoint inhibitors (BTNL2 and TNFRSF25) and three immune factors (CCR7, CCR10, CCL18) were significantly higher in the low-risk groups (Fig. 7A–C). Furthermore, the risk score significantly negatively correlated with the IPS score (Fig. 7D), indicating that patients in the low-risk group were more likely to benefit from immunotherapy.

Fig. 7

The risk score can predict the effectiveness of immunotherapy. A–C Expression of seven immune checkpoints, 14 ICIs and 32 immune factors in the high- and low-risk groups. D IPS scores in the high- and low-risk groups. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

The risk score model can predict the drug sensitivity between high- and low-risk groups

Considering that the risk score model could be used to predict the prognosis of BLCA patients, we further investigated the relationship between risk score and chemotherapy resistance. The results revealed that there were significant differences in the IC50 values of 114 drugs between the high- and low-risk groups. (Additional file 8: Table S8), and we further analyze the differences in drug sensitivity between the high- and low-risk groups for 10 drugs, such as parthenolide, sunitinib, cyclopamine, bexarotene, etc. (Additional file 9: Fig. S9). Therefore, we believe that this risk score model can be utilized to predict the sensitivity of BLCA patients to chemotherapy drugs, which could improve the effectiveness of chemotherapy.

IP6K2 and PLA2G2F were downregulated in bladder cancer tissues and could regulate bladder cancer cell proliferation in vitro

IP6K2 and PLA2G2F expression in bladder cancer and adjacent tissues was detected using RT-qPCR. Expression levels of IP6K2 and PLA2G2F were significantly downregulated in cancer tissues compared with their expression levels in adjacent tissues (Fig. 8A).

Fig. 8

IP6K2 and PLA2G2F were downregulated in bladder cancer tissues and could regulate bladder cancer cell proliferation in vitro. A IP6K2 and PLA2G2F expression in bladder cancer and matched adjacent normal tissues detected by RT-qPCR. B Efficiency of IP6K2 and PLA2G2F knockdown verified by RT-qPCR. C The proliferative ability of IP6K2 and PLA2G2F knockdown cells was determined using CCK-8 assay. D The colony-forming ability of IP6K2 or PLA2G2F knockdown cells was determined using colony formation assay. *P < 0.05, **P < 0.01, ***P < 0.001

Effect of IP6K2 and PLA2G2F on cell proliferation in vitro was explored to determine their role in bladder cancer cell proliferation. Efficiency of IP6K2 and PLA2G2F knockdown was verified using RT-qPCR (Fig. 8B). CCK-8 assays demonstrated that IP6K2 knockdown significantly promoted cancer cell proliferation in 5637 and UM-UC-3 cells (Fig. 8C). Moreover, the colony-forming ability of 5637 and UM-UC-3 cells enhanced after IP6K2 or PLA2G2F knockdown (Fig. 8D).