Prognostic model development and molecular subtypes identification in bladder urothelial cancer by oxidative stress signatures

Data download and extraction

From the TCGA-BLCA database, we acquired the transcriptome profiles of 414 BLCA samples and 19 samples of normal bladder tissue. After excluding bladder cancer cases with no clinical information, we ultimately retrieved transcriptome profiling for 406 samples with complete information. The TCGA database provided clinical annotations including clinicopathological characteristics and survival rate. To validate externally, we employed information from the GSE31684 group, acquired from the GEO repository, comprising microarray data from 93 instances of bladder cancer that were examined using the Affymetrix Human Genome U133 Plus 2.0 Array (Platform GPL570). When collecting OS-related genes from the Genecards database and published literature (Supplementary Table 1), a sum of 275 potential genes was incorporated.

Selection of differentially expressed genes (DEGs) related to OS

To identify differentially expressed genes (DEGs) associated with oxidative stress (OS), we employed the ‘limma’ R package to analyze the gene expression levels in bladder cancer in comparison to normal bladder tissues. To identify significant DEGs, the criteria included an FDR below 0.05 and an absolute log2 FC exceeding 1. In this context, a log2FC above 1 indicated gene upregulation, whereas a log2FC below -1 indicated gene downregulation in tumor samples.

Molecular subtype identification using non-negative matrix factorization (NMF) algorithm

To perform molecular subtype identification, we employed the NMF algorithm using the ‘NMF’ R package. The decision to use NMF was motivated by its ability to identify and group the inherent characteristics present in tumor samples, which is crucial for precise molecular subtyping. Additionally, this method aided in the retrieval of biologically significant correlation coefficients, offering a more profound understanding of the fundamental biology of tumors. To ensure the reliability and accuracy of the subtype classification, the selection of the best k value, which ranged from 2 to 10, was guided by finding a balance between stability and clustering performance. This involved evaluating the covariance and the sum of squared residuals (RSS).

Definition and comparison of tumor-infiltrating immune cells

In BLCA samples, we employed the MCP-counter to examine immune cells that have infiltrated the tumor microenvironment. The process included normalizing gene expression data, utilizing an MCP-counter to measure immune and stromal cells, and comparing these scores to detect important tumor-infiltrating cells.

Construction and validation of OS-related genes prognosis features

To identify DEGs with potential prognostic significance, we utilized Cox regression analyses in R with the assistance of the ‘survival’ and ‘survminer’ packages [15]. By utilizing this method, we were able to identify possible genes that could be examined in more detail. Afterward, we utilized the Lasso technique, employing the ‘glmnet’ software in R, to build a predictive model. This model integrated these candidate genes, focusing on their potential predictive power in relation to patient outcomes [16]. The risk score was calculated utilizing the equation [17].

$$Riskscore={\displaystyle \sum _{i=1}^n(coef\text{}i*Expi)}$$

In our study, we utilized a particular mathematical formula to determine the risk score. This formula involves multiplying the expression levels of each chosen candidate gene by their respective regression coefficients obtained from the multivariate Cox regression analysis. Adding up these products provides the risk score for every individual. Patients were categorized into high-risk and low-risk groups based on these scores, using the median risk score as the dividing factor.

To validate our prognostic model, we trained it using data from TCGA and then tested it on data from the GEO database. Confirmation of the model’s stability and effectiveness was achieved by analyzing the Receiver Operating Characteristic (ROC) curve. In order to establish the risk score as a crucial element in forecasting prognosis, we conducted both univariate and multivariate Cox regression analyses. Furthermore, the R package ‘heatmap’ was employed to visually depict and examine the correlation between the risk score and different clinical factors.

Gene set enrichment analysis (GSEA)

In order to explore the functional annotations linked to high-risk and low-risk groups in our investigation, we employed the Gene Set Enrichment Analysis (GSEA) program, specifically version 4.1. The examination was performed utilizing the c2.cp.kegg.v7.4.symbols gene set compilation [18]. Results were deemed statistically significant if the p-value was below 0.05. To ensure effective visualization and interpretation, we have chosen and presented the top eight significant results obtained from this analysis.

Establishment and verification of the nomogram

The nomogram was developed using the ‘rms’ package in R, and its visualization was performed with the ‘regplot’ package [19]. The model integrated the risk scores from the nomogram and the clinical characteristics of the patients. The nomogram offers a numerical instrument to forecast the overall lifespan within 1, 3, and 5-year timeframes. Calibration curves were created and examined to evaluate the precision and effectiveness of the nomogram. Verifying the effectiveness and applicability of the nomogram in clinical prognosis was of utmost importance in this process.

Correlation between TMB and risk scores

Somatic mutation data was obtained from the TCGA-BLCA cohort and then assessed using the R program ‘maftools’ [20] for somatic non-synonymous point mutations. TMB refers to the aggregate count of genetic coding mistakes, base substitutions, gene insertions, or deletions identified per million base pairs.

Correlation of risk score with tumor immune infiltration microenvironment (TIME) characterization

To evaluate the presence of immune cells in the tumor microenvironment [21], we employed a total of 7 techniques: XCELL, TIMER, QUANTISEQ, MCP-counter, EPIC, CIBERSORT, and CIBERSORT-ABS. To visualize the aforementioned results, the package ‘ggpubr’ was utilized. Every one of these techniques provides a distinct strategy for evaluating the presence of immune cells within the tumor microenvironment. XCELL and TIMER offer information on the prevalence and positioning of immune cells, QUANTISEQ and MCP-counter concentrate on measuring immune cell populations, EPIC expands these evaluations to encompass a wider array of cell types, whereas CIBERSORT and CIBERSORT-ABS employ gene expression profiles to deduce the makeup of cell mixtures in tumor samples. By adopting this all-encompassing method, a meticulous assessment of TIME is guaranteed, thereby augmenting comprehension of its correlation with the risk score in BLCA.

Gene set variation analysis (GSVA)

We obtained two gene sets, the c7.immunesigdb HALLMARK and c2.cp.kegg.v7.4, from the Molecular Signatures Database for our analysis. Furthermore, we have incorporated a compilation of 47 genes renowned for their correlation with immune checkpoint blockade. By utilizing these sets of genes, we conducted GSVA in a standardized manner for every sample in order to measure the comparative level of activity in different pathways, with particular emphasis on immunological characteristics and immune regulatory points. This analysis allowed us to ascertain the functional relevance of these pathways and checkpoints in the context of bladder cancer.

Anticancer drug sensitivity prediction

In clinical practice, acknowledging chemotherapy and immunotherapy as the main treatment approaches for bladder cancer, our objective was to compare the half-maximal inhibitory concentration (IC50) of medications on BLCA samples between groups categorized as high-risk and low-risk. We used the “pRRophetic” package, which predicts drug sensitivity based on gene expression profiles from TCGA and cell line expression profiles from Genomics of Cancer Drug Sensitivity. By utilizing this method, we were able to evaluate the possible effectiveness of these therapies in various categories of risk, which aids in making personalized decisions regarding treatment. Immune checkpoint inhibitors are a form of therapy used to treat advanced bladder cancer. The Immunopheno Score (IPS) was utilized to evaluate the immunogenicity of tumors in BLCA populations categorized as high-risk and low-risk. The Immune Profiling Score (IPS) was computed using genes linked to immune checkpoint inhibition, utilizing gene information obtained from the Genecards repository.

Statistical analysis

In our study, we utilized the Kruskal-Wallis test to compare more than two groups and the Wilcoxon test to compare two groups. To assess differences in survival curves, we utilized the Kaplan-Meier log-rank test, a highly effective method for evaluating survival data. The chi-square analysis played a crucial role in investigating the connections between risk group categorizations and the burden of somatic mutations. Furthermore, we employed Spearman’s rank correlation analysis to compute correlation coefficients, offering valuable insights into the associations among different study variables. In order to conduct a more in-depth analysis of the data, specifically regarding the makeup of immune cells, the CIBERSORT algorithm was utilized, with outcomes deemed statistically significant if the p-value was below 0.05. R software (version 4.1.1) was utilized for all statistical analyses, incorporating specific R packages like ‘survival’ for Kaplan-Meier analysis and ‘ggplot2’ for data visualization [22]. A significance level of 0.05 was used to determine statistical significance based on two-sided p-values.

Data availability

All data used in the study were from the publicly available The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). Informed consent forms are not required for patient data extracted from public databases.

Identification of DEGs between normal bladder and BLCA tissues

Figure 1 displays the flowchart of analyses. Transcriptomic data was acquired from the TCGA-BLCA database, consisting of 414 BLCA samples and 19 samples of adjacent non-cancerous tissue. Finally, the transcriptome data of 406 BLCA patients were obtained after excluding cases with clinical information. Our research primarily centered around a collection of 275 genes associated with OS, where we extensively examined and compared their expression patterns. In order to guarantee the integrity and precision of the data, our initial focus was on resolving batch effects in the TCGA dataset, which was then followed by the normalization of the data. Table 1 provides a detailed account of the 43 differentially expressed genes (DEGs) that were identified between BLCA and normal bladder tissues through this procedure. Illustrated are the clinic features of BLCA individuals obtained from the TCGA and GEO databases. It is worth mentioning that most of the individuals in both groups are over the age of 65 and primarily male. Additionally, the data reveals a significant percentage of patients in advanced phases (III-IV) and with tumors of high grade, specifically within the TCGA-BLCA group. The differentially expressed genes (DEGs), which exhibited either increased or decreased expression, were visually depicted using volcano plots (as shown in Figure 2A). This graphical representation offered a straightforward and succinct visualization of their expression patterns.

