An Oxidative Stress-Related Genes Signature for Predicting Survival in Bladder Cancer: Based on TCGA Database and Bioinformatics

Introduction

Bladder cancer is the tenth most prevalent kind of cancer globally, with a morbidity rate of 3.0% and a death rate of 2.1%.¹ Approximately 75% of bladder cancer cases are classified as non-muscle invasive bladder cancer (NMIBC), and 25% are classified as muscle invasive bladder cancer (MIBC).² Although NMIBC can be treated by surgery, it often relapses further after resection, and approximately 15% of cases deteriorate to MIBC, which is related to more aggressive characteristics and poor prognosis.³ The traditional pathological grade and clinical TNM stage are currently utilized to evaluate the prognosis of BC patients. BC with early metastasis is associated with poor prognosis and poor therapeutic outcomes owing to the lack of reliable screening techniques for early diagnosis. Some studies have reported that the mutation of tumor suppressor genes to oncogenes, such as HER-2, H-Ras and Bcl-2, promotes the occurrence of BC.⁴ In recent years, the advent of novel molecularly targeted medicines opens up new avenues for the treatment of advanced BC. Targeted molecular therapy has progressively emerged as a viable option in cancer treatment. Many research has concentrated on developing more effective prediction signatures that might better elucidate the variables that contribute to the prognosis and progression and even more specific immune landscape as well efficacy of immunotherapy of BC, with the goal of providing more evidence for individual molecular therapy methods. Although BC has been extensively studied, few validated biomarkers or screening tools general exist in clinical. Hence, there exists an urgency to develop a robust biomarker to predict prognosis and guide therapeutic treatment, particularly in the areas of targeted therapy and immunotherapy for BC patients.

Oxidative stress is a pathological phenomenon that occurs when the synthesis of oxidants and antioxidants is imbalanced, leading in the buildup of high levels of reactive oxygen species (ROS), which are a potentially critical component driving carcinogenesis and cancer development.^5–7 ROS include oxygen ions, free radicals and peroxides.⁸ Previous research has been that when the scavenging capacity decreases, excessive ROS can injure the DNA, causing genotoxicity and eventually inducing genomic alterations that begin cancer.^9–12 Modest ROS level might play a crucial role in tumor development by disrupting DNA sequences.¹³ Many key genes and gene products have been involved in this system and processes of oxidative stress. However, there has been no comprehensive investigation aiming at identifying particular oxidative stress-related genes (OSRGs) that connect with tumorigenesis, advance and prognosis of BC. During the last several years, database-based bioinformatic expression profile analysis has been frequently used to identify target biomarkers of diagnostic and prognostic significance in different tumors.^14,15 Studies have shown that OSRGs can be used as a predictor of gastric cancer¹⁶ and pancreatic cancer.¹⁷ In urologic neoplasms, oxidative stress pathway genes were reported to be associated with the prognosis of prostate cancer.¹⁸ Nonetheless, the role of OSRGs in BC prognosis prediction remains mainly unclear, and the underlying processes need to be validated further.

This paper concentrates on exploring the possible roles of OSRGs in cancer development and progression as well tumor immunology in BC. Within those differentially expressed OS-related genes (DEOSGs), we identified eight hub OSRGs with non-negligible correlation with BC prognosis. Furthermore, we constructed a prognostic OSRGs signature and established a novel nomogram to better predict the overall survival of BC individuals. We also performed immune cell infiltration and tumor microenvironment analyses in different high- and low-risk group. In addition, ssGSEA analysis was used for immune function evaluation, and drug sensitivity analyses was also carried out. As far as we know, this is the first OS-associated risk signature for prognostic prediction in BC patient, and it may provide fresh insight into the pathophysiology of BC, allowing for more customized treatment and improving BC patients’ prognosis.

Materials and Methods

Data Acquisition and Differentially Expressed OS Genes (DEOSGs)

