Construction and validation of a metabolism-associated gene signature for predicting the prognosis, immune landscape, and drug sensitivity in bladder cancer
BMC Medical Genomics volume 16, Article number: 264 (2023)
Tumor Metabolism is strongly correlated with prognosis. Nevertheless, the prognostic and therapeutic value of metabolic-associated genes in BCa patients has not been fully elucidated. First, in this study, metabolism-related differential expressed genes DEGs with prognostic value in BCa were determined. Through the consensus clustering algorithm, we identified two molecular clusters with significantly different clinicopathological features and survival prognosis. Next, a novel metabolism-related prognostic model was established. Its reliable predictive performance in BCa was verified by multiple external datasets. Multivariate Cox analysis exhibited that risk score were independent prognostic factors. Interestingly, GSEA enrichment analysis of GO, KEGG, and Hallmark gene sets showed that the biological processes and pathways associated with ECM and collagen binding in the high-risk group were significantly enriched. Notely, the model was also significantly correlated with drug sensitivity, immune cell infiltration, and immunotherapy efficacy prediction by the wilcox rank test and chi-square test. Based on the 7 immune infiltration algorithm, we found that Neutrophils, Myeloid dendritic cells, M2 macrophages, Cancer-associated fibroblasts, etc., were more concentrated in the high-risk group. Additionally, in the IMvigor210, GSE111636, GSE176307, or our Truce01 (registration number NCT04730219) cohorts, the expression levels of multiple model genes were significantly correlated with objective responses to anti-PD-1/anti-PD-L1 immunotherapy. Finally, the expression of interested model genes were verified in 10 pairs of BCa tissues and para-carcinoma tissues by the HPA and real-time fluorescent quantitative PCR. Altogether, the signature established and validated by us has high predictive power for the prognosis, immunotherapy responsiveness, and chemotherapy sensitivity of BCa.
The prevalence and mortality rates of bladder cancer (BCa) are high. It is estimated that currently > 1.6 million individuals have BCa worldwide, with approximately 550,000 new cases annually . It is anticipated that the incidence of BCa will be extremely rapid with increasing life expectancy and the increase of risk factors such as smoking in the developing countries. Despite research in recent years has made remarkable progress in our understanding of cancer, the burden of it is increasing.
Tumor cells reprogram the metabolic pathway to meet the bioenergetic, biosynthetic, and redox demands of malignant cells . The Warburg effect shows that a unique metabolic phenotype was identified for cancer cells, featuring that compared to normal tissues, they are characterized by a high rate of glycolytic metabolism. Recent evidence suggests that metabolic deregulation may not only be a characteristic feature of human cancers but also a potential underlying factor . The model of metabolic genes have been reported in pancreatic cancer , prostate cancer , hepatocellular carcinoma , colon cancer , and other cancer. However, the effects of metabolic changes on BCa progression are unknown. Recently, the new biomarkers used in drug screen of immunotherapy were introduced into clinical practice, which has shown favorable prospects in treating to some degree. However, the reaction to immunotherapy remains poorly elucidated.
In the current study, we first analyzed metabolism-associated prognostic DEG genes in bladder cancer, then constructed and validated a new signature across multiple datasets, and finally explored the prognostic value of the model and its relationship with the clinicopathological features, immune microenvironment, and the therapeutic effect of chemotherapy agent and immunotherapy of patients with BCa. We also confirmed the mRNA and protein expression of interested model-related genes from the Human Protein Atlas (HPA) and 10 paired BCa tissues collected by us. The most valuable aspect of the research is that our developed model provides more information for the selection of clinical treatment options and prognostic prediction.
Materials and methods
From September 2021 to October 2022, a total of 10 cases of BCa and its matched adjacent non-tumor tissues were collected in the Second Affiliated Hospital of Tianjin Medical University (China). The final diagnosis of each patient was confirmed by histopathology. In this study, the written informed consent of each patient or their guardian was obtained following the guidelines approved by the Medical Ethics Committee of the Second Hospital of Tianjin Medical University.
Data collection and processing
The datasets used to construct or validate metabolism-relevant gene signatures for BCa come from five different platforms: Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/repository) and Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). According to the TCGA database, 406 BCa samples were downloaded, including mRNA expression profile data and clinical information. Gene annotation is made possible through the Ensemble database. The microarray expression data of GEO datasets such as GSE13507, GSE32894, GSE48075, GSE31684, GSE111636, and GSE176307 were quantified and positiveized. IMvigor210 dataset, a cohort of 348 MIBC patients treated with atezolizumab (PDL1 inhibitor), is also used to verify metabolism-related gene signatures.
Notably, in our single-arm phase 2 trial (term_id, TRUCE-01; registration number, NCT04730219), tislelizumab (200 mg) combined with low-dose nab-paclitaxel (200 mg) also preliminary confirmed clinical benefits and safety in the therapy of muscle-invasive urothelial bladder carcinoma patients. Therein, tislelizumab is a novel humanized monoclonal antibody programmed death receptor-1 (PD-1) inhibitor and shows a predictable and manageable safety/tolerability profile in patients with PD-L1+ UC.
Extraction of metabolism-related prognostic genes, identification of BCa subclasses, and its association with infiltration of immune cells
We used cancer versus normal differential analysis to find metabolism-associated genes differentially expressed (DEGs) in tumor tissue and tumor surrounding tissues, and set |log(fold change, FC)|> 1, false-discovery rate (FDR) < 0.05 as the threshold to construct the heat map. Then, we performed a univariate Cox regression analysis on these DEGs.
NMF was used to identify subtypes of molecules and obtained the molecular subtypes C1, and C2 for BCa, and the log-rank test was applied. Kaplan–Meier survival analysis was used to compare the overall survival curves (OS), progression-free survival (PFS), disease-free survival (DFS), and disease-specific survival (DSS) between C1 and C2.
To analyze the immune status of each sample, Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) was used to calculate the relative proportions of 22 immune cells in BCa patients, which include seven T cell types, naïve and memory B cells, plasma cells, NK cells, and myeloid subsets. CIBERSORT was applied to convert mRNA data into the proportions of infiltrating non-tumor cells in the tumor microenvironment using standard annotation files to organize gene expression characteristics. Afterward, correlation analysis was conducted between the above molecular typing and immune cell infiltration estimated by CIBERSORT.
Establishment and validation of the metabolism-related prognostic prediction model
Then, the TCGA BCa cohort was randomly divided into a training set (n = 285) and a testing set (n = 121) at a ratio of 7:3. Chi-squared tests were employed to compare baseline characteristics between the training set and external validation set. Ten prognostic DEGs were screened out using univariate Cox regression and the Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis. Genes with P < 0.05 have prognostic significance. Then we constructed a risk prediction model with the following formula: risk score = 0.267*PLOD1 + 0.167*SERPINB7 + 0.160*TSPAN7 + 0.291*HSD17B1 + 0.132*PYCR1 + 0.266*PDGFRA + 0.107*EGR1 + 0.537*AHCY + 0.121*ATP6V1B1 + 0.202*SCD. Next, we divided BCa patients into high-risk and low-risk groups according to the median risk score. Across TCGA-training set, internal TCGA-testing and multiple external GEO validation (including GSE13507, GSE32894, GSE48075, GSE31684, and IMvigor210), Kaplan–Meier (KM) survival and time-dependent receiver operational feature (ROC) curves were plotted through “survival”, and “timeROC” R packages to estimate the discrimination of the metabolism-related 10-gene signature. Then, calibration curves were also generated to assess the stability of the models. The metabolism-related gene model was established by referring to previous studies . Subsequently, the some further analytic strategy on this model, such as the pathway enrichment analysis, and immune infiltration-related analysis were also similar to a previously published researchs [8, 9].
Correlations between the model and clinical features
Kaplan–Meier survival analyses for BCa patients of the high- and low-risk groups were carried out in the different clinical subgroups. In addition, univariate and multivariate Cox regression analyses confirmed that risk score was an independent prognostic predictor of OS. A combined nomogram model was further constructed to predict the 1-, 3- and 5-year OS rates of BCa patients. This nomogram included the risk scores and common clinical information such as grade, gender, age, and stage.
Predictive combined nomogram model construction
Using the “rms” package of the R software, a calibration chart was generated to verify the predictive power of the nomogram. Nomogram accuracy was evaluated using ROC curve analysis. A decision curve analysis (DCA) was then performed to assess the clinical usefulness of the predictive model. Furthermore, we compared the predictive power of our risk model with the previous two human models, including the tumor immune dysfunction and exclusion (TIDE) TIDE and the tumor inflammation signature (TIS).
GO, KEGG, and Hallmark genesets enrichment analyses via GSVA method based on TCGA, GSE13507, and GSEA32894
We explored the GO, KEGG, and Hallmark enrichment differences between the high- and low-risk groups by GSVA algorithm based on the TCGA dataset [10,11,12,13]. Of them, GO analysis includes molecular function (MF), cellular component (CC), and biological process (BP). Pathways with p < 0.05 were considered significantly enriched. Additionally, two additional different datasets (i.e., GSE13507 and GSEA32894) were interrogated to identify enriched functional and pathway terms.
Correlation analyses between the metabolism-related gene model and the tumor microenvironment and expression of immune checkpoint genes
We also used ESTIMATE’s algorithm to evaluate immune and stromal fractions to reflect the ratio of immune cells to stromal cells. Pearson’s correlation was used to analyze the correlation between the risk score of the signature and above immune or stromal scores. We evaluated also the expression differences of the multiple immune checkpoint genes, such as PD-1, PD-L1, CTLA4, etc., between high- and low-risk groups. Immune cell infiltration of tumor samples was identified using TIMER 2.0 (cistrome.shinyapps.io/timer/) via the MCPCOUNTER, CIBERSORT, QUANTISEQ, TIMER, CIBERSORT-ABS, EPIC, and XCELL algorithms. Wilcox rank sum test was used to analyze the difference in the level of each immune cell between the two risk groups. Spearman rank correlation analysis was applied to analyze the correlation between risk scores and immune infiltrating cell scores.
Association analyses of the model with additional immunological characterization
The tumor immune dysfunction and exclusion (TIDE) algorithm  based on the simulation of tumor immune evasion mechanism was used to analyze the expression level of relatively more important immunological checkpoint genes (including PD1, PDL1, VISTA, TIGIT, and CTLA4, etc.), which have been identified and developed for the treatment of cancer. Moreover, the TIDE was employed to analyze multiple signatures to estimate tumor immune evasion, such as exclusion, or CAF signatures. The immunotherapy response of cancer patients was predicted by the above multiple TIDE scores via the Wilcox test. PD1 and CTLA4 immune therapy agents can be predicted using the immunophenoscore (IPS). In addition, we also analyzed associations of the model with tumor mutation burden (TMB) and immunophenoscore (IPS), which are all superior predictors of the response to immune checkpoint blockades (anti-PD1, anti-CTLA4, etc.). Hence, we could predict the immunotherapy response of each BCa patient using our model.
Screening potential chemotherapy agents based on our prediction model
First, we downloaded the drug sensitivity information from the CellMiner database to explore the connection between modeled gene expression levels and drug susceptibility. Drugs approved by the FDA or clinical trials were selected for analysis. Next, the “pRRophetic” package was applied to analyze and compare the half-maximal inhibitory concentration (IC50) between two risk score groups to Doxorubicin, Docetaxel, Paclitaxel, and Vinblastine.
Correlation analysis of ten modeled genes with clinical immunotherapy efficacy across four independent cohorts
We drew a boxplot of model-related gene expression levels in the responder vs. non-responder groups based on Imvigor210, GSE111636, GSE176307, and Truce01 datasets. Notely, we also conducted a comparison of model-related gene expression differences after vs. before the treatment using a paired Wilcoxon test based on our Truce01 sequencing data. The statistical difference was set at p ≤ 0.05.
The expression validation of the model genes by Human Protein Atlas (HPA) and quantitative real-time PCR (qRT-PCR)
The protein expression levels of the ten model genes were verified through the HPA databases (HPA: https://www.proteinatlas.org/), and the consistency between the transcriptome described above and proteome levels was observed.
RNA extraction, quantitative reverse-transcription PCR, and human protein atlas (HPA) analysis. Total RNA was extracted from 10 paired BCa tumors and adjacent tissues via E.Z.N.A.TM Hp total RNA Kit (OMEGA); and was reversed into cDNA using RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Rockford, IL, USA). Quantitative reverse-transcription PCR (qRT-PCR) was used to determine the relative expression of EGR1, PLOD1, and PYCR1 mRNAs in 10 BCa tissues compared to matched paracancerous tissues. The EGR1 primer sequences were: forward, 5’-GGTCAGTGGCCTAGTGAGC-3’; reverse, 5’-GTGCCGCTGAGTAAATGGGA-3’. PLOD1 primer: forward, 5’-AGACCAAGTATCCGGTGGTGT-3’; reverse, 5’- CTTGAGCACGACCTCATCCAA-3’. PYCR1 primer: forward, 5’-CTTCACAGCAGCAGGCGTC-3’; reverse, 5’-TCTCCTTGTTGTGGGGTGTC-3’. GAPDH was used as a control gene: forward primer, 5’-CGGAGTCAACGGATTTGGTC-3’; reverse primer, 5’-TTCCCGTTCTCAGCCTTGAC-3’. The final results were analyzed using the 2−ΔΔCT method.
R software version 4.1.0 was used for all statistical analyses. Using Wilcox’s test, we compared the two groups’ variables. Chi-square tests have been used to examine the association of risk groups with clinicopathological features. The survival data were assessed using the Kaplan–Meier curve. The ROC analysis was conducted using the R package time. An univariate and multivariate Cox regression analysis was conducted to evaluate independent prognostic factors. P ≤ 0.05 was considered to be statistically significant. Furthermore, p-value summaries were as follows: P > 0.05 (ns); P ≤ 0.05 (*); P < 0.01 (**); P < 0.001 (***); P < 0.0001 (****).
Differentially expressed metabolism-related prognostic genes in BCa
This study's general step-by-step process is illustrated in Fig. 1. We collected 1314 metabolism-associated genes from TCGA, and 406 metabolism-associated genes with differential expression were identified using the R project “Limma” (| logFC |> 1, FDR < 0.05) (Fig. 2A). Subsequently, a total of 68 OS-related genes were detected via the univariate Cox regression (P < 0.05).
Molecular typing based on metabolism-related prognostic DEGs
BCa patients were clustered into two distinct subtypes using the ‘NMF’ algorithm, which is an effective method for reducing dimension in cancer subtype identification (Fig. 2C). The optimal number of clusters was determined as k = 2 after the NMF rank survey (Fig. 2B; Supplementary Figure S1). The results showed that patients with the C1 subtype had a more favorable prognosis of OS, DSS, and PFS than those with the C2 subtype (Fig. 2D-G). To explore the differences in immune cell infiltration between the two subtypes, we also evaluated 22 kinds of immune cells by the CIBERSORT algorithm. There were significant differences in infiltration levels between the two groups for many types of immune cells. For example, cluster 1 presented higher infiltration of naive B cells, plasma cells, CD8 + T cells, follicular helper T cells, Tregs, and activated dendritic cells, and cluster 2 exhibited greater infiltration of memory B cells, CD4 + memory resting T cells, M0 Macrophages, M2 Macrophages, resting dendritic cells and neutrophils (Fig. 2H, I).
Construction and validation of a prognostic prediction model of metabolism-related genes
We randomly divided the TCGA BCa cohort into training and internal testing sets at a ratio of 7:3 (n = 286, and n = 120, respectively). Based on the chi-square test, there were no significant differences between the training and testing sets (Table 1). Then, the prognostic model of ten metabolism-related genes (PLOD1, SERPINB7, TSPAN7, HSD17B1, PYCR1, PDGFRA, EGR1, AHCY, ATP6V1B1, and SCD) was successfully constructed using the lasso and multivariate cox regression based on TCGA dataset (Fig. 3A, B). It can be seen that the C2 subtype had a higher risk score than the C1 subtype (Fig. 3C). BCa patients with high risk survived significantly shorter times than low-risk patients, regardless of whether in a training group, a test group, or an entire group (Fig. 3D). The performance of the KM survival curves was further validated in the external validation cohorts (GSE13507, and GSE32894) (Fig. 3G, J). Through the use of ROC curves, the model’s performance was evaluated. The area under the ROC curve (AUCs) of 1-, 3-, and 5-year were 0.752, 0.659, and 0.671, respectively, in the TCGA-entire set; 0.775, 0.687, and 0.705, respectively, in the TCGA-training set; 0.717, 0.629, and 0.623, respectively, in the TCGA-testing set; 0.750,0.670 and 0.648, respectively, in the GSE13507 validation set; and 0.723, 0.819, and 0.809, respectively, in the GSE32894 validation set (Fig. 3E, H, K). Additionally, the calibration curves of the prediction model were close to the standard curves in TCGA-entire, TCGA-training, TCGA-testing, and two external geo-validation cohorts, separately (Fig. 3F, I, L).
Correlations of the prognostic model established by 10 metabolism-related genes with clinical characteristics
Subgroup analysis results showed that the model had good predictive efficacy in different age, stages, gender, or stage_T groups (P < 0.05) (Supplementary Figure S2). Our univariate and multivariate Cox regression analysis exhibited that the risk score of the model can serve as an independent prognostic factor (p < 0.0001; Table 2). According to the TCGA cohort, the risk scores for each BCa sample were ranked from low to high, with a higher risk score indicating a higher risk of death (Fig. 4A, B). Expression levels of ten model genes showed a trend of increase with the risk score from low to high (Fig. 4C). To investigate the relationship between the gene signature and clinical features, we compared the risk scores between different clinical features using chi-square or Wilcoxon nonparametric tests across the TCGA and GEO datasets. As illustrated in Fig. 4D, patients older than 65 years have a higher risk score than patients younger than 65 years of age. Moreover, patients with high grades exhibited higher risk scores than patients with low grades. BCa patients with M1, N1-3, or T3-4 displayed a higher risk score than M0, N0, or T1-2, respectively. Furthermore, immunotyping of wound healing (C1) and interferon-gamma dominant (C2) may happen in patients with the high-risk score, while the inflammatory (C3), and lymphocyte depleted (C4) are more likely to appear in those with low-risk score (p < 0.05, Fig. 4D).
The risk scores, survival status, and heatmap of model gene expression were also performed using the GSE13507 and GSE32894 datasets, and the results were almost consistent with the above analysis findings of TCGA (Fig. 4E-J). The correlation results between risk score and clinicopathological parameters in GSE13507 (Fig. 4K) or GSE32894 datasets (Fig. 4L) were nearly consistent with the above analysis from the TCGA datasets. Additionally, we also compared the risk score among Sjodahl et al.’s five molecular classifications and found that the SCC-like group had a higher risk score (Fig. 4L).
Based on risk group, age and stage, we established a combined nomogram model for predicting the 1-, 3-, and 5-year OS incidences of BCa patients (Fig. 5A). We also utilized calibration plots to assess and confirm the concordance between the nomogram-predicted and actual 1-, 3-, and 5-year OS (Fig. 5B). The concordance index (C-index) for the nomogram was 0.698 (95% CI: 0.658–0.7372). Furthermore, Decision Curve Analysis (DCA) of 1-, 3-, and 5-year OS for the model further demonstrated our expectations. It was found that compared to a single clinical factor, the model and combined nomogram showed the highest clinical net benefit (Fig. 5C). Multivariate time-dependent ROC curve according to age, gender, grade, stage, nomogram, and risk score, is presented in Fig. 5D. Notely, compared to TIDE, Exclusion, and TIS prediction models from previous studies, the signature constructed by us presented a better sensitivity (Fig. 5E).
Function enrichment analysis of the model
To explore the underlying molecular mechanism of the metabolism-related gene model, we next conducted GSEA for GO, KEGG, and HALLMARK gene sets, and further identified differences in biological function and pathways between the high- and low-risk groups based on the TCGA, GSE13507, or GSE32894 cohort (Fig. 6, Supplementary Table S1). For GO analysis results, ‘immunoglobulin complex circulation’, ‘complement activation’, ‘immunoglobulin receptor binding’, ‘complement activation’, and ‘spliceosomal snRNP assembly’ were enriched in the high-risk group (Fig. 6A, D, G). KEGG pathways enrichment analysis revealed that ‘ECM receptor interaction’, ‘systemic lupus erythematosus’, ‘focal adhesion’, ‘melanoma’, ‘arrhythmogenic right ventricular cardiomyopathy arvc’, ‘neuroactive ligand-receptor interaction’, ‘regulation of actin cytoskeleton’, ‘leukocyte transendothelial migration’ and ‘taste transduction’ were particularly significant in the high-risk group (Fig. 6B, E, H). Among them, the most significantly enriched HALLMARK terms (‘epithelial-mesenchymal transition’, ‘angiogenesis’, ‘apical junction’, ‘myogenesis’, ‘hypoxia’, ‘inflammatory response’, etc.) is shown in Fig. 6C, F, I.
The external validation of the metabolism-associated gene signature based on additional three cohorts
To evaluate the validity and accuracy of the 10 gene prognostic model in clinical practice, in addition to the GSE13507 and GSE32894 datasets mentioned above, we further conducted a series of analyses for three additional external validation datasets, including GSE48075, GSE31684, and IMvigor210 (Figs. 7 and 8. Based on GSE48075 and GSE31684 cohorts, the findings from the KM, ROC, and calibration curve analysis were highly congruent with the analyses of the model performance described above (Fig. 7A-C, F–H). No correlation was detected between risk groups and several clinical traits (chi-square test, P > 0.05; Fig. 7D, I). It could be seen that patients with positive lymph nodes had a higher risk score than those with lymph node metastasis negative via the Wilcox test, which meant that a higher risk score indicated a worse prognosis; whereas risk scores were not significantly different between the presence and absence of FGFR3, P53, RB1, or Ras mutation, distant metastasis, high grade (Fig. 7E, J).
Notably, the IMvigor210 cohort, which is a cohort of 348 MIBC patients treated with Atezolizumab (PD-L1 inhibitor), was included to further evaluate the predictive power of the model in the BCa immunotherapy cohort. For the IMvigor210, the OS of high-risk patients is worse than that of low-risk patients according to the Kaplan–Meier curve (Fig. 8A, P < 0.05). Calibration curves for predicting the 10- and 20-month OS showed good performance (Fig. 8B). The predictive performance of OS risk scores was determined using a time-dependent ROC curve with an AUC of 0.580 at 10 months and 0.602 at 20 months (Fig. 8C). We also found that the risk score, survival status, and heatmap of model-gene expression in the validation set were nearly consistent with those in the training set (Fig. 8D). The chi-square test and Wilcox rank test were used to evaluate the correlation of the risk score group with existing clinical variables (Fig. 8E). There was a significant correlation between risk score/groups and treatment response or TCGA_subtype. It is worth mentioning that almost all patients with high risk were in the PD/SD cohort (83%) (Fig. 8F), and patients with PD/SD had higher risk scores than those with CR/PR (Fig. 8G). The analysis results of all these validation datasets also suggested that the 10 gene signature is an important prognostic predictor for BCa patients who received or did not receive immunotherapy.
The correlation analysis between risk score and immune score, stromal score, expression of immune checkpoint genes, immune-cell infiltration, or the immunotherapy response of BCa patients
Based on TCGA cohorts, we demonstrated that the risk score was positively associated with stromal score (p = 1.6e-12), and so is the immune score (p = 0.0072) (Fig. 9A). The results are also verified in GSE13507 and GSE32894 datasets (Fig. 9B, C). To further understand better the relationship between the risk score and immune response, we also explored the correlation between risk score and immune checkpoint gene expression. The results showed that in the high-risk group, CD27, CD80, CD86, TNFRSF9, PD-1, and CD48 checkpoint genes were all highly expressed in TCGA, GSE13507, and GSE32894 datasets (Fig. 9D, E, F).
The degree of immune cell infiltration in the tumor microenvironment affects tumor occurrence, progression, and therapeutic effect, especially immunotherapy. Based on the 7 immune infiltration algorithms (TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC), we investigated the correlation between risk scores and immune cell infiltration. Heat map of 45 significantly different immune cells in high-risk vs. low-risk groups, including CD8 + T cell, Neutrophil, Myeloid dendritic cell, M2 Macrophage, Tregs, Endothelial cell, Cancer-associated fibroblast, etc., was plotted in Fig. 10A (Wilcoxon test, P < 0.01). According to the Tumor Immune Dysfunction and Exclusion (TIDE)  tool, the results exhibited that the Exclusion score, CAF score, PD-1, CTLA4, IDO1, and TIGIT score of the high-risk group were higher than those of the low-risk group (Fig. 10B). To further strengthen this result, we applied the IPS method to predict the response of BCa patients to immune checkpoint blockade (PD-1 and CTLA4 therapy). The IPS score of the low-risk group was significantly higher than that of the high-risk group (Fig. 10E). Meanwhile, we performed a correlation analysis between risk scores and immune infiltrating cells using the Spearman test (Fig. 10C, D); this result is consistent with that of the above.
Screening potential chemotherapeutics drugs
To study the relationship between these model genes and the sensitivity of chemotherapeutic agents, we determined their expression level and IC50 of the drug across the NCI-60 cell lines. Next, we examined the association of 10 model gene expressions with IC50 for each drug type (Supplementary Table S2) and screened out the top 25 most significant model-gene-related drugs (Fig. 11A). Furthermore, several commonly used clinical drugs (including oxaliplatin, doxorubicin, cisplatin, docetaxel, gemcitabine, and paclitaxel) were analyzed and presented in Fig. 11B. It can be seen that high expression of TSPAN7 has a significant association with the resistance to Nelarabine, Pipobroman, Cytarabine, Asparaginase, Fludarabine, Dexamethasone Decadron, Thiotepa, and Bendamustine (p < 0.001); whereas, PDGFRA is sensitive indicators to Aliopurinol (p < 0.001). Similarly, the high expression of PYCR1 is significantly associated with Fluorouracil and Seliciclib resistance; on the other hand, its high expression exhibited a sensitivity to the drugs Chelerythrine and Depsipeptide. The high expression of SERPINB7 is a resistant indicator of Rebimastat and Simvastatin; conversely, it is a sensitive marker for Mithramycin and Homoharringtonine. AHCY is the resistance factor of Dabrafenib, Fludarabine, and XL-147. Besides, the high expression of PLOD1 is sensitive index to Dexrazoxane, Palbociclib, Oxaliplatin and Nelarabine; HSD17B1 is sensitive index to Docetaxel and Eribulin mesylate; EGR1 is sensitive to Palbociclib and AFP464; ATP6V1B1 is sensitive to Ixazomib citrate and 3-Bromopyruvate (acid). SCD is sensitive to Panobinostat, while it is insensitive to Salinomycin and SR16157.
Furthermore, through the ‘pRRophetic’ algorithm, we identified that the IC50 of four common chemotherapeutics drugs (Docetaxel, Vinblastine, Paclitaxel, and Doxorubicin) differed between the high- and low-risk groups of the model (p < 0.05). The high-risk group compared to the low-risk group had higher IC50 values for Docetaxel, Vinblastine, and Doxorubicin (P < 0.05, Wilcox test; Fig. 11C-E); on the contrary, Paclitaxel IC50 values in the high-risk group were higher than low-risk group (Fig. 11F). In addition, we performed a Spearman correlation analysis between the risk scores and IC50 of the above chemotherapy agents. These results verify the above analysis of Wilcox variance (P < 0.05; Fig. 11C-F).
Prediction value of the model for clinical immunotherapeutic efficacy among four BCa immunotherapy cohorts
The association analysis between model-gene expression and immunotherapy effectiveness was performed based on four independent real datasets (i.e., Imvigor210, GSE111636, GSE176307, and Truce01) (Fig. 12). We found that the expression levels of PYCR1 were higher and EGR1 was lower in the response group than in non-response across Imvigor210 and GSE176307 datasets. Moreover, TSRAN7 had higher expression levels in BCa patients with non-response according to Imvigor210 and Truce01; and AHCY was highly expressed in the non-response group from GSE111636 and Truce01. Additionally, ATP6V1B1 and SERPINB7 had low expression (P = 0.03 and 0.13, respectively); whereas PLOD1 and SCD were high expression (P = 0.067 and 0.2, respectively) in the response group from GSE111636 cohort.
In addition, in our ongoing single-arm phase II clinical study (TRUCE01), we also found that in BCa cases that responded to tislelizumab in combination with nab-paclitaxel, expression levels of PYCR1, AHCY, SCD were significantly decreased after treatment, and expression levels of TSPAN7, PDGFRA were significantly elevated after treatment. In non-responsive cases, the expression level of TSRAN7 and SERPINB7 was increased after treatment (Fig. 13). Of the other model genes, no other significant correlations were found (Supplementary Figure S3). Together, these results suggest that the expression of some model genes can help predict immunotherapy response.
Validation of model-related gene expression by HPA database and qRT-PCR
We then analyzed the protein expression levels of these ten prognostic genes using the Human Protein Atlas (HPA). Compared with normal bladder samples, the protein expression of PLOD1, ATP6V1B1, HSD17B1, SERPINB7, and PYCR1 in bladder cancer tissues was upregulated, while EGR1 in cancer tissues was down-regulated (Fig. 14A-J). Subsequently, mRNA expression levels of EGR1, PLOD1, and PYCR1 were detected in 10 pairs of BCa tissues and corresponding adjacent normal tissues. The qRT-PCR results of these genes were similar to the protein expression results described above (Fig. 14K-M).
In this study, we identified two molecular subtypes from metabolism-related prognostic DEGs through NMF algorithm based on TCGA datasets. And, we found that the two subgroups of patients had different prognosis (C1 was better than C2). Next, we constructed and verified the signature of 10 metabolism-related genes (including PLOD1, SERPINB7, TSPAN7, HSD17B1, PYCR1, PDGFRA, EGR1, AHCY, ATP6V1B1 and SCD) through LASSO, and multivariate Cox regression analyses using TCGA, GSE13507, GSE32894, GSE48075, GSE31684, and Imvigor210 datasets. Furthermore, we systematically investigated and analyzed the correlation of the model with clinical traits, immune infiltration, and sensitivity prediction of chemotherapeutic drugs or immunotherapy. Notely, most genes of the model are closely related to the metabolism of malignant tumors and cancer development.
In January 2021, Wang Z et al. discovered that hypoxia-induced PLOD1 overexpression promotes the malignant phenotype of glioblastoma through NF-κB signaling . In October 2022, Qiongjing Zeng et al. established a prediction model for cervical cancer with five signals including SERPINB7 indicating that SERPINB7 is highly expressed in cervical cancer patients . As early as September 2012, Wuttig D determined that TSPAN7 is a promising prognostic marker for clear cell renal cell carcinoma and indicated that patients with higher expression of the TSPAN7 gene or those with TSPAN7-positive blood vessels in both cores of tissue microarray studies had significantly longer DFS and tumor-specific survival (TSS) . Yu X et al. observed the downregulation of TSPAN7 in BCa tissue samples and cell lines and found that this downregulation was associated with relatively high tumor staging and tumor grade. Western blotting shows that overexpression of TSPAN7 activates Bax, and cleaves caspase-3 and PTEN, but inactivates Bcl-2, p-PI3K, and p-AKT, thereby inhibiting BCa cell growth via the PTEN/PI3K/AKT pathway . HSD17B1 is a steroid-metabolizing enzyme. Previous studies have shown that it is closely related to the occurrence and development of breast and endometrial cancer. PYCR1 is a key enzyme in proline synthesis. In August 2022, Kay EJ et al. discovered that cancer-associated fibroblasts need to synthesize proline by PYCR1 to deposit prototrophic extracellular matrix . Numerous studies have shown that PYCR1 is closely related to the prognosis of various cancers such as bladder cancer, pancreatic ductal adenocarcinoma, renal adenocarcinoma, gastric cancer, prostate cancer, hepatocellular carcinoma, colorectal cancer, etc. Platelet-derived growth subunit A (PDGFA) plays a key role in the development of glioblastoma (GBM). Guo L , Gao Z , and others also support the high correlation between PDGFRA and the prognosis of patients with bladder cancer. Tao T et al. showed that AHCY showed high potential as a prognostic factor for bladder cancer as a core gene of the co-expression network of lncRNA/mRNA and circRNA/mRNA constructed by WGCNA . The study of Bai S et al. showed that ATP6V1B1 showed excellent prognostic value as one of the nine genes in the hypoxia prognosis model for colorectal cancer.
The BCa samples in the TCGA dataset were then divided into high-risk and low-risk groups based on the median risk score. By analyzing the clustering and risk score, we could find that C2 had a higher risk score and C1 had a lower risk score. K-M survival analysis showed that the OS of the low-risk group was better than the low-risk group. According to the multivariate Cox regression results, the histogram of TNM stage risk score including age and sex was established. The ROC curve was shown to provide favorable predictive performance, and the calibration indicated good agreements between prediction and observation. The nomogram showed that the model was a valid and accurate tool. Our results show that this model can well distinguish BCa patients and predict prognosis, thus helping to formulate the best treatment plan based on risk score. Notably, the signature is well validated in several external BCa datasets (GSE13507, GSE32894, Imvigor210, GSE48075, and GSE31684). To further investigate the underlying mechanisms in the differences, we then performed pathway enrichment analyses for GSEA in each group.
GO analysis showed immunoglobulin receptor binding enrichment in the high-risk group. KEGG pathway enrichment analysis showed that ECM receptor interaction, focal adhesions, actin cytoskeletal regulation, and Leukocytes transendothelial migration were particularly significant enrichment in high-risk groups. Hallmark’s GSEA analysis shows epithelial-mesenchymal transition, angiogenesis, apical junction, myogenesis, hypoxia, inflammatory response, etc. gathering. Although immunoglobulins are about 90% similar at the amino acid level, each subclass has a unique way of antigen binding and immune complex formation. Triggering FcγR-expressing cells leads to a wide range of reactions, includes phagocytosis, antibody-dependent cell-mediated cytotoxicity, and complement activation. Tumor cell-derived IgG may hinder the cytotoxicity of antigen-dependent cells by binding antigens while lacking the ability to activate complement. These findings recommend tumor cell-derived IgG as a potential therapeutic target . CCL2 binds to the homologous receptor CCR2, a signaling pair that has been shown to have a variety of protumorigenic effects, from mediating tumor growth and angiogenesis to recruiting and usurping host stromal cells to support tumor progression . Wei Y et al. elucidate the molecular and cellular basis of KIR3DL3 inhibitory function, demonstrating that the KIR3DL3-HHLA2 pathway is a potential cancer immunotherapy target . However, immunoglobulin binding studies for bladder cancer need to be further studied.ECM plays an important role in tumor shedding, adhesions, degradation, motility, and hyperplasia . Pathology is characterized by abnormal neovascularization and diffuse infiltration of tumor cells. Their role in other cancers has been proven. The expression of ECM is upregulated in prostate cancer tissues , and it is involved in the process of tumor invasion and metastasis in gastric cancer . ECM in colorectal cancer promotes epithelial-mesenchymal transformation of cancer cells . Glioblastoma is the deadliest adult brain tumor. The interaction between the ECM and the glioblastoma microenvironment is important in this development . Focal adhesion (FA) is a group of macromolecular proteins that connect specialized actin ends with the extracellular matrix (ECM) and achieve cell migration, which is essential for the process of tumor metastasis . Studies have shown that lncRNA ITGB8-AS1 as ceRNA promotes colorectal cancer growth and migration through integrin-mediated plaque signaling . The hippo component YAP promotes adhesion plaque and tumor aggressiveness by transcriptionally activating THBS1/FAK signaling in breast cancer . Abnormal actin cytoskeletal dynamics have been implicated in a variety of diseases, including cancer . These may be potential targets for future treatments.
Immune cells in the tumor microenvironment (TME) play an important role in tumorigenesis, tumor progression, and metastasis. The tumor microenvironment consists of a heterogeneous population, including the cancer cells themselves, infiltrating immune cells, and stromal cells such as fibroblasts . Correlation analysis showed that risk scores were positively correlated with stromal scores and immune scores. Moreover, the risk scores in the TCGA and several GEO cohorts were positively correlated with PD-L1, POLE2, FEN1, MCM6, MSH6, MSH2, FAP, and TAGLIN.
PD-L1 expression in colorectal cancer tissues was negatively correlated with FOXP3 + cell density, suggesting that PD-L1-expressing cancer cells may affect regulatory T cells in the tumor microenvironment . Based on the results so far, we can define the expression of PD-L1 as a prognostic factor for immunotherapy and a predictor for pembrolizumab . FGFR3 disrupts PD-L1 via NEDD4 to control T cell-mediated immune surveillance of bladder cancer . Overexpression of PD-L1 and PD-1 on tumor cells and tumor-infiltrating lymphocytes, respectively, is associated with poor prognosis in certain human cancers . POLE2 promotes the malignant phenotype of glioblastoma by promoting AURKA-mediated stabilization of FOXM1 . None in bladder cancer. Upregulation and downregulation studies have shown that MCM6 regulates the cell cycle, proliferation, metastasis, immune response, and maintenance of DNA replication systems. MCM6 can also regulate downstream signals, such as MEK/ERK, to promote carcinogenesis . In many cancer cells, MCM6 expression is enhanced and can be used as a therapeutic target. Such as hepatocellular carcinoma , breast cancer  and gastric cancer , and so on. The RR for any cancer was 3.3 (95% CI 2.9 to 3.7) and 2.5 (95% CI 1.7 to 3.2) for path_MSH2 and path_MSH6 carriers, respectively. Older path_MSH2 carriers had a particularly high incidence of urinary tract and prostate cancer. Compared with path_MLH1 and path_MSH2 carriers, we found that path_MSH6 carriers had a lower risk of early-onset cancer and a lower risk of late-onset cancer in addition to an intermediate risk of urinary tract or prostate cancer . MSH6 and bladder cancer (OR, 5.63 [95% CI, 2.75–11.49]) . CircLIFR can interact with MSH2 to positively regulate CDDP sensitivity in bladder cancer through the MutSα/ATM-p73 axis. CircLIFR and MSH2 may be promising therapeutic targets in CDDP-resistant bladder cancer . There is evidence that MSH2 protein levels may contribute to the chemotherapy resistance observed in muscle-invasive bladder cancer. MSH2 has potential as a biomarker to predict response to platinum-based therapy . CAFs express the IL-6 cytokine, and its receptor IL-6R was found in RT4 bladder cancer cells. CM iCAF culture of RT4 bladder cancer cells resulted in significantly enhanced cell growth, migration, and invasion. Importantly, inhibition of CAFs-secreted IL-6 by neutralizing antibodies significantly reversed the IL-6-induced EMT phenotype, suggesting that this cytokine is essential for CAF-induced EMT in human bladder cancer progression. IL-6 expression is upregulated in invasive bladder cancer and correlates with the CAF marker ACTA2 . TAGLN is an antitumor gene in the human bladder. The expression level in normal bladder epithelial cells is higher than that in cancer cells .
In addition, immune interactions between tumors and TME play a key role in tumorigenesis and can be used as therapeutic targets for BCa . The composition and abundance of immune cells in TME influence tumor progression and the efficacy of immunotherapy . CD8 + T cells have been recognized as major effector cells of cell-mediated anti-tumor immunity, which kill tumor cells by releasing perforin . Interestingly, we found that CD8 + T cells were significantly higher in the high-risk group. We speculate that this may be related to its dysfunction in a depleted state . A recent study has indicated that the lack of innate inflammatory signaling in tumors leads to the inability to induce transcription factor-regulated functional effector differentiation, further impairing effector function and induced TOX expression and multiple other negative regulators of T cell signaling and function via persistent tumor antigen/TCR stimulation and/or other negative regulatory signals, ultimately leading to dysfunction of tumor-specific CD8 + T cells even before undergoing cell division . Studies have shown that immunosuppressive factors such as Tregs, and cancer-associated fibroblasts evade surveillance and clearance by the immune system through different mechanisms . In peripheral lymphoid organs, Treg cells are classified into resting and effector types, with effector Treg cells secreting IL-10 as an important characteristic [57, 58]. Treg cell physiology is dependent on the expression of GATA3 during inflammation . CD39 and CD73, which are highly expressed on the surface of Treg cells, increase intracellular AMP levels by breaking down ATP into adenosine and binding to adenosine receptor A2A (ADORA2A), activating the CREB pathway, thereby resulting in an anti-inflammatory milieu . Extensive research has been conducted on the role of cancer-associated fibroblasts (CAFs) in solid tumors, particularly in relation to their production of soluble factors such as IL-1α, IL-1β, CXCL1, CXCL12, G-CSF, and IL-6. These factors play a crucial role in recruiting monocytes and myeloid cells, thereby influencing the polarization of immune cells within the tumor microenvironment. Notably, the presence of CAFs has been found to induce the transformation of macrophages into an IL10-mediated tumor-promoting M2 phenotype . Neutrophils recruited to the site of inflammation promote cancer development primarily by increasing DNA damage, angiogenesis, and immunosuppression . Neutrophils can sustain tumor growth through different mechanisms, including inhibiting T cell activation and promoting the proliferation of genetically unstable tumor cells, angiogenesis, and metastasis. These mechanisms include the induction of genetic instability through the production of reactive oxygen species (ROS) and the release of microparticles containing microRNAs miR-23A and miR-155, etc. . Macrophages contribute to tumor progression at different stages, from initiation to distant metastasis formation. Evidence suggests that tumor cells induce macrophages to produce iconic acid and that acidosis in TME promotes immune escape . The research demonstrated by Nixon et al.  the immunosuppressive role of pTAMs, which is linked to their capacity to present tumor-associated antigens to CD8 + T cells and induce T cell exhaustion. A recent narrative review  proposed that mononuclear macrophages in tumor parenchyma and peritumoral regions of human hepatocellular carcinoma specimens exhibit heightened glycolytic activity, implying that TAMs’ glucose uptake facilitates tumor advancement. Several metabolites produced in glucose and lipid metabolism pathways, as well as those derived from amino acids, can also function as signaling molecules, promoting scavenging and anti-inflammatory functions in tumor-associated macrophages (TAMs). The pro-tumor activities exhibited by TAMs hinder patient responses to conventional chemotherapy, radiotherapy, and immunotherapy. Cancer cells can promote the formation and survival of endothelial cell tubes, at least in part through the PI3KAkt signaling pathway, thereby altering the microenvironment in favor of tumor growth . This is consistent with our findings of abundant CD8 + T cells, Neutrophils, Myeloid dendritic cells, M2 macrophages, Tregs, Endothelial cells, and Cancer-associated fibroblasts in patients in the BCa high-risk group. Therefore, we speculate that metabolic prognostic models may influence survival outcomes for BCa by altering ECM and immunosuppressive cells.
Despite this work were also similar to the previous study [68, 69], thereby lacking a certain amount of innovation, the model built by us was more streamlined and validated by multiple datasets, and its association analyses from multiple angles were performed. Although, a series of integrative analysis of multiple datasets from open databases (i.e., TCGA, GEO, TIDE, CellMiner, and TIMER) and our mRNA sequencing data (TRUCE01), as well as validation by immunohistochemistry and qRT-PCR were carried out, the main limitation of the research is that it lacks some functional experiments in vivo and in vitro to clarify the relevant molecular mechanisms of these modeled genes. Furthermore, further prospective studies are required to validate the clinical value of this metabolism-based molecular subtype and its signature.
Notably, our study showed unique immune landscapes, immune checkpoint gene expression, and immunotherapy responses between the high-risk and low-risk groups. In addition, we also calculated IC50 values to explore the chemotherapy drug’s sensitivity for BCa and screened candidate small molecule compounds. Furthermore, the TIDE and IPS algorithms were all used to predict the immunotherapy response of our model. Compared with the low-risk group, the high-risk BCa patients may be non-responsive and advanced for the immunotherapy. Additionally, the research also found that these model genes may act as a promising biomarker for predicting the efficacy of immunotherapy in BCa patients based on four independent real immunotherapy datasets, including IMvigor210, GSE111636, GSE176307, and our Truce01. These findings may provide suitable treatment options for patients with BCa. The proteins and mRNA expression of EGR1, PLOD1, and PYCR1 were also detected by the HPA database and qRT-PCR.
In summary, we developed and validated a new signature based on metabolism-related genes that may serve as a predictor for BCa prognosis, chemotherapy, or immunotherapy sensibility. Therefore, there are direct implications for guiding clinical oncology practice.
Availability of data and materials
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors. This study generated sequencing data that is available at TCGA (https://portal.gdc.cancer.gov) or NCBI GEO (https://www.ncbi.nlm.nih.gov/geo/).
The cancer genome atlas
Gene expression omnibus
Gene set enrichment analysis
Kyoto encyclopedia of genes and genomes
Receiver operating characteristic
Immune checkpoint inhibitor
Programmed cell death-ligand 1
Programme death 1
Overall survival curves
Progression free survival
Disease free survival
Disease specific survival
Least absolute shrinkage and selection operator
Area under the ROC curve
Decision curve analysis
The tumor immune dysfunction and exclusion
The tumor inflammation signature
Differentially expressed genes
Semi-maximum inhibitory concentrations
Richters A, Aben K, Kiemeney L. The global burden of urinary bladder cancer: an update. World J Urol. 2020;38(8):1895–904.
DeBerardinis RJ, Chandel NS. Fundamentals of cancer metabolism. Sci Adv. 2016;2(5):e1600200.
Gyamfi J, Kim J, Choi J. Cancer as a metabolic disorder. Int J Mol Sci. 2022;23(3):1155.
Zuzčák M, Trnka J. Cellular metabolism in pancreatic cancer as a tool for prognosis and treatment (Review). Int J Oncol. 2022;61(2):93.
Pardo JC, Ruiz de Porras V, Gil J, Font A, Puig-Domingo M, Jordà M. Lipid metabolism and epigenetics crosstalk in prostate cancer. Nutrients. 2022;14(4):851.
Yang C, Huang X, Liu Z, Qin W, Wang C. Metabolism-associated molecular classification of hepatocellular carcinoma. Mol Oncol. 2020;14(4):896–913.
Zuo D, Li C, Liu T, Yue M, Zhang J, Ning G. Construction and validation of a metabolic risk model predicting prognosis of colon cancer. Sci Rep. 2021;11(1):6837.
Jiang S, Ren X, Liu S, Lu Z, Xu A, Qin C, et al. Integrated analysis of the prognosis-associated RNA-binding protein genes and candidate drugs in renal papillary cell carcinoma. Front Genet. 2021;12: 627508.
Han Q, Zhang X, Ren X, Hang Z, Yin Y, Wang Z, et al. Biological characteristics and predictive model of biopsy-proven acute rejection (BPAR) after kidney transplantation: evidences of multi-omics analysis. Front Genet. 2022;13: 844709.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587–587D592.
Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417–25.
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550–8.
Wang Z, Shi Y, Ying C, Jiang Y, Hu J. Hypoxia-induced PLOD1 overexpression contributes to the malignant phenotype of glioblastoma via NF-κB signaling. Oncogene. 2021;40(8):1458–75.
Zeng Q, Jiang H, Lu F, Fu M, Bi Y, Zhou Z, et al. Prediction of the immunological and prognostic value of five signatures related to fatty acid metabolism in patients with cervical cancer. Front Oncol. 2022;12:1003222.
Wuttig D, Zastrow S, Füssel S, Toma MI, Meinhardt M, Kalman K, et al. CD31, EDNRB and TSPAN7 are promising prognostic markers in clear-cell renal cell carcinoma revealed by genome-wide expression analyses of primary tumors and metastases. Int J Cancer. 2012;131(5):E693-704.
Yu X, Li S, Pang M, Du Y, Xu T, Bai T, et al. TSPAN7 exerts anti-tumor effects in bladder cancer through the PTEN/PI3K/AKT pathway. Front Oncol. 2020;10: 613869.
Kay EJ, Paterson K, Riera-Domingo C, Sumpton D, Däbritz J, Tardito S, et al. Cancer-associated fibroblasts require proline synthesis by PYCR1 for the deposition of pro-tumorigenic extracellular matrix. Nat Metab. 2022;4(6):693–710.
Guo L, Wu Q, Ma Z, Yuan M, Zhao S. Identification of immune-related genes that predict prognosis and risk of bladder cancer: bioinformatics analysis of TCGA database. Aging (Albany NY). 2021;13(15):19352–74.
Gao Z, Chen C, Gu P, Chen J, Liu X, Shen J. The tumor microenvironment and prognostic role of autophagy- and immune-related genes in bladder cancer. Cancer Biomark. 2022;35(3):293–303.
Tao T, Yuan S, Liu J, Shi D, Peng M, Li C, et al. Cancer stem cell-specific expression profiles reveal emerging bladder cancer biomarkers and identify circRNA_103809 as an important regulator in bladder cancer. Aging (Albany NY). 2020;12(4):3354–70.
Kdimati S, Mullins CS, Linnebacher M. Cancer-cell-derived IgG and its potential role in tumor development. Int J Mol Sci. 2021;22(21):11597.
Xu M, Wang Y, Xia R, Wei Y, Wei X. Role of the CCL2-CCR2 signalling axis in cancer: mechanisms and therapeutic targeting. Cell Prolif. 2021;54(10):e13115.
Wei Y, Ren X, Galbo PM Jr, Moerdler S, Wang H, Sica RA, et al. KIR3DL3-HHLA2 is a human immunosuppressive pathway and a therapeutic target. Sci Immunol. 2021;6(61):eabf9792.
Verschueren E, Husain B, Yuen K, Sun Y, Paduchuri S, Senbabaoglu Y, et al. The immunoglobulin superfamily receptome defines cancer-relevant networks associated with clinical outcome. Cell. 2020;182(2):329-44.e19.
Andersen MK, Rise K, Giskeødegård GF, Richardsen E, Bertilsson H, Størkersen Ø, et al. Integrative metabolic and transcriptomic profiling of prostate cancer tissue containing reactive stroma. Sci Rep. 2018;8(1):14269.
Yan P, He Y, Xie K, Kong S, Zhao W. In silico analyses for potential key genes associated with gastric cancer. PeerJ. 2018;6:e6092.
Rahbari NN, Kedrin D, Incio J, Liu H, Ho WW, Nia HT, et al. Anti-VEGF therapy induces ECM remodeling and mechanical barriers to therapy in colorectal cancer liver metastases. Sci Transl Med. 2016;8(360):360ra135.
Cui X, Morales RT, Qian W, Wang H, Gagner JP, Dolgalev I, et al. Hacking macrophage-associated immunosuppression for regulating glioblastoma angiogenesis. Biomaterials. 2018;161:164–78.
Legerstee K, Houtsmuller AB. A layered view on focal adhesions. Biology (Basel). 2021;10(11):1189.
Lin X, Zhuang S, Chen X, Du J, Zhong L, Ding J, et al. lncRNA ITGB8-AS1 functions as a ceRNA to promote colorectal cancer growth and migration through integrin-mediated focal adhesion signaling. Mol Ther. 2022;30(2):688–702.
Shen J, Cao B, Wang Y, Ma C, Zeng Z, Liu L, et al. Hippo component YAP promotes focal adhesion and tumour aggressiveness via transcriptionally activating THBS1/FAK signalling in breast cancer. J Exp Clin Cancer Res. 2018;37(1):175.
Lappalainen P, Kotila T, Jégou A, Romet-Lemonne G. Biochemical and mechanical regulation of actin dynamics. Nat Rev Mol Cell Biol. 2022;23(12):836–52.
Lei X, Lei Y, Li JK, Du WX, Li RG, Yang J, et al. Immune cells within the tumor microenvironment: Biological functions and roles in cancer immunotherapy. Cancer Lett. 2020;470:126–33.
Masugi Y, Nishihara R, Yang J, Mima K, da Silva A, Shi Y, et al. Tumour CD274 (PD-L1) expression and T cells in colorectal cancer. Gut. 2017;66(8):1463–73.
Fabrizio FP, Trombetta D, Rossi A, Sparaneo A, Castellana S, Muscarella LA. Gene code CD274/PD-L1: from molecular basis toward cancer immunotherapy. Ther Adv Med Oncol. 2018;10:1758835918815598.
Jing W, Wang G, Cui Z, Xiong G, Jiang X, Li Y, et al. FGFR3 destabilizes PD-L1 via NEDD4 to control T-cell-mediated bladder cancer immune surveillance. Cancer Res. 2022;82(1):114–29.
Ohaegbulam KC, Assal A, Lazar-Molnar E, Yao Y, Zang X. Human cancer immunotherapy with antibodies to the PD-1 and PD-L1 pathway. Trends Mol Med. 2015;21(1):24–33.
Zhang P, Chen X, Zhang L, Cao D, Chen Y, Guo Z, et al. POLE2 facilitates the malignant phenotypes of glioblastoma through promoting AURKA-mediated stabilization of FOXM1. Cell Death Dis. 2022;13(1):61.
Zeng T, Guan Y, Li YK, Wu Q, Tang XJ, Zeng X, et al. The DNA replication regulator MCM6: an emerging cancer biomarker and target. Clin Chim Acta. 2021;517:92–8.
Blanc V, Riordan JD, Soleymanjahi S, Nadeau JH, Nalbantoglu I, Xie Y, et al. Apobec1 complementation factor overexpression promotes hepatic steatosis, fibrosis, and hepatocellular cancer. J Clin Invest. 2021;131(1):e138699.
Issac M, Yousef E, Tahir MR, Gaboury LA. MCM2, MCM4, and MCM6 in breast cancer: clinical utility in diagnosis and prognosis. Neoplasia. 2019;21(10):1015–35.
Wang Y, Chen H, Liu W, Yan H, Zhang Y, Cheung A, et al. MCM6 is a critical transcriptional target of YAP to promote gastric tumorigenesis and serves as a therapeutic target. Theranostics. 2022;12(15):6509–26.
Møller P, Seppälä TT, Bernstein I, Holinski-Feder E, Sala P, Gareth Evans D, et al. Cancer risk and survival in path_MMR carriers by gene and gender up to 75 years of age: a report from the Prospective Lynch Syndrome Database. Gut. 2018;67(7):1306–16.
Zeng C, Bastarache LA, Tao R, Venner E, Hebbring S, Andujar JD, et al. Association of pathogenic variants in hereditary cancer genes with multiple diseases. JAMA Oncol. 2022;8(6):835–44.
Zhang H, Xiao X, Wei W, Huang C, Wang M, Wang L, et al. CircLIFR synergizes with MSH2 to attenuate chemoresistance via MutSα/ATM-p73 axis in bladder cancer. Mol Cancer. 2021;20(1):70.
Goodspeed A, Jean A, Costello JC. A whole-genome CRISPR screen identifies a role of MSH2 in cisplatin-mediated cell death in muscle-invasive bladder cancer. Eur Urol. 2019;75(2):242–50.
Goulet CR, Champagne A, Bernard G, Vandal D, Chabaud S, Pouliot F, et al. Cancer-associated fibroblasts induce epithelial-mesenchymal transition of bladder cancer cells through paracrine IL-6 signalling. BMC Cancer. 2019;19(1):137.
Tsui KH, Lin YH, Chang KS, Hou CP, Chen PJ, Feng TH, et al. Transgelin, a p53 and PTEN-upregulated gene, inhibits the cell proliferation and invasion of human bladder carcinoma cells in vitro and in vivo. Int J Mol Sci. 2019;20(19):4946.
Kamoun A, de Reyniès A, Allory Y, Sjödahl G, Robertson AG, Seiler R, et al. A consensus molecular classification of muscle-invasive bladder cancer. Eur Urol. 2020;77(4):420–33.
Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–509W514.
T G S. Innate and adaptive immune cells in tumor microenvironment. Gulf J Oncolog. 2021;1(35):77–81.
Kersten K, Hu KH, Combes AJ, Samad B, Harwin T, Ray A, et al. Spatiotemporal co-dependency between macrophages and exhausted CD8+ T cells in cancer. Cancer Cell. 2022;40(6):624-38.e9.
Rudloff MW, Zumbo P, Favret NR, Roetman JJ, Detrés Román CR, Erwin MM, et al. Hallmarks of CD8+ T cell dysfunction are established within hours of tumor antigen encounter before cell division. Nat Immunol. 2023;24(9):1527–39.
Hanahan D, Coussens LM. Accessories to the crime: functions of cells recruited to the tumor microenvironment. Cancer Cell. 2012;21(3):309–22.
Cretney E, Kallies A, Nutt SL. Differentiation and function of Foxp3(+) effector regulatory T cells. Trends Immunol. 2013;34(2):74–80.
Smigiel KS, Richards E, Srivastava S, Thomas KR, Dudda JC, Klonowski KD, et al. CCR7 provides localized access to IL-2 and defines homeostatically distinct regulatory T cell subsets. J Exp Med. 2014;211(1):121–36.
Wohlfert EA, Grainger JR, Bouladoux N, Konkel JE, Oldenhove G, Ribeiro CH, et al. GATA3 controls Foxp3+ regulatory T cell fate during inflammation in mice. J Clin Invest. 2011;121(11):4503–15.
Antonioli L, Pacher P, Vizi ES, Haskó G. CD39 and CD73 in immunity and inflammation. Trends Mol Med. 2013;19(6):355–67.
Caligiuri G, Tuveson DA. Activated fibroblasts in cancer: Perspectives and challenges. Cancer Cell. 2023;41(3):434–49.
Xiong S, Dong L, Cheng L. Neutrophils in cancer carcinogenesis and metastasis. J Hematol Oncol. 2021;14(1):173.
Hedrick CC, Malanchi I. Neutrophils in cancer: heterogeneous and multifaceted. Nat Rev Immunol. 2022;22(3):173–87.
Locati M, Curtale G, Mantovani A. Diversity, mechanisms, and significance of macrophage plasticity. Annu Rev Pathol. 2020;15:123–47.
Nixon BG, Kuo F, Ji L, Liu M, Capistrano K, Do M, et al. Tumor-associated macrophages expressing the transcription factor IRF8 promote T cell exhaustion in cancer. Immunity. 2022;55(11):2044-58.e5.
Zhang X, Ji L, Li MO. Control of tumor-associated macrophage responses by nutrient acquisition and metabolism. Immunity. 2023;56(1):14–31.
Cheng HW, Chen YF, Wong JM, Weng CW, Chen HY, Yu SL, et al. Cancer cells increase endothelial cell tube formation and survival by activating the PI3K/Akt signalling pathway. J Exp Clin Cancer Res. 2017;36(1):27.
Qiu T, Chen Y, Meng L, Xu T, Zhang H. Identification of a metabolism-related gene signature predicting overall survival for bladder cancer. Genomics. 2022;114(4): 110402.
Shen C, Liu J, Wang L, Liang Z, Niu H, Wang Y. Identification of metabolism-associated genes and construction of a prognostic signature in bladder cancer. Cancer Cell Int. 2020;20(1):538.
We sincerely thank all participants in the study.
This study was supported by Tianjin Municipal Health Industry Key Project (No. TJWJ2022XK014), Scientific Research Project of Tianjin Municipal Education Commission (No. 2022ZD069 and No. 2020KJ169), The Natural Science Youth Cultivated Foundation of Second Hospital of Tianjin Medical University (No. 2022ydey15 and No. 2017ydey19), The Talents Cultivated Project of Department of Urology, the Second Hospital of Tianjin Medical University (No. MNRC202313), and Tianjin Key Medical Discipline (Specialty) Construction Project (No. TJYXZDXK-032A).
Ethics approval and consent to participate
This research has been conducted by the Declaration of Helsinki and has been approved by the ethics committee of the 2nd Affiliated Hospital of Tianjin Medical University (Ethics code: KY2021K003). All participants in the study have signed the informed consent.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
To achieve a preferable clustering performance, we set the k cluster range from 2 to 10 using the “NMF” R package. Supplementary Figure S2. Kaplan–Meier analyses of OS in different clinicopathological subgroups stratified by (A) age, (B) gender, (C) stage, and (D) stage_T. Supplementary Figure S3. Based on (A) Imvigor210, (B) GSE111636, (C) GSE176307, and (D) our Truce01 cohorts, the partial genes expression of the model were not statistically different between response and non-response group. (E-H) The four model genes (including HSD17B1, EGR1, PLOD1, and ATP6V1B1) showed no difference in expression before and after immunotherapy combined with nab-paclitaxel. Res., Response; Non-Res., Non-Response.
Based on (A) TCGA, (B) GSE13507, and (C) GSE32894 datasets, GO, KEGG, or HALLMARK enrichment analysis were compared between the high-risk and low-risk cohorts using the GSEA.
Drug sensitivity analysis of metabolism-associated model genes based on the CellMiner database.
About this article
Cite this article
Shen, C., Bi, Y., Chai, W. et al. Construction and validation of a metabolism-associated gene signature for predicting the prognosis, immune landscape, and drug sensitivity in bladder cancer. BMC Med Genomics 16, 264 (2023). https://doi.org/10.1186/s12920-023-01678-6