Figure 1. Graphical abstract of construction of a prognostic index associated with oxidative stress signatures in bladder urothelial cancer.

Figure 2. (A) The Volcano plot of differentially expressed genes from the BLCA patients in the TCGA database. (B) Non-negative matrix factorization (NMF) algorithm-based Consensus map clustered. (C) The stability of the cluster obtained from NMF was determined by the cophenetic correlation coefficient. (D) The clustering performance of the model was measured by RSS. (E, F) The statistical differences between C1 and C2 in overall survival (OS) and progression-free survival (PFS). (G) The mapping relationship of C1 and C2 to molecular subtypes was exhibited in the Alluvial plot.

NMF algorithm-based molecular subtypes of OS-related genes

By employing covariance and RSS, we identified two clusters, namely C1 and C2, where k=2 was determined to be the optimal clustering parameter considering stability and performance, as depicted in Figure 2B–2D. Figure 2E, 2F depicted that the C2 cluster exhibited enhanced overall survival (OS) and progression-free survival (PFS) in comparison to C1, as indicated by the Kaplan-Meier analysis. Furthermore, Figure 2G showcases the allocation of DEGs across the traditional immune subtypes C1-C6, encompassing wound repair, IFN-gamma predominant, inflammatory, lymphocyte-depleted, immunologically silent, and TGF-beta predominant categories. The different subcategories are associated with particular types of immune cells and indicators, which have an impact on the patient’s prognosis and how they respond to treatment. Figure 3 shows that there were noticeable differences in immune and stromal scores between C1 and C2, highlighting the grouping of OS-related genes and their connection to various immune cell clusters.

Figure 3. The statistical difference in the immune cells (A–H) and stromal cells (I, J) from the tumor microenvironment (TME) between C1 and C2.

Developing and validating a prognostic prediction model for OS-related genes

We conducted a univariate Cox regression analysis using 43 differentially expressed genes (DEGs) (see Supplementary Table 2 and Figure 4A). In order to avoid overfitting, we applied Lasso regression to these candidate genes (Figure 4B). Cross-validation was also used to determine the ideal value of the penalty parameter, and 13 genes associated with OS were identified (Figure 4C). Following the implementation of multivariate Cox regression, we ultimately recognized four genes (RBPMS, CRYAB, P4HB, PDGFRA) that are crucial for constructing the model. CRYAB, P4HB, and PDGFRA were regarded as unfavorable prognostic markers (HR > 1), while RBPMS was regarded as a favorable prognostic marker (HR < 1) (Supplementary Table 3).

Figure 4. (A) The survival showed as Hazard ratio was determined by univariate Cox regression analysis in Forest plot. (B) Ten-fold cross-validation for tuning parameter selection in lasso regression. Vertical lines are drawn from the best data according to the minimum criterion and 1 standard error criterion. The vertical lines on the left represent the final 4 genes identified. (C) LASSO coefficient profiles of 13 OS-related genes. The 10-fold cross-validation value is marked by a vertical line. (D) The statistical difference exists in the overall survival between high-risk and low-risk groups as shown in the Kaplan-Meier curve. (E–H) The different overall survival between low and high CRYAB, P4HB, PDGFRA, and RBPMS expression was shown in the Kaplan-Meier curve respectively.

The risk score for each sample was determined using the formula: Risk Score = -(0.2379 × RBPMS) + (0.1077 × CRYAB) + (0.4696 × P4HB) + (0.2422 × PDGFRA).

The tumor samples were categorized into high-risk and low-risk groups based on the median risk score (1.02) serving as the threshold. According to the findings from the Kapan Meier curve (Figure 4D), individuals categorized as high-risk experienced shorter overall survival durations compared to those classified as low-risk. A notable distinction was observed in overall survival (OS) between the two groups of patients based on the varying levels of expression of the four genes (Figure 4E–4H).

Furthermore, we demonstrated RS in the TCGA group (Figure 5A), the survival condition (Figure 5C), and the distribution of polygenic model RS among BLCA individuals (Figure 5E). Subsequently, we validated these findings in the external GEO group (Figure 5B, 5D, 5F).