After downloading the data including RNA sequencing of 414 BC and 19 normal bladder tissues and the clinical information of BC patients from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) on April 14, 2021. We extracted 1399 OS protein domains all of whose relevance scores were higher than 7 (Supplement File Table 1) from GeneCards (https://www.genecards.org) to identify DEOSGs. According to the numerous previous studies, this study defined those DEOSGs with log2 fold change (FC)|≥1 and false discovery rate (FDR)<0.05 as differently expressed one.¹⁹ Meanwhile, genes with an average count value of < 1 were eliminated. This yielded 379 DEOSGs for further analyses. Besides, we established a validation cohort (GSE15307) contained 165 patients with RNA sequencing and corresponding clinical data from GEO database (https://www.ncbi.nlm.nih.gov/geo/). All detailed patient characteristics are summarized in Table 1.

Table 1 Characteristics of Patients from TCGA and GSE13507 Cohort

GO and KEGG Pathway Enrichment Analyses

GO enrichment and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis^20–22 were perform to explore the effect of those identified DEOSGs with the help of several packages including “clusterProfiler”, “ ggplot2 ”, “enrichplot” and “ org.Hs.eg.db ”in R. P and FDR values < 0.05 represent statistically significant.

Protein–Protein Interaction (PPI) Network Construction and Critical Module Screening

To further explore the Interaction of all DEOSGs we constructed PPI network based on the STRING database (http://www.string-db.org/; medium confidence=0.400; hide disconnected nodes in the network),²³ and subsequently, their interactions were imported into the Cytoscape 3.7.2 software to construct and visualize a PPI network. The Molecular Complex Detection (MCODE) plug-in was applied to identified those critical modules in the PPI network with both MCODE scores and node counts >5.²⁴ P value were less than 0.05 indicating a significant difference.

Prognostic Signature Construction and Efficacy Evaluation

Via the survival R package, we integrated those selected DEOSGs in the PPI network into univariate Cox regression analysis to investigate their impact to patient’s survival outcome and survival time, and the candidate genes were identified as prognosis related DEOSGs (P<0.005). Subsequently, these candidate OS genes were added into the least absolute shrinkage and selection operator (LASSO) regression²⁵ and multivariate Cox regression analyses together to identified the hub OS genes that were used to construct a risk score in BC. Patients were divided into two subgroups (low-risk group and high-risk one) according to the median risk score. Next, we performed Kaplan-Meier test and Log-Rank test to analyze the overall survival of two group. Additionally, using the Survival ROC and timeROC package, we carried out a receiver operating characteristic (ROC) curve analysis to evaluate the overall predictive performance of the model.²⁶ We then applied univariate and multivariate Cox regression analysis to investigate the correlation and independence of clinical characteristics and risk score. At last, we developed a nomogram based on OS-associated risk signature as a prognostic tool for patients with BC.²⁷ All above methods were also performance to the GSE15307 cohort to assess the prognostic discrimination and accuracy of this model.

Gene Set Enrichment Analyses

Via GSEA software (version 4.0.3; https://www.gsea-msigdb.org/gsea/), gene set enrichment analyses were performed to the high- and low-risk group to explore the enriched functional pathway of those selected genes to infer their impact. Gene sets with nominal P value < 0.05 and the FDR q value < 0.05 were taken as significant ones.

Hub Gene Evaluation

To clarify the differential expression of eight hub DEOSGs, all of them were verified in TCGA cohort. In addition, we observed differences in the expression levels of eight prognosis-associated hub OS genes between patients with different stages and grade included in the TCGA and GEO cohort. The survival analyses about these eight hub DEOSGs was assess by using GEPIA2 database (http://gepia2.cancer-pku.cn/#index).

Evaluation of Immune Cell Infiltration and Generation of Tumor Microenvironment (TME) Score

With the aid of CIBERSORT, a full-fledged immune cell estimating analysis method, we assessed the proportions of 22 immune cell subtypes based on TCGA cohort. To explore the degree of immune cell infiltration in different risk groups, ESTIMATE algorithm, of which R script was attained from the website (https://sourceforge.net/projects/estimateproject/) helped us to generate the estimate scores, immune scores, and stromal scores for further analyzing the TME.²⁸ The results with P < 0.05 were used for further analyses. Spearman correlation was utilized to test the correlation between risk score and TME scores. The association between risk score and immune infiltration subtype was tested by 2-way ANOVA analyses.

Immune Checkpoint, m6A Genes and Single-Sample Gene Set Enrichment Analyses

ssGSEA helped quantify several tumor-infiltrating immune cell subcohorts between the two subgroups and assess their immune function. Potential immune checkpoint and m6A genes were also retrieved from previous literature, and were performed difference analyzing in the two risk groups. Applying “GSVA” R package, we performed ssGSEA analysis to acquire the infiltration scores of 16 immune cells and evaluate the activities of 13 immune-related pathways of both the high- and low-risk groups.

Tumor Stem Cell Analyses and Immune Subtype Analyses

We investigated the stem cell-like features of tumor by analyzing the tumor stem cell features of extracted TCGA tumor samples in both transcriptome and epigenetics.²⁹ We performed Spearman correlation test to research the interrelationship of tumor stemness and risk score. In addition, immune subtype analyses including C1, C2, C3 and C4 subtype were also performed.

Chemotherapy Sensitivity Analyses

We explored the NCI-60 database (https://discover.nci.nih.gov/cellminer) which contained 60 different cancer cell lines of 9 different types of tumors using the Cell Miner interface. We carried out correlation analyses to potency of 216 selected drugs which attained confirmation of FDA or were in the stages of clinical testing. We also performed Pearson correlation analysis to research the interrelationship of expression levels of these hub OSRGs and drug susceptibility.

Results

Differential Expression Analyses and Functional Enrichment Analyses

The workflow of this study was showed in Figure 1. A total of 379 OS genes, comprising 166 downregulated and 213 upregulated genes, met the screening criteria and were identified as DEOSGs. The expression proportion of these DEOSGs was displayed in the volcano plot and heatmap (Figure 2A and B). We sought to further investigate the underlying molecular mechanisms by which these identified DEOSGs their functional effects on progression of BC using GO and KEGG analyses (Figures 3A and B; 4A and B). Enrichment analyses results revealed that these DEOSGs were mainly enriched in tumor-related pathways, such as the MAPK signaling pathway, TNF-signaling pathway, P13K-Akt signaling pathway and oxidative stress pathway.

Figure 1 The flowchart of this study.

Figure 2 Identification of differently expressed oxidative stress genes (DEOSGs). (A) Volcano plot of DEOSGs between bladder cancer and normal bladder samples. Dots in green represent down-regulated genes, dots in red represent up-regulated genes, and dots in black represent unchanged genes. (B) Heatmap of DEOSGs. Abbreviations: N, normal; T, tumor.

Figure 3 Functional enrichment analysis of up-regulated DEOSGs. (A) Top 10 classes of GO enrichment terms in biological process (BP), cellular component (CC), and molecular function (MF). Circle diagram which enriched in GO analysis. (B) Top 30 classes of KEGG enrichment terms. Circle diagram which enriched in KEGG analysis.

Figure 4 Functional enrichment analysis of down-regulated DEOSGs. (A) Top 10 classes of GO enrichment terms in biological process (BP), cellular component (CC), and molecular function (MF). Circle diagram which enriched in GO analysis. (B) Top 30 classes of KEGG enrichment terms. Circle diagram which enriched in KEGG analysis.

Prognosis-Related Candidate Genes Screening and Construction of OSRGs Signature for BC

Using Cytoscape software and data from the STRING database, the PPI network and top 1 critical modules from PPI network were constructed (Figure 5A and B). A total of 26 genes were identified both in PPI network and univariate Cox regression analyses as candidate OS genes with P < 0.005 (Figure 6A). Subsequently, eight prognosis hub genes were selected by the LASSO analyses and multivariate Cox regression analyses (Figure 6B-E, Table 2). Significant variance between the two groups was shown in the PCA plot (Figure 6F). The risk score of this OSRGs signature was calculated using the following linear formula: Riskscore=∑Expgenei*βi, where Expgene represents the relative expression value of OS genes, and β represents the regression coefficient. Patients were categorized into low- and high-risk groups using the medium risk score.

Table 2 Eight Prognosis-Associated Oxidative Stress Genes with BC in the TCGA Dataset Were Identified

Figure 5 PPI network and modules screening. (A) PPI network of DEOSGs; (B) critical module from PPI network. Green circles represent down-regulated genes, and red circles represent up-regulated genes.

Figure 6 Construction of prognostic model in the TCGA and GSE13507 cohort. (A) Univariate Cox regression analysis for identification prognosis-associated OS genes. (B–D) LASSO analysis and multivariate Cox regression analysis of prognosis-associated OS genes together to constructing the prognosis prediction model. (E) Histogram of prognostic key genes. (F) Principal component analysis (PCA) of the total RNA expression profile of high group (red) and low group (blue).

Hub Genes Evaluation Analyses

The boxplot analyses showed that all these eight OS genes were statistically different in BC and normal samples (Figure 7). In addition, the high expression of these genes was related to the poor prognosis of BC patients, except for AGER (Figure 8A-H). Expressive analyses revealed that LDLR and MSRB3 had higher expression in stage III–IV and high grade; and at the same time, RUNX2 was higher expression in stage III–IV, AGER was higher expression in Low grade, which mean AGER might be a good prognostic factor (Supplement File Figure S1A-C). Taken together, these hub OS genes might play a very important role in the diagnosis and prognosis of BC.

Figure 7 The mRNA expression pattern of oxidative stress-associated OS genes in TCGA cohort. (*p < 0.05; **p < 0.01; ***p < 0.001).

Figure 8 GEPIA Survival analysis of NCAM1 (A), FGF1 (B), AGER (C), RUNX2 (D), MSRB3 (E), LDLR (F), HSPG2 (G), HSPB2 (H).

Evaluation and External Verification of Prognostic Risk Signature Performance

The risk chart shows that the survival rate and survival time of patients in the high-risk group are lower, while the survival rate and survival time of patients in the low-risk group are opposite (Figure 9A-J). Similar results were verified in GSE13507 cohort. These results further proved that our signature effectively distinguished between high- and low-risk patients. Survival curve results indicate that patients in the high-risk group had worse survival outcomes (Figure 10A and C). When time-dependent ROC curves were utilized to assess the risk signature’s discrimination, the risk signature demonstrated good accuracy in predicting 7-year overall survival (Figure 10B). Similar result was observed in the GSE13507 cohort by reaching an AUC of 0.695 at 7 years (Figure 10D), which indicated that this prognostic signature could achieve a more balanced sensitivity and specificity. In both the TCGA and GSE13507 cohorts, the risk score was found to be an independent predictive risk factor in multivariate Cox regression analysis (Figure 10E-H). The signature analyzed using the ROC curve over 5 years showed good predictability in the TCGA and GSE13507 cohort (Figure 10I and J). Moreover, the correlation between risk score and each clinicopathological characteristic was evaluated. We found that the risk score of high-grade tumor or tumor stage III–IV was significantly higher than that of low-grade tumor or tumor stage I–II. Patients with lower tumor grade had significantly superior survival rate than those with higher tumor grade in the TCGA and GSE13507 cohort (Figure 10K-M). The heatmap also revealed that the grade factors were significant differences in TCGA and GSE13507 cohort (Figure 10N and O). Taken together, the OSRGs signature exhibited a good predictive ability in different aspects and might be used as a new predictive tool for BC patients.

Figure 9 Prognostic significance of the oxidative stress-related genes (OSRGs) signature in bladder cancer. (A and B) Patient survival in TCGA and GSE13507 cohort. (C and D) Correlations between the prognostic key genes included in the risk model. Positive and negative correlations are indicated in red and green, respectively. (E and F) Survival time in the high- and low-risk patient groups in TCGA and GSE13507 cohort. (G and H) Patient risk scores in TCGA and GSE13507 cohort. (I and J) Heatmaps of gene expression levels in the risk model for the high- and low-risk groups in TCGA and GSE13507 cohort.

Figure 10 Efficacy evaluation of constructed prognostic model. (A) Survival curve of TCGA cohort. (B) TimeROC curves for forecasting overall survival in TCGA cohort. (C) Survival curve of validation cohort. (D) TimeROC curves for forecasting overall survival in validation cohort. Univariate (E) and multivariate (F) Cox regression analysis of the clinicopathological features in TCGA cohort. Univariate (G) and multivariate (H) Cox regression analysis of the clinicopathological features in GSE13507 cohort. ClinicalROC curves of five-years for forecasting overall survival in TCGA (I) and GSE13507 (J) cohort. (K and L) The relationship between the risk scores and different stage as well as different grade in TCGA cohort. The relationship between the risk scores and different grade (M) in GSE13507 cohort. The heatmap shows the distribution of clinicopathological features and OS genes expression in TCGA (N) and GSE13507 (O) set. Columns in blue represent down-regulated genes, columns in red represent up-regulated genes, and columns in white represent unchanged genes. (*p < 0.05; **p < 0.01; ***p < 0.001).

Performance and Verification of Nomogram

A nomogram was constructed to predict 1-, 3- and 5-year BC overall survival by combining the results of the OSRGs signature with clinical attributes (Figure 11A and B). We assessed the prognostic value of this nomogram using time-dependent ROC analyses. The results demonstrated that the area under the ROC curve (AUC) for 5 years OS was 0.714 and 0.778, respectively (Figure 11C and D), which showed the nomogram plot was accuracy and convenient. From the calibration plots, good calibration was observed in both TCGA and GSE13507 cohorts (Figure 11E-H). The calibration curves and ROC curves demonstrated that the nomogram had great potential to be applied in clinical practice. Our analyses showed that the nomogram was well developed and showed a well accurate value for BC outcomes prediction through discrimination capability evaluation.

Figure 11 Construction of nomogram based on the expression of 8 hub OS genes. Nomogram of risk score and other clinical factors for predicting BC 1-, 3-, and 5-year overall survival in TCGA (A) and GSE13507 (B) cohort. ROC curves showed the predictive efficiency of the risk signature in TCGA (C) and GSE13507 (D) cohort. (E and F) The calibration plot of the nomogram in TCGA cohort. (G and H) The calibration plot of the nomogram in GSE13057 cohort.

Immune Cell Infiltration and Tumor Microenvironment Analyses

The results of OS gene pathway enrichment analyses in the TCGA database show that inflammation, immunity and other factors related pathways are enriched. Based on these findings, we evaluated the rate of immune cell infiltration in different risk groups. The immune infiltration histogram showed that plasma cells, T cells regulatory (Tregs), macrophages M2 and neutrophils were significant differences in the high- and low-risk group (Figure 12A). Figures 12B-E showed the expression abundance of these four cell types in different groups and the correlation between these cell types and risk scores. Patients with high-risk scores had higher proportions of macrophages M2 and neutrophils, and patients with low-risk scores had lower proportions of plasma cells, T cells regulatory (Tregs), confirming the roles of oxidative stress-related mRNA in the regulation of tumor immune infiltration. Furthermore, the correlation analyses results suggest that both immune score and stromal score had a positive correlation with risk score, and their scores in the high-risk group were higher than those in the low-risk group, which might indicate that higher immune activity was in the high-risk group (Figure 13A-D). In summary, these outcomes revealed that OSRGs signature may significantly inhibit or enhance the expression of specific immune cell types, which may affect the response to immunotherapy.

Figure 12 Infiltration of oxidative stress-related immune cells. (A) Heat map of oxidative stress risk and immune cell infiltration using in TCGA databases. (B–E) Scatter plots showing the correlations of the expression of eight genes from TCGA and GEO databases with immune cell regulation, showing differences in the expression levels between the different oxidative risk groups (p < 0.05). The blue line in each plot is a fitted linear model indicating the relationship between significantly Immune cells and the risk of oxidative. Pearson coefficients were used to assess the correlation between the two factors. The box plots show the Immune cells whose infiltration is significantly associated with the risk of oxidative (p < 0.05).

Figure 13 Analysis of immune microenvironment. The relationship between risk score of ORSGs signature and immune score (A and C) as well as stromal score (B and D).

Immune Status, Tumor Stemness, Immune Checkpoint and m6A Genes Analyses

To investigate the correlation of risk scores and immune state. We implemented ssGSEA for the evaluation of enrichment scores of several immune cell subgroups, associated features, and pathways. We found that the content of adc, iDCs, pDCs, APC co-inhibition, APC co-stimulation, HLA and MHC class I antigen presentation process in the TCGA cohort was increased in the high-risk group. In addition, in comparison to the low-risk group, all of these immune cells, related functions and pathways were higher in high-risk group, except for type II IFN response, showing the differences in T cell regulation between high- and low-risk groups (Figure 14A and B). In view of significance immune checkpoint inhibitor-based immunotherapies, investigation of the differential expression of immune checkpoints between the high- and low-risk group were performed. We found a substantial difference in the expression of all these immune checkpoints in the different risk groups (Figure 14C). As expected, the expression levels of PD-L1, PD-1, CTLA4, and LAG3 were considerably greater in the high-risk group than in the low-risk group, as predicted, and the expression levels of these immunological checkpoints were positively associated with the risk score (Supplement File Figure S2A-H). In addition, we found the significant difference of m6A genes expression level of LRPPRC, METTL3, YTHDC1, RBMX, FTO ALKBH5, IGFBP3 and IGFBP2 in the two subgroups above (Figure 14D). Tumor stemness can be measured with RNA stemness score based on mRNA expression (RNAss) and DNA stemness based on DNA methylation pattern (DNAss).³⁰ Spearman correlation analyses between risk score and tumor stem cells showed a risk score was considerable correlation of DNAss and RNAss (Figure 15A and B). We also examined the association between risk scores and immune infiltration to better understand the relationship between risk scores and immunological components. Six types of immune invasion have been identified in human tumors, corresponding to tumor promotion and tumor suppression,³¹ namely C1 (wound healing), C2 (INF-g dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet) and C6 (TGF-b dominant). Due to the number of BC samples were very small, we were not performed correlative analyses at C5 and C6 subtype. We analyzed the immune infiltration of BC in the TCGA cohort and correlated it with the risk score, showing that the high-risk score was significantly linked with C2 and the low-risk score with C3 (Figure 15C).

Figure 14 Immune status between different risk groups and the association between risk score and tumor microenvironment (A and B). (C) Expression of immune checkpoints among high and low BC risk groups. (D) The expression of m6A-related genes between high and low BC risk group. (*p < 0.05; **p < 0.01; ***p < 0.001).

Figure 15 The relationship between risk score and DNAss (A), RNAss (B). Comparison of the risk score in different immune infiltration subtypes (C).

Biological Function and Pathway Analyses

KEGG pathway enrichment analyses shows that the bladder cancer pathway is significantly enriched in the high-risk group. Besides, five GO functions were found to enrich in the high-risk group, the error detection rate < 0.05. GSEA using TCGA data of the Hallmarks gene sets indicated that hypoxia, inflammatory response and TNF-signaling pathway were statistically significant programs (Figure 16A-C).

Figure 16 Gene set enrichment analysis of biological functions and pathways. (A) Hallmark gene set. (B) GO, Gene Ontology. (C) KEGG, Kyoto Encyclopedia of Genes and Genomes.

Prognostic Gene Expression and Cancer Cell Sensitivity to Chemotherapy

We investigated the expression of the OS gene in the prognostic center of the NCI-60 cell line, as well as the connection between its level of expression and drug sensitivity. The finding showed that all prognostic genes are linked to the sensitivity of certain chemotherapeutic treatments (Figure 17 and Supplement File Table 2). For example, the increase in HSPG2, FGF1, LDLR and HSPB2 expression is correlated with tumor cells’ response to ifosfamide, oxaliplatin, dexazolidine, etoposide, Belinostat, oxaliplatin, actinomycin D, Mitriamycin, the resistance of methamphetamine, etoposide, etc. is related to the increase. Conversely, enhanced expression of AGER and FGI seems to contribute to the increased drug sensitivity of cancer cells to chemotherapeutic drugs such as Nelarabine, Dexamethasone Decadron, Fludarabine and Zoledronate. The above results suggested that the expression of hub OS genes might affect the sensitivity of tumors to chemotherapy drugs, which provided a new perspective on the treatment of chemotherapy drugs.

Figure 17 Scatter plot of relationship between prognostic hub OS gene expression and drug sensitivity.

Discussion

With the advancement of high-throughput sequencing and bioinformatics, increasingly bioinformatics methods have been used to identify innovative biomarkers involved in tumor progression,³² survival³³ or prognosis³⁴ of malignancies, which might greatly aid the early diagnosis and evaluation of prognosis in malignant tumors. According to recent studies, OS plays an important function in several phases of carcinogenesis and cancer development.^35,36 However, little is known regarding OS’s possible prognostic relevance in BC, and the OS gene signature has not been described as a prognostic marker for BC.

This research analyzed the expression of 1399 OS genes in BC and investigated their impact on the overall survival of patients. In total, 379 DEOSGs were discovered from the TCGA cohort, and related biological pathways and PPI networks for these genes were comprehensively identified. The KEGG pathway enrichment analyses revealed that our identified DEOSGs were not only highly related with bladder cancer prognosis, but also played a critical role in a variety of other tumors, including pancreatic cancer, breast cancer, prostate cancer, melanoma, and colorectal cancer, prompting us to investigate the potential role of OS genes in other tumors. Interestingly, KEGG analyses also indicated that the down-regulated DEGs were significantly related to COVID-19. Furthermore, these DEOSGs were found to be considerably enriched in cytokine activity, aging, hypoxia response, and oxidative stress response, among other biological processes. All of the biological processes identified in the enrichment analyses have been shown to be associated with carcinogenesis and progression. As a result, these genes were found to be significantly linked to BC progression and are crucial for a complete understanding of the mechanism of these DEOSGs. In addition, the prognosis-related DEOSGs were screened using univariate Cox, LASSO, and multivariate Cox regression analyses, and a total of seven hub OS genes: RUNX2, FGF1, LDLR, HSPB2, HSPG2, MSRB3, and NCAM1 have all been linked to a poor prognosis in BC patients with the exception of AGER. Using TCGA expression data, the expression patterns of these 8 hub OS genes on the mRNA level revealed that AGER, RUNX2, and FGF1 were overexpressed in BC tissues, while LDLR, NCAM1, HSPB2, HSPG2, and HSPB2 were downregulated. These results could be related to the previously described modulating effects of these genes on BC metastasis and growth. We found that differential expression of eight hub OSRGs is connected with patients’ overall survival and may play a role in the development of BC.

The PCA analyses revealed that these hub genes have good discriminative power. And then we discovered that being in the high-risk group was significantly associated with a poor outcome, being older, having a high grade, being in the advanced T stage, being in stage III–IV, and having a shorter overall survival period, which was consistent with clinical characteristics’ predictive implications on overall survival. Except for AGER, Kaplan-Meier analyses demonstrated that these hub OS genes were inversely related with overall survival in BC patients. This seemed to mean that the higher the overall expression level of these hub OS genes, the worse the patient’s prognosis. In terms of the relationship between risk score and clinical characteristics, a high-risk score was strongly correlated with tumor high grade or tumor stage III–IV, indicating that a high-risk score is unquestionably connected with a poor prognosis. The risk score obtained from this signature was found to be an independent prognostic factor for BC patients in both univariate and multivariate Cox models. Furthermore, time-dependent ROC curve tests demonstrated the risk score prognostic signature’s strong predictive potential as an independent prognostic factor. And the predictive nomogram based these eight OSRGs shows robust performance in the individual predicting PC prognosis compared with the clinical characteristics in BC. We also validated the applicability of the oxidative stress-related risk signature in the GSE13507 cohort. The bioinformatic analyses results showed that the risk signature constructed by eight OSRGs could be used for overall survival prognosis, with good prognostic value and provides new insights into the molecular mechanism of bladder cancer for follow-up studies.

Notably, the intricate interplay between tumor cells and the tumor microenvironment not only plays a critical role in tumor formation, but it also has a substantial impact on immunotherapeutic effectiveness and overall survival.³⁷ Besides, the tumor microenvironment is crucial in the progression of cancer and its response to immunotherapy. Tumor microenvironment is becoming an increasingly popular topic and could affect tumor progression and recurrence.³⁸ Tumor heterogeneity is dependent on cancer cells and the tumor microenvironment, which is made up of many cell types such as immune cells, mesenchymal stem cells, endothelial cells, and fibroblasts.^39,40 In cancer progression, cancer cells generally induce an immunosuppressive microenvironment to fight against anticancer immunity and it is also the premise and guarantee of tumor immune escape. Cells and molecules in the tumor microenvironment are dynamic, reflecting the essence of tumor microenvironment evolution, with the end result being that a substantial percentage of immunosuppressive cells, such as myeloid-derived suppressor cells, aggregate in the tumor microenvironment to regulate immune escape, tumor growth, and metastasis, and to influence tumor outcome.^28,41,42 A functional enrichment analysis revealed that OSRGs are predominantly engaged in immunological pathways in this research. The quantity of four immune cells, comprising macrophages M2, neutrophils, Tregs, and plasma cells, differed considerably between high- and low-risk groups, validating that OSRGs were contributory in regulating tumor immune infiltration. It is worth noting that Macrophages, a prominent component of the tumor microenvironment, contribute to carcinogenesis by increasing angiogenesis, immunosuppression, invasion, and metastasis.^43,44 Previous research has linked an increase in tumor-associated macrophages, neutrophils, and Treg cells to a poor prognosis of patients with BC on account of their function in immunological invasion. In conclusion, the four survival-related cells discovered above have the potential to be key roles in immune infiltration as well as BC immunotherapy, validating the validity of the cell-based analyses of oxidative stress-related genes. Based on the GSEA analyses, tumor-related signal pathways such as Hypoxia, Inflammatory response and TNFα- signaling pathway were significantly enriched, providing new idea of therapy.

In order further explore the relationship between risk score and immune components, a relevant analysis was evaluated. Interestingly, we discovered a substantial correlation between high-risk score and C2 while a low-risk score was clearly related to C3, implying that C2 promotes tumor incidence and progression and C3 is a good preventive factor.

Immune checkpoint blockade (ICB), which includes anti-PD1/PD-L1 and anti-CTLA-4/LAG3 treatments, has proven to be extremely effective in the treatment of human malignancies, particularly solid tumors. Combining immune checkpoint inhibition with immunotherapies such as CTLA-4, LAG3, PD-1, and PD-L1 inhibitors is a potential strategy for treating a range of cancers, and an activated tumor immune microenvironment is associated with excellent immune checkpoint inhibitor treatment outcomes.^45,46 Previous research has indicated that anti-PD1/PD-L1 antibody treatments are effective in malignancies with a high mutational load, such as Hodgkin’s lymphoma, melanoma, renal cell carcinoma, non-small lung cancer carcinoma, and urothelial bladder carcinoma.^47–52 In our study, we observed enhanced expression of PD-L1 in in the high-risk group, suggesting that patients with high-risk scores may benefit from anti-PD-L1 immunotherapy, the same to other genes signatures and corresponding immune checkpoints. Furthermore, the expression levels of immune checkpoints in the high-risk cohort were substantially greater than those in the low-risk cohort, being agree with our discovery that the positive association between the risk score and the expression level of these four immune checkpoints. And it helps offer novel understanding of tumor immunotherapy. Therefore, the use of biomarkers from OS genes to guide treatment decisions along with development of immunotherapy may benefit the patient. Moreover, increased activities of 16 kinds of immune cells and 12 kinds of immune functions in the high-risk group indicated that patients in the high-risk group were in an immune-active state. As a result, it is fair to conclude that the high-risk group’s anti-tumor immunity is compromised, which may be a contributing factor to its poor prognosis. Researches focused on the N6-methyladenosine (m6A) RNA methylation level and its impact on tumor were increasingly reported over the last few years. In this research, analyzed the different expression level of m6A-related genes in high-risk groups and the low ones and the results indicated that LRPPRC, METTL3, YTHDC1, RBMX, FTO, ALKBH5, IGFBP3 and IGFBP2 were significance.

Cancer stem-like cells (CSCs) can develop from a variety of sources, including long-lived stem or progenitor cells or by dedifferentiation from non-stem cancer cells that convert to CSCs via dysregulation of associated signaling pathways.³⁰ CSCs may possess capability for self-renewal and multipotent differentiation, which leads to drug resistance.^53–55 In our study, the risk score was positively correlated with cancer stem-cell-like features (RNAss and DNAss), with the particular mechanism warranting further investigation. According to ESTIMATE algorithm, both immune score and stromal score were positively linked to the risk score, which suggested that the tumor tissue in the high-risk group is heavily invaded by immune cells.

We discovered that higher expression of certain prognostic genes was related with enhanced treatment resistance for a variety of FDA-approved chemotherapeutic medicines by analyzing data from NCI-60 cell lines. Up-regulated expression level of HSPG2, FGF1, LDLR and HSPB2 had been linked to increased drug resistance of cancer cells to Ifosfamide, Oxaliplatin, Dexrazoxane, Etoposide, Belinostat, Oxaliplatin, Actinomycin D, Mitriamycin, Eribulin mesylate, Etoposide, etc. The hub OS genes found for each drug were significantly diverse, suggesting that drug resistance genes are extremely drug-specific. These findings show that some prognostic genes can be utilized as therapeutic targets to overcome treatment resistance or adjuvant medication sensitivity.

The present study aimed to identify whether OSRGs play predictive roles in BC patients. Therefore, all pathologically confirmed bladder cancer samples were included in the analysis without distinguishing between MIBC and NMIBC. There were 165 bladder cancer samples with available clinical data in the GEO cohort, including 103 (62.4%) NMIBC samples and 62 (37.6%) MIBC samples. For the TCGA cohort, there were only four NMIBC samples (1%) in this study. If only MIBC or NMIBC samples were included in the analysis, the patient numbers would have been too small for reliable statistics. Therefore, larger-scale gene sequencing validation studies on NMIBC alone or MIBC alone are needed in the future. It is worth noting that the risk group in the present study was significantly correlated with MIBC or NMIBC status in both the GSE13507 (p = 0.001) and TCGA (p = 0.005) cohorts (Table 3).

Table 3 Patient Characteristics in Different Risk Groups

Nonetheless, this study has some limitations that should be noted. First, because this study was intended as a retrospective analysis, further prospective investigations should be conducted to validate these findings. Second, there is no in vitro or in vivo testing to corroborate the validity of the mechanistic analyses. Third, since we do not have our own clinical data of BC, although the results in this study have certain credible, they still need to be corroborated with our own clinical data. As a result, a variety of tests will be carried out in the future to establish the molecular links between these genes and BC development.

Conclusion

In summary, we constructed an OSRGs signature for BC to predict prognosis. We established an effective nomogram including eight hub OS genes. Importantly, the OSRGs signature generated in our study might be associated with immune infiltration levels, tumor microenvironment, drug sensitivity and even the efficacy of tumor immunotherapy.