Figure 5. Validation the prognostic predictive value of risk scores in the TCGA cohort (A) and GEO cohort (B). Survival status and duration of BLCA patients in the TCGA cohort (C) and GEO cohort (D). The polygenic model risk score distribution in the TCGA cohort (E) and GEO cohort (F).

Construction of prognostic nomogram

In Figure 6A, the analysis of the Receiver Operating Characteristic (ROC) curve for our predictive model revealed AUC values of 0.656, 0.668, and 0.696 for the 1-year, 3-year, and 5-year rates of overall survival, respectively. In order to confirm the importance of the risk score as a significant predictor among different clinicopathological features, we included variables such as sex, tumor grade, and TNM stages in the AUC analysis for 1-year (Figure 6B), 3-year (Figure 6C), and 5-year (Figure 6D) overall survival. The combined factors resulted in lower AUC values compared to the risk score alone.

Figure 6. Risk signatures prognosis efficiency validation. (A) The predictive value of prognostic features was validated in the ROC analysis. The area under the curve (AUC) of the risk score to predict the overall survival at 1, 3, and 5 years and (B–D) corresponding other clinicopathological factors. (E) Univariate and (F) Multivariate Cox regression results for overall survival. (G) The predicted survival of BLCA patients by Nomogram, and (H) 1-, 3-, and 5-year nomogram calibration curves.

In terms of clinical aspects, prognosis is influenced by factors including age, gender, stage, and grade. Significant disparities were observed in the overall survival outcomes between univariate and multivariate Cox regression, as depicted in Figure 6E, 6F. In order to assess the accuracy of the model’s predictions, we created a predictive nomogram that included the risk score and clinicopathological characteristics (Figure 6G). As an example, a male who is 73 years old and diagnosed with BLCA, classified as T3N0M0 and placed in the low-risk category, as per the model, is estimated to have survival probabilities of 89.9%, 66.5%, and 55.6% at 1, 3, and 5 years respectively. This estimation is based on a total nomogram score of 414. However, the calibration plot suggested a potential underestimation of the 3-year survival percentage (Figure 6H).

Functional analysis of OS-related genes

All samples were classified into low or high-expression groups based on the median expression of the hub genes. Following that, an analysis called Gene Set Enrichment Analysis (GSEA) was performed, which unveiled the presence of functional enrichment within these genes. According to the KEGG pathway analysis, an increase in CRYAB expression was linked to particular signaling pathways, particularly the interaction between cytokines and cytokine receptors, as well as the chemokine signaling pathways (Figure 7A). Likewise, heightened P4HB expression was found to be associated with the activation of the cytosolic DNA sensing pathway, along with other pathways (Figure 7B). Two distinct pathways, namely the chemokine signaling pathway and the cytokine-receptor interaction system (Figure 7C), were associated with elevated expression levels of PDGFRA. In addition, increased RBPMS expression was observed to be linked to the activation of pathways including cytochrome P450-mediated xenobiotic metabolism and retinol metabolism (Figure 7D). The results highlight the intricate interaction between the expression of these central genes and their participation in different crucial signaling pathways within the cellular environment.

Figure 7. GSEA for samples with high and low expression of 4 central genes. (A) The enriched gene set obtained from KEGG for samples with high CRYAB expression. (B) The enriched gene set obtained from KEGG for samples with high P4HB expression. (C) The enriched gene set obtained from KEGG for samples with high PDGFRA expression. (D) The enriched gene set obtained from KEGG for samples with high RBPMS expression.

Association of risk characteristics with clinicopathological variables

Using clinicopathological features, we created a plot (Figure 8A) to investigate the relationship between risk and clinicopathological variables. Gender (Figure 8B), tumor grade (Figure 8C), tumor stage (Figure 8D), stage T (Figure 8E), stage N (Figure 8F), and stage M (Figure 8G) exhibited variations that aligned with clinical results. The results of these findings, along with the outcomes of univariable and multivariable regression analysis, suggest that our risk score serves as a reliable prognostic indicator, regardless of other clinical factors.

Figure 8. The clinicopathologic feature of the prognostic risk characteristics. (A) The clinicopathologic feature and corresponding risk scores were shown in the heatmap. The clinicopathologic feature incidence (B) Gender, (C) Tumor Grade, (D) Tumor stage, (E) Stage T, (F) Stage N, and (G) Stage M in the low and high-risk groups.

Association of risk signature with TMB

In order to investigate the correlation between risk score and gene mutations, our initial analysis focused on examining the prevalence of gene alterations within both high-risk and low-risk groups. Next, our attention shifted towards visualizing the waterfalls, which represent the top 20 genes exhibiting the highest occurrence of somatic mutations in each group (Figure 9A, 9B). Notably, the mutation profiles of significantly mutated genes (SMGs) indicated a higher prevalence of TP53 mutations in the high-risk group (54% as opposed to 39%), whereas FGFR3 mutations were found to be more frequent in the low-risk group (22% compared to 5%). This discrepancy is attributed to TP53’s known association with an increased cancer risk.

Figure 9. Correlation between risk score and TMB. High risk score. High-risk score (A) and low-risk score (B) were shown in the established oncoPrint. (C) High TMB and low TMB groups were shown in the Kaplan-Meier curve. (D) Kaplan-Meier curve stratification based on TMB and risk signature.

The results of the Kaplan-Meier analysis showed a positive correlation between elevated TMB levels and extended overall survival (OS) (Figure 9C). The TMB risk of patients was used to categorize them into four groups. The plot obtained indicated that TMB did not impact the predictive precision of the risk score (Figure 9D). The results confirm that risk scores are able to accurately forecast the influence of immunotherapy and function as a dependable indicator of prognosis.

Risk signature in TIME context

In order to examine the correlation between the risk score and tumor immune microenvironment (TIME), we utilized seven distinct approaches, offering a thorough examination of the infiltration of immune cells and stromal cells in connection with the risk score (Figure 10A). The XCELL examination indicated a favorable association between the microenvironment rating and the levels of CRYAB and PDGFRA expression (Figure 10B, 10C). This suggests that the microenvironment score is directly related to the expression of these genes.

Figure 10. Tumor-infiltrating cells abundance estimation. (A) The significantly stronger correlation exists in the high-risk group as shown by the Spearman correlation analysis. (B, C) The correlation between the microenvironment and CRYAB and PDGFRA expression. (D) The TME score in the low and nigh risk groups.

Moreover, the ESTIMATE technique demonstrated a significant rise in scores among patients at high risk. Significant variations in ESTIMATE scores, immune scores, and stromal scores were observed between the high-risk and low-risk groups, as depicted in Figure 10D. The observed pattern indicates a robust correlation between elevated risk scores and heightened participation of immune and stromal cells in the tumor microenvironment, emphasizing the potential importance of these discoveries in comprehending and controlling bladder cancer.

Biological function and signal pathway enrichment analysis

The GSVA yielded valuable information about the functional traits of the high-risk cohort. Our findings indicated that individuals in this group exhibited enhanced activity in several key biological pathways. Significantly, there was a rise in apical junction and epithelial-mesenchymal transition pathways. Furthermore, there was evidence of activation in both the MAPK pathway and the chemokine signaling pathway, as depicted in Figure 11A, 11B. The findings indicate a heightened and active tumor microenvironment in the high-risk category, marked by notable pathway activations that may impact the advancement of the disease and the response to treatment.

Figure 11. Enrichment pathways for GSVA. The correlation between pathways (A), immune hallmarks (B) and risk scores in the KEGG database. (C) The correlation between immune checkpoint blockade genes and risk scores to predict immunotherapy response. (D) The chemotherapy efficiency estimation in the IPS score distribution map. (E) Cisplatin, (F) Docelaxel and (G) methotrexate sensitivity in the low and high groups.

Sensitivity of immunotherapy and chemotherapy prediction

Additionally, our efforts encompass forecasting the responsiveness of immunotherapy. A total of 47 genes associated with checkpoint blockade were obtained from Supplementary Table 4. Figure 11C indicated that the risk scores showed a positive correlation with several genes associated with checkpoint blockade, including TNFSF4, PDCD1LG2, LAIR1, CTLA4, and CD86. According to our predictive model, patients belonging to the CTLA4 positive and PD-1 negative groups (Figure 11D) exhibit a higher IPS score when they are classified as low-risk. This suggests that individuals with a minimal risk might be suitable candidates for receiving CTLA4 inhibitors.

Furthermore, we assessed the effectiveness of BLCA when used alongside conventional chemotherapy medications like cisplatin, docetaxel, and methotrexate. Figure 11E, 11F demonstrated a notable decrease in IC50s among high-risk patients for both Cisplatin and Docelaxel. This suggests that individuals with a greater susceptibility are more receptive to the chemotherapy medications paclitaxel and cisplatin compared to those with a lower susceptibility. Figure 11G shows that patients classified as low-risk had a more favorable response to methotrexate therapy compared to those classified as high-risk. This suggests that patients with a lower risk profile exhibit a greater response to methotrexate compared to patients with a higher risk profile.