Identification of molecular subtypes and prognostic risk model of glucocorticoid-related lncRNAs in bladder cancer to evaluate prognosis and immunological characteristics

Introduction

Bladder cancer (BC) ranks among the most frequently observed malignancies impacting the urinary tract, distinguished by its elevated relapse rates and significant mortality, accounting for nearly 550,000 new diagnoses and 170,000 deaths each year (1,2). Despite progress in surgical techniques and treatment strategies, patients diagnosed with advanced BC face a generally poor prognosis, with high recurrence rates and limited response to immunotherapy in certain patient subsets (3). Consequently, the discovery of reliable biomarkers to improve diagnostic accuracy, prognosis prediction, and targeted intervention is critically important.

Long non-coding RNAs (lncRNAs), which consist of structurally diverse RNA transcripts exceeding 200 nucleotides and lacking protein-coding potential, fail to code for proteins but exert important effects on regulating gene expression at a range of levels, thereby influencing fundamental biological processes (BPs), namely cell proliferation, differentiation, apoptosis, and immune responses (4). Mounting evidence consistently demonstrates that dysregulated lncRNAs contribute intimately to the initiation, progression, and metastasis of a wide array of human cancers, including BC (5,6). For instance, specific lncRNAs have been reported to act as oncogenes or tumor suppressors in BC, impacting disease prognosis and contributing to therapeutic resistance (7,8).

Glucocorticoids play crucial roles in modulating inflammation, cellular stress responses, and cancer-associated signaling pathways comprising PI3K-Akt and cytokine-cytokine receptor interaction pathways often dysregulated in BC (9-11). While immune-related lncRNA signatures have demonstrated predictive power for both prognosis and immune infiltration, the prognostic and immunological impact of lncRNAs linked to glucocorticoid biology has not been fully investigated (12). Elucidating these glucocorticoid-related lncRNAs (GR-lncRNAs) may reveal new prognostic biomarkers and refine our understanding of BC’s tumor-immune landscape.

In this investigation, our objective was to systematically determine GR-lncRNAs in BC and develop a robust prognostic risk model. We further investigated the association between our constructed risk model, identified molecular subtypes, immune cell infiltration patterns, tumor mutational burden (TMB), and drug sensitivities. Our findings provide an in-depth insight into the profile of GR-lncRNAs in BC, introducing an innovative tool for predicting patient outcomes and a foundation for developing more personalized and optimal therapeutic interventions. We present this article in accordance with the TRIPOD reporting checklist (available at https://tau.amegroups.com/article/view/10.21037/tau-2025-528/rc).

Methods

Data acquisition

We accessed BC messenger RNA (mRNA) expression, mutation, and clinical data from The Cancer Genome Atlas (TCGA) repository (https://portal.gdc.cancer.gov/). Samples exhibiting overall survival under 30 days were removed from the analysis. A combined total of 335 tumor samples and 32 normal samples were considered in the following analysis. A random partition of the BC dataset allocated 70% of the samples to the training set and 30% to the validation set. Additionally, glucocorticoid-related genes (GRGs) data were acquired from GeneCards (https://www.genecards.org/) by inputting “Glucocorticoid” and including only those genes with scores greater than 2. This yielded 939 GRGs. Normal samples correspond to TCGA solid tissue normal specimens (sample-type code 11), representing adjacent non-tumor bladder tissue. Information on neoadjuvant chemotherapy (NAC) was not uniformly available in the TCGA-BLCA cohort; therefore, NAC status could not be incorporated into the present analysis. In the TCGA-BLCA dataset, tumor specimens are annotated only as ‘Primary Tumor’ without specifying whether they were obtained from transurethral resection of bladder tumor (TURBT) or radical cystectomy; therefore, the exact surgical source of the samples could not be determined.

Identification of differentially expressed GRGs (DEGRGs) and GR-lncRNAs

The limma package in R was applied to detect differentially expressed genes (DEGs) within normal and tumor BC samples. Differential expression was determined based on |log₂ fold change (FC)| >0.585 and false discovery rate (FDR)-adjusted P<0.05. The intersection of these DEGs with the collected GRGs identified 230 DEGRGs. Subsequently, differentially expressed lncRNAs (DElncRNAs) were identified utilizing the limma package with thresholds of |log₂ FC| >0.585 and adjusted P value <0.05. To identify potential GR-lncRNAs, Spearman correlation coefficients were calculated within the expression profiles of DEGRGs and DElncRNAs. LncRNAs with |cor| >0.4 and P value <0.05 were considered GR-lncRNAs.

Establishment of a prognostic model relying on GR-lncRNAs

The identified GR-lncRNAs were subjected to univariate Cox regression analysis utilizing the survival R package to assess for prognostically significant lncRNAs (P value <0.05). To minimize the risk of overfitting, the least absolute shrinkage and selection operator (LASSO) regression analysis was applied for the candidate lncRNAs employing the glmnet package. The optimal penalty parameter λ was determined via cross-validation, reducing model complexity by removing highly correlated lncRNAs. Subsequently, using the survival package, a multivariate Cox regression was carried out on candidate lncRNAs chosen by LASSO to develop the prognostic model.

A risk score was determined individually per patient employing the expression patterns of the selected lncRNAs and their relevant risk coefficients. Risk-based classification separated patients into high- and low-risk groups employing the median risk score as a cutoff. Kaplan-Meier survival analyses were carried out via the survival package to assess differences in overall survival across the high- and low-risk groups. Receiver operating characteristic (ROC) curves were illustrated using the timeROC package, and the area under the curve (AUC) values were estimated to investigate the model’s predictive accuracy. Risk score distribution plots and survival status distribution plots for high- and low-risk groups were also generated. The prognostic model was additionally validated utilizing the validation set, with corresponding survival curves, ROC curves, risk score distribution, and survival status distribution plots. Kaplan-Meier survival curves for each of the seven characteristic lncRNAs in the training set, validation set, and the entire cohort were also generated.

Construction of a nomogram for independent prognostic analysis

To investigate the potential of the risk score from the prognostic model to independently predict patient outcomes, both univariate and multivariate Cox regression analyses were carried out on the training cohort by integrating clinical variables and the risk score. A nomogram was constructed using the rms R package to predict the overall survival probabilities in BC patients. To determine the accuracy of the nomogram’s predictions, calibration curves were plotted comparing predicted versus observed survival outcomes. The decision curve analysis (DCA) curves were plotted to evaluate the practical value of the nomogram in guiding clinical decision-making.

Immune infiltration analysis

Immune cell infiltration was quantified employing single-sample gene set enrichment analysis (ssGSEA) (29 immune features) and CIBERSORT algorithms. Comparisons were made between the risk groups for ESTIMATE scores (calculated using the estimate R package), immune checkpoint (ICP) expression, and immunophenoscore (IPS). The classification of tumor samples included six immune subtypes: C1 (wound healing), C2 (IFN-γ dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-β dominant). TCGA BC subtypes were downloaded from a previous publication (13). Tumor Immune Dysfunction and Exclusion (TIDE) scores were obtained from the TIDE database (http://tide.dfci.harvard.edu/) and compared between the high- and low-risk groups.

Enrichment analysis

Gene set enrichment analysis (GSEA) was carried out using GSEA software to identify enriched pathways within the high- and low-risk groups based on their prognostic scores. DEGs between high- and low-risk groups were recognized applying the limma package (|log FC| >0.585, adjusted P value <0.05). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were then accomplished on these DEGs via the clusterProfiler and GOplot R packages.

Construction of GR-lncRNA BC subtypes

Based on the lncRNAs incorporated in the prognostic model, BC samples were clustered into two subtypes through the non-negative matrix factorization (NMF) R package. Immune infiltration analysis was performed for these subtypes. DEGs across the two subtypes were recognized via the limma package, followed by GO and KEGG enrichment analysis.

TMB analysis

TMB scores were assessed for individual BC samples applying TCGA mutation information. Wilcoxon rank-sum tests were performed to compare TMB values across high- and low-risk groups. Waterfall plots illustrating the 20 most frequently mutated genes within both risk groups were constructed by applying the GenVisR R package.

Drug sensitivity prediction

To discover possible therapeutic targets and efficient anti-tumor drugs, the CellMiner database (https://discover.nci.nih.gov/cellminer/) served as a resource for screening anti-tumor compounds whose sensitivity was markedly linked to the levels of the prognostic lncRNAs. The half-maximal inhibitory concentration (IC₅₀) of various drugs within high- and low-risk groups was forecasted via the pRRophetic R package. A lower IC₅₀ indicates greater drug efficacy in inhibiting tumor growth.

Statistical analysis

All statistical analyses were performed using R software (version 4.4.1) and its associated packages. The Wilcoxon test was used to determine statistical differences between two groups. Kaplan-Meier survival curves were analyzed using the log-rank test to assess differences in survival between groups. ROC curves were generated using the timeROC package. Data visualization was primarily conducted with the ggplot2 package. Correlations between two variables were evaluated using Spearman’s and Pearson’s correlation analyses, as appropriate. The median value was used as the cutoff point for group stratification. Statistical significance was defined as a P value <0.05. The levels of statistical significance were summarized as follows: ****, P<0.0001; ***, 0.0001<P<0.001; **, 0.001<P<0.01; *, 0.01<P<0.05; ns, P>0.05.

Ethical statement

This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Results

Construction of prognostic model with differentially expressed GR-lncRNAs in BC patients

Initially, 230 DEGRGs were determined through the intersection of GRGs and DEGs across tumor and normal samples (Figure 1A). GO enrichment analysis of DEGRGs demonstrated notable enrichment in BPs, including response to steroid hormone, gland development, and epithelial cell proliferation (Figure 1B). KEGG enrichment analysis indicated that DEGRGs were highly presented in the PI3K-Akt signaling pathway, cytokine-cytokine receptor interaction, and human T-cell leukemia virus 1 infection pathways (Figure 1C).

Figure 1 Enrichment analysis of DEGRGs in BC. (A) Upset plot showing the intersection of GRGs and DEGs. (B) GO enrichment analysis of DEGRGs. (C) KEGG enrichment analysis of DEGRGs. BC, bladder cancer; BP, biological process; CC, cellular component; DEGs, differentially expressed genes; DEGRGs, differentially expressed glucocorticoid receptor-related genes; GO, Gene Ontology; GCG, glucocorticoid-related gene cluster; GRGs, glucocorticoid receptor-related genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

We then identified 248 DElncRNAs between tumor and normal samples (table available at https://cdn.amegroups.cn/static/public/tau-2025-528-1.xlsx). A co-expression analysis between DEGRGs and DElncRNAs resulted in the identification of 151 GR-lncRNAs (table available at https://cdn.amegroups.cn/static/public/tau-2025-528-2.xlsx). Univariate Cox regression analysis was carried out on these 151 GR-lncRNAs, identifying those strongly correlated with overall survival. Subsequently, LASSO regression analysis was applied, which selected 15 characteristic lncRNAs (Figure 2A,2B). Multivariate Cox regression analysis was then utilized for these 15 lncRNAs, finally identifying seven characteristic lncRNAs to construct a robust prognostic model (Figure 2C):

Riskscore=−0.236∗AL390728.6−0.224∗AL355353.1−0.261∗AC011477.2−0.581∗AL357033.4+0.138*AC073210.3+0.321∗AC105942.1+0.775∗FRMD6.AS2[1]

Figure 2 Construction and identification of prognostic GR-lncRNAs. (A) Coefficient profiles from LASSO regression analysis. (B) Cross-validation curve for optimal lambda selection. (C) Multivariate Cox regression analysis identifying the seven final characteristic genes. (D) ROC curves for the prognostic model in the training set, validation set, and entire set. (E) Kaplan-Meier survival curves for the prognostic model in the training set, validation set, and entire set. (F) Scatter plots showing risk score distribution and survival status in the training set, validation set, and entire set. AUC, area under the curve; CI, confidence interval; GR-lncRNAs, glucocorticoid receptor-related long non-coding RNAs; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

ROC analysis demonstrated that the prognostic model exhibited stable and robust predictive performance across the TCGA training cohort, validation cohort, and the entire cohort (Figure 2D). Survival analysis further demonstrated that patients within the low-risk group consistently exhibited markedly improved overall survival relative to the high-risk group across all datasets (Figure 2E). The distribution scatter plot was also generated for the training, validation, and full cohorts (Figure 2F).

Kaplan-Meier survival curves for the individual lncRNAs showed that high expression of AL355353.1, AL357033.4, AL390728.6, and AC011477.2 was associated with improved survival outcomes, whereas low expression of AC073210.3, AC105942.1, and FRMD6-AS2 correlated with better prognosis (see Figure S1). Analysis of expression patterns between tumor and normal tissues indicated that AL390728.6, AL355353.1, AC073210.3, and AC011477.2 were significantly upregulated in tumor samples, while AL357033.4, AC105942.1, and FRMD6-AS2 were predominantly expressed in normal tissues (see Figure S2).

Enrichment analysis

The KEGG pathway enrichment analysis was carried out applying GSEA software to explore the functional differences across the high- and low-risk groups. In the high-risk group, GSEA revealed significant accumulation of pathways connected with arrhythmogenic right ventricular cardiomyopathy (ARVC), dilated cardiomyopathy, ECM-receptor interaction, and focal adhesion (Figure 3A). Conversely, the low-risk group presented positive overrepresentation in several metabolic and biosynthetic pathways, including glycerophospholipid metabolism, linoleic acid metabolism, peroxisome, and ribosome pathways (Figure 3B). Differential expression analysis between the two risk groups (|log FC| >1, adjusted P value <0.05) identified 199 upregulated genes and 8 downregulated genes. KEGG and GO enrichment analyses were then carried out on these DEGs. KEGG pathway enrichment analysis exhibited that the main enrichment of upregulated genes occurred in the cytoskeleton in muscle cells and focal adhesion pathways (Figure 3C), while the majority of downregulated genes were associated with PPAR signaling pathway and arachidonic acid metabolism pathways (Figure 3D). GO analysis of upregulated genes revealed enrichment in external encapsulating structure organization, extracellular matrix organization, and extracellular structure organization in BP; collagen-containing extracellular matrix, endoplasmic reticulum lumen, and collagen trimer in cellular component (CC); and extracellular matrix structural constituent, collagen binding, and glycosaminoglycan binding in molecular function (MF) (Figure 3E). For downregulated genes, GO analysis indicated enrichment in fatty acid metabolic process, icosanoid metabolic process, and fatty acid derivative metabolic process (BP); Golgi cisterna, apical plasma membrane, and basolateral plasma membrane (CC); and tetrapyrrole binding, monooxygenase activity, and aromatase activity (MF) (Figure 3F).

Figure 3 Functional enrichment analysis of DEGs between high- and low-risk groups. (A) GSEA results for the high-risk group. (B) GSEA results for the low-risk group. (C) KEGG enrichment of upregulated genes in high- and low-risk groups. (D) KEGG enrichment of downregulated genes in high- and low-risk groups. (E) GO enrichment results for upregulated genes across BP, CC, and MF processes. (F) GO enrichment results for downregulated genes, across BP, CC, and MF processes. BP, biological process; CC, cellular component; ECM, extracellular matrix; GSEA, gene set enrichment analysis; MF, molecular function.

Construction of a nomogram for independent prognostic analysis

The baseline clinicopathological characteristics of the TCGA-BLCA cohort are summarized in Table 1. Univariate Cox regression analysis was performed by integrating the risk score with clinical parameters. The evidence demonstrated that stage, T stage, M stage, N stage, and risk score demonstrated a notable association with overall survival (Figure 4A). According to the multivariate Cox regression results, the risk score was validated as an independent predictor of prognosis (Figure 4B). A nomogram was constructed by integrating the risk score and clinical features to assess 1-, 3-, and 5-year overall survival probabilities (Figure 4C). DCA demonstrated a favorable net clinical benefit for predicting survival at 1, 3, and 5 years (Figure 4D). Calibration plots demonstrated strong concordance across predicted and observed survival outcomes across 1, 3, and 5 years, indicating excellent model calibration (Figure 4E). Furthermore, comparison of risk scores across clinical stages revealed that a notably higher risk score was observed in stage III–IV patients compared to those in stage I–II (Figure 4F).

Figure 4 Evaluation of the independent prognostic value of the risk model. (A) Forest plot of univariate Cox regression analysis for clinical characteristics and risk score. (B) Forest plot of multivariate Cox regression analysis confirming risk score as an independent prognostic factor. (C) Nomogram predicting 1-, 3-, and 5-year overall survival based on clinical factors and risk score. (D) DCA for 1-, 3-, and 5-year overall survival, evaluating the clinical utility of the model. (E) Calibration plots demonstrating the agreement between predicted and observed survival rates at 1, 3, and 5 years. (F) Boxplot showing significantly higher risk scores in stage III–IV patients compared to stage I–II patients. ****, P<0.00001. f-1, 1-year survival probability; f-3, 3-year survival probability; f-5, 5-year survival probability. DCA, decision curve analysis; M, metastasis; N, node; T, tumor.

Immune infiltration analysis

The ESTIMATE algorithm revealed that patients in the high-risk group had markedly higher ESTIMATE scores, immune scores, and stromal scores, along with decreased tumor purity, contrasted with the low-risk group (Figure 5A). IPS analysis demonstrated that the low-risk group exhibited notably higher IPS scores related to the high-risk group (Figure 5B). Using the CIBERSORT algorithm, we compared immune cell infiltration across high- and low-risk groups. The observations indicated that most immune cell types, such as dendritic cells, T cells, follicular helper, and CD8 T cells, were more highly expressed within the low-risk group (Figure 5C,5D). To further explore immune status, we conducted ssGSEA analysis, which demonstrated that significantly higher levels of immune cell infiltration were observed in the high-risk group, including macrophages M0, macrophages M1, macrophages M0, and neutrophils (Figure 5E). We additionally evaluated the differential expression of ICP-related genes across the high- and low-risk groups. A large proportion of ICP genes exhibited significantly increased expression in the high-risk group (Figure 5F). Additionally, TIDE scores were significantly higher within the high-risk group, suggesting an elevated potential for immune evasion among high-risk patients (Figure 5G).

Figure 5 Comprehensive immune landscape analysis between high- and low-risk groups in BC. (A) ESTIMATE algorithm results showing differences in immune score, stromal score, ESTIMATE score, and tumor purity in high- and low-risk groups. (B) Violin plot of IPS indicating higher immunotherapy sensitivity in the low-risk group. (C) Heatmap of immune cell infiltration levels based on the CIBERSORT algorithm. (D) Boxplot of immune cell infiltration levels based on the CIBERSORT algorithm. (E) ssGSEA analysis displaying immune cell types and immune-related functions. (F) Differential expression of ICP genes between the two groups. (G) TIDE scores demonstrating greater immune evasion potential in the high-risk group. *, P<0.05; **, P<0.001; ***, P<0.0001; ****, P<0.00001. BC, bladder cancer; ICP, immune checkpoint; IPS, immunophenoscore; ssGSEA, single-sample gene set enrichment analysis; TIDE, tumor immune dysfunction and exclusion.

TMB and drug sensitivity analysis

TMB was calculated using somatic mutation data from TCGA-BC. Subsequently, Waterfall diagrams were employed to present the most commonly mutated 20 genes identified per group. In the high-risk group, TP53, TTN, and ARID1A showed the highest mutation frequencies (Figure 6A), while TTN, TP53, and MUC16 represented the most frequently mutated genes within the low-risk group (Figure 6B). To further explore therapeutic relevance, we used the CellMiner database to assess correlations across the expression of model-constructed lncRNAs and drug sensitivity. Among the seven lncRNAs, FRMD6-AS2 was negatively correlated with sensitivity to estramustine, cordycepin, and acetalax, but positively correlated with midostaurin (Figure 6C). Detailed results are provided in Table S1. Drug sensitivity prediction exhibited that significantly lower IC₅₀ values were observed in the low-risk group for docetaxel and cisplatin, suggesting higher sensitivity. Conversely, the high-risk group exhibited markedly increased IC₅₀ values for 5-fluorouracil, doxorubicin, imatinib, and sorafenib (Figure 6D).

Figure 6 TMB and drug sensitivity analysis in high- and low-risk groups. (A) The top 20 mutated genes are shown in waterfall plots for high-risk groups. (B) The top 20 mutated genes are shown in waterfall plots for low-risk groups. (C) Correlation of FRMD6-AS2 expression with drug sensitivity identified using the CellMiner database. (D) Predicted drug IC₅₀ values showing differential responses to chemotherapy across risk groups. *, P<0.05; **, P<0.001; ***, P<0.0001. FRMD6-AS2, FERM domain containing 6 antisense RNA 2; IC₅₀, half-maximal inhibitory concentration; TMB, tumor mutational burden.

Subtype clustering analysis based on seven lncRNAs

Based on the seven characteristic lncRNAs from the prognostic model, BC samples were categorized into two separate molecular subtypes, Cluster 1 and Cluster 2, using the NMF package (Figure 7A,7B). Analysis of survival data indicated that Cluster 1 was associated with significantly greater overall survival rates than Cluster 2 (Figure 7C). Consistent with the survival differences, Cluster 2 exhibited markedly higher risk scores than Cluster 1 (Figure 7D). The boxplots were subsequently generated to visualize the expression profiles of the seven signature lncRNAs across the two molecular subtypes (Figure 7E). The results revealed that a large portion of these lncRNAs revealed significant differences in expression across the two subtypes.

Figure 7 Identification of molecular subtypes based on GR-lncRNAs in BC. (A) CDF curve for determining optimal clustering number using NMF. (B) Consensus heatmap displaying two distinct molecular subtypes (Cluster 1 and Cluster 2) based on lncRNA expression. (C) Kaplan-Meier survival curves showing that Cluster 1 is associated with better overall survival than Cluster 2. (D) Boxplot comparing risk scores between the two clusters, with Cluster 2 showing significantly higher risk scores than Cluster 1. (E) Boxplots illustrating significant differences in the expression of most lncRNAs between Cluster 1 and Cluster 2. *, P<0.05; **, P<0.001; ***, P<0.0001; ****, P<0.00001; ns, not significant. BC, bladder cancer; CDF, cumulative distribution function; Exp, expression; GR-lncRNAs, glucocorticoid receptor-related long non-coding RNAs; lncRNAs, long non-coding RNAs; NMF, nonnegative matrix factorization.

Molecular subtype immune infiltration analysis

We first compared the IPS across the two molecular subtypes. The observations revealed that significantly higher IPS scores were observed in Cluster 1 relative to Cluster 2 (Figure 8A). Evaluation utilizing the ESTIMATE algorithm further demonstrated that a significant increase in immune score was observed in Cluster 1 relative to Cluster 2 (Figure 8B). The ssGSEA results showed that most immune-related functions and immune cell types were more highly enriched in Cluster 1 (Figure 8C,8D). Besides, the CIBERSORT analysis evidenced that plasma cells, macrophages M0, and naive B cells were markedly enriched in Cluster 2, while resting memory CD4⁺ T cells were more abundant in Cluster 1 (Figure 8E). Next, we investigated the expression of ICP genes across the subtypes. The observations revealed that a large portion of ICP-related genes were highly expressed in Cluster 1 (Figure 8F). We then analyzed the pattern of the six characterized immune subtypes (C1: wound healing, C2: IFN-γ dominant, C3: inflammatory, C4: lymphocyte depleted, C5: immunologically quiet, C6: TGF-β dominant) within the high- and low-risk groups, together with the two molecular subtypes, showing a clear correspondence between molecular clustering and clinical risk stratification (Figure 8G).

Figure 8 Comprehensive immunological of molecular subtypes in BC based on GR-lncRNAs. (A) Violin plot of IPS scores for the molecular subtypes. (B) Estimate scores for the molecular subtypes. (C) Box plots of immune infiltration differences in molecular subtypes. (D) Heatmap of ssGSEA immune infiltration differences in molecular subtypes. (E) Box plots of immune cell types with significant differential distribution in molecular subtypes. (F) Box plots of ICP expression differences in molecular subtypes. (G) Sankey diagram illustrating the relationship between TCGA immune subtypes, high/low-risk groups, and the two identified subtypes. *, P<0.05; **, P<0.001; ***, P<0.0001; ****, P<0.00001; ns, not significant. BC, bladder cancer; GR-lncRNAs, glucocorticoid receptor-related long non-coding RNAs; ICP, immune checkpoint; IPS, immunophenoscore; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas.

Furthermore, enrichment analyses for GO and KEGG pathways were carried out based on the DEGs between the two molecular subtypes. Analysis of KEGG pathways revealed that significant enrichment of upregulated genes was observed in ECM-receptor interaction, focal adhesion, and small cell lung cancer pathways (Figure 9A), while the majority of downregulated genes participated in cornified envelope formation and metabolism of xenobiotics by cytochrome P450 (Figure 9B). Functional enrichment analysis using GO exhibited that the main enrichment of upregulated genes occurred in response to virus, defense response to virus, and leukocyte migration pathways (Figure 9C), whereas downregulated genes were associated with epidermis development, skin development, and amide metabolic processes (Figure 9D).

Figure 9 Functional characterization of molecular subtypes in BC based on GR-lncRNAs. (A) KEGG enrichment for upregulated genes. (B) KEGG enrichment for downregulated genes. (C) GO enrichment for upregulated genes. (D) GO enrichment for downregulated genes. BC, bladder cancer; GO, Gene Ontology; GR-lncRNAs, glucocorticoid receptor-related long non-coding RNAs; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Discussion

BC is a highly heterogeneous malignancy, and its diverse clinical outcomes necessitate the identification of reliable biomarkers for accurate prognosis prediction and individualized therapeutic approaches (14). In the current study, we conducted a comprehensive investigation into the role of GR-lncRNAs in BC, leading to the construction of a novel prognostic risk model and the identification of distinct molecular subtypes.

Our initial screening successfully identified 230 DEGRGs. Functional enrichment analyses, including GO and KEGG, revealed that these DEGRGs were primarily associated with response to steroid hormone, gland development, epithelial cell proliferation, and key signaling pathways such as PI3K-Akt. This highlights the intricate involvement of glucocorticoid signaling in BC pathogenesis, potentially influencing tumor growth, differentiation, and progression (15). Subsequent co-expression analysis with DElncRNAs led to the identification of 151 GR-lncRNAs, underscoring the extensive and complex regulatory network between lncRNAs and glucocorticoid pathways in BC.

A pivotal finding of our study is the development of a robust seven GR-lncRNA prognostic signature. Using this model, BC patients were distinctly separated into high- and low-risk groups, demonstrating significantly different overall survival rates. The model’s stability and strong predictive capabilities were rigorously confirmed across the TCGA training, validation, and entire cohorts, indicating its robust potential for clinical application. The individual roles of these seven lncRNAs in BC, namely AL390728.6, AL355353.1, AC011477.2, AL357033.4, AC073210.3, AC105942.1, and FRMD6-AS2. AL390728.6, significantly upregulated in tumors related to normal samples. Although AL390728.6 remains largely uncharacterized in the literature, its consistent co-expression with GRGs hints at a potential role in hormonal and immunoregulatory signaling. Glucocorticoid signaling has been reported to exert context-dependent effects in cancers, including immune modulation, induction of apoptosis, and suppression of inflammatory responses (16,17). It is plausible that AL390728.6 enhances the anti-inflammatory or immune-sensitizing effects of glucocorticoids in BC, thereby contributing to improved immune surveillance and tumor control. Besides, Tang et al. have reported that AL390728.6 was predicted to regulate the hypoxic response during tumorigenesis via a hypoxia-miRNA-mRNA axis in HCC (18). AL390728.6 also emerged in acute myeloid leukemia studies, where it was one of ten lncRNAs incorporated into an immune-related prognostic signature (19). AL355353.1 also showed higher expression within tumors related to normal samples. Recent research demonstrated AL355353.1’s significant upregulation in BC tumor tissues and exosomal urine samples, supporting its clinical relevance and potential as a non-invasive biomarker (20). AC011477.2, another favorable prognostic marker, was recognized as one of five protective lncRNAs in a cuproptosis-associated prognostic signature in lung adenocarcinoma (21). AL357033.4, downregulated in tumors but associated with improved survival. Recent work by Zhao et al. established a redox-related lncRNA prognostic signature in BC, and AL357033.4 emerged as one of eight key lncRNAs incorporated into a robust risk scoring model (22). AC073210.3 and AC105942.1, which are primarily downregulated in tumors, were associated with poor prognosis when highly expressed, suggesting a paradoxical oncogenic role. Such lncRNAs may promote tumorigenesis by interfering with GRGs-mediated immunosuppressive signaling or by activating pro-tumor metabolic pathways such as glycolysis (23). FRMD6-AS2, the most notable risk-associated lncRNA, has highlighted its tumor-suppressive role across multiple cancer types. In prostate adenocarcinoma, FRMD6-AS2 was identified as part of an immune-related prognostic signature, where elevated expression of it correlated with favorable survival outcomes and an immune-active tumor microenvironment (24). Accumulating evidence from recent studies consistently highlights the significant prognostic value of lncRNA signatures across various cancer types, providing valuable tools for patient risk stratification (25,26). Our work extends this understanding by specifically focusing on the under-explored realm of GR-lncRNAs in BC.

Furthermore, our detailed immune infiltration analysis unveiled profound contrasts in the immune landscape across high- and low-risk BC patients. High-risk patients exhibited significantly higher ESTIMATE, Immune, and Stromal scores, alongside increased expression of ICP genes and elevated TIDE scores. This pattern suggests a more active, yet often dysfunctional and immunosuppressive, tumor microenvironment in high-risk tumors, indicative of enhanced immune evasion mechanisms and potential resistance to immunotherapy (27). The significantly lower IPS within the high-risk group further supports this notion, since a higher IPS generally correlates with improved response to ICP blockade (28). Conversely, CIBERSORT analysis revealed a higher infiltration of certain beneficial immune cell types, such as activated dendritic cells, follicular helper T cells, and CD8⁺ T cells, within the low-risk group. These cell types are well known for promoting anti-tumor immunity and facilitating immune-mediated tumor eradication (29-31). These findings suggest that individuals classified as low-risk may derive greater benefit from combination therapies designed to overcome immune suppression or enhance the efficacy of ICP inhibitors.

While our bioinformatics analysis is comprehensive, it relies on publicly available TCGA data, necessitating external validation in independent clinical cohorts with diverse ethnic backgrounds to confirm generalizability. In addition, treatment-related information, particularly regarding the administration of NAC in muscle-invasive bladder cancer (MIBC) patients, was not uniformly available in TCGA, which represents an important limitation of our study. Moreover, because the TCGA dataset does not specify whether samples were derived from TURBT or cystectomy, we could not assess the surgical source of the specimens. Nevertheless, if validated on TURBT material, our model may have important translational relevance in non-muscle-invasive bladder cancer (NMIBC), potentially identifying high-risk patients who might benefit from early cystectomy. Furthermore, treatment data, including BCG response, were not available in TCGA, limiting our ability to evaluate whether the GR-lncRNA signature could predict BCG efficacy. Rigorous in vitro and in vivo experimental validation is also crucial to fully elucidate the functional roles and precise molecular mechanisms of the identified GR-lncRNAs in BC. Future research should prioritize prospective clinical trials to confirm the prognostic and predictive power of this GR-lncRNA signature in real-world settings. Further mechanistic studies are vital for uncovering the exact pathways through which these lncRNAs impact glucocorticoid signaling, immune responses, and BC pathogenesis, setting the stage for novel therapeutic targets and innovative treatment strategies.

Conclusions

To summarize, two separate molecular subtypes of BC based on GR-lncRNA expression were recognized in this study and we have successfully constructed a novel 7 GR-lncRNA prognostic risk model. This model effectively stratifies BC patients into high- and low-risk groups with significantly different overall survival rates and distinct molecular and immunological characteristics. Our comprehensive analysis further revealed the intricate interplay between GR-lncRNAs, immune cell infiltration, TMB, and differential drug sensitivities, offering valuable insights into the profound heterogeneity of BC. These observations underscore the critical function of GR-lncRNAs in BC pathogenesis and provide potential biomarkers for improved prognostic assessment and the progression of more precise, individualized treatment approaches for BC patients.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tau.amegroups.com/article/view/10.21037/tau-2025-528/rc

Peer Review File: Available at https://tau.amegroups.com/article/view/10.21037/tau-2025-528/prf

Funding: This work was supported by the Zhejiang Provincial Medical and Health Science and Technology Planning Project (No. 2025KY1973).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tau.amegroups.com/article/view/10.21037/tau-2025-528/coif). L.Y. reports receiving funding from the Zhejiang Provincial Medical and Health Science and Technology Planning Project (No. 2025KY1973). The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. van der Heijden AG, Bruins HM, Carrion A, et al. European Association of Urology Guidelines on Muscle-invasive and Metastatic Bladder Cancer: Summary of the 2025 Guidelines. Eur Urol 2025;87:582-600. [Crossref] [PubMed]
2. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
3. Powles T, Bellmunt J, Comperat E, et al. ESMO Clinical Practice Guideline interim update on first-line therapy in advanced urothelial carcinoma. Ann Oncol 2024;35:485-90. [Crossref] [PubMed]
4. Yang J, Ariel F, Wang D. Plant long non-coding RNAs: biologically relevant and mechanistically intriguing. J Exp Bot 2023;74:2364-73. [Crossref] [PubMed]
5. Zou Y, Chen B. Long non-coding RNA HCP5 in cancer. Clin Chim Acta 2021;512:33-9. [Crossref] [PubMed]
6. Liu Q. The emerging roles of exosomal long non-coding RNAs in bladder cancer. J Cell Mol Med 2022;26:966-76. [Crossref] [PubMed]
7. Chen X, Xie R, Gu P, et al. Long Noncoding RNA LBCS Inhibits Self-Renewal and Chemoresistance of Bladder Cancer Stem Cells through Epigenetic Silencing of SOX2. Clin Cancer Res 2019;25:1389-403. [Crossref] [PubMed]
8. Rahmani F, Safavi P, Fathollahpour A, et al. The interplay between non-coding RNAs and Wnt/ß-catenin signaling pathway in urinary tract cancers: from tumorigenesis to metastasis. EXCLI J 2022;21:1273-84. [Crossref] [PubMed]
9. Chen JB, Zhang M, Zhang XL, et al. Glucocorticoid-Inducible Kinase 2 Promotes Bladder Cancer Cell Proliferation, Migration and Invasion by Enhancing β-catenin/c-Myc Signaling Pathway. J Cancer 2018;9:4774-82. [Crossref] [PubMed]
10. Rezaei S, Nikpanjeh N, Rezaee A, et al. PI3K/Akt signaling in urological cancers: Tumorigenesis function, therapeutic potential, and therapy response regulation. Eur J Pharmacol 2023;955:175909. [Crossref] [PubMed]
11. Wang C, Li K, Huang R, et al. Urine proteomics-based analysis identifies CHI3L1 as an immune marker and potential therapeutic target for bladder cancer. BMC Cancer 2025;25:271. [Crossref] [PubMed]
12. Hou J, Liang S, Xie Z, et al. An immune-related lncRNA model for predicting prognosis, immune landscape and chemotherapeutic response in bladder cancer. Sci Rep 2022;12:3225. [Crossref] [PubMed]
13. Thorsson V, Gibbs DL, Brown SD, et al. The Immune Landscape of Cancer. Immunity 2018;48:812-830.e14. [Crossref] [PubMed]
14. Qin Y, Zu X, Li Y, et al. A cancer-associated fibroblast subtypes-based signature enables the evaluation of immunotherapy response and prognosis in bladder cancer. iScience 2023;26:107722. [Crossref] [PubMed]
15. Kettunen K, Mathlin J, Lamminen T, et al. Profiling steroid hormone landscape of bladder cancer reveals depletion of intratumoural androgens to castration levels: a cross-sectional study. EBioMedicine 2024;108:105359. [Crossref] [PubMed]
16. Oakley RH, Cidlowski JA. The biology of the glucocorticoid receptor: new signaling mechanisms in health and disease. J Allergy Clin Immunol 2013;132:1033-44. [Crossref] [PubMed]
17. Khadka S, Druffner SR, Duncan BC, et al. Glucocorticoid regulation of cancer development and progression. Front Endocrinol (Lausanne) 2023;14:1161768. [Crossref] [PubMed]
18. Tang Y, Zhang H, Chen L, et al. Identification of Hypoxia-Related Prognostic Signature and Competing Endogenous RNA Regulatory Axes in Hepatocellular Carcinoma. Int J Mol Sci 2022;23:13590. [Crossref] [PubMed]
19. Qin L, Li B, Wang S, et al. Construction of an immune-related prognostic signature and lncRNA-miRNA-mRNA ceRNA network in acute myeloid leukemia. J Leukoc Biol 2024;116:146-65. [Crossref] [PubMed]
20. Tang D, Li Y, Tang Y, et al. Recognition of Glycometabolism-Associated lncRNAs as Prognosis Markers for Bladder Cancer by an Innovative Prediction Model. Front Genet 2022;13:918705. [Crossref] [PubMed]
21. Di H, Zhao J, Zhu X, et al. A novel prognostic signature for lung adenocarcinoma based on cuproptosis-related lncRNAs: A Review. Medicine (Baltimore) 2022;101:e31924. [Crossref] [PubMed]
22. Zhao F, Xie H, Guan Y, et al. A redox-related lncRNA signature in bladder cancer. Sci Rep 2024;14:28323. [Crossref] [PubMed]
23. Liu C, Li H, Chu F, et al. Long non coding RNAs: Key regulators involved in metabolic reprogramming in cancer Oncol Rep 2021;45:54. (Review). [Crossref] [PubMed]
24. Liang L, Xia W, Yao L, et al. Long non-coding RNA profile study identifies an immune-related lncRNA prognostic signature for prostate adenocarcinoma. Int Immunopharmacol 2021;101:108267. [Crossref] [PubMed]
25. Yu Z, Lu B, Gao H, et al. A New Prognostic Signature Constructed with Necroptosis-Related lncRNA in Bladder Cancer. J Oncol 2022;2022:5643496. [Crossref] [PubMed]
26. Wu L, Chen W, Cao Y, et al. A novel cuproptosis-related lncRNAs signature predicts prognosis in bladder cancer. Aging (Albany NY) 2023;15:6445-66. [Crossref] [PubMed]
27. Cao R, Yuan L, Ma B, et al. Tumour microenvironment (TME) characterization identified prognosis and immunotherapy response in muscle-invasive bladder cancer (MIBC). Cancer Immunol Immunother 2021;70:1-18. [Crossref] [PubMed]
28. Zhang X, Zhang Y, Zhao L, et al. Exploitation of tumor antigens and construction of immune subtype classifier for mRNA vaccine development in bladder cancer. Front Immunol 2022;13:1014638. [Crossref] [PubMed]
29. Wang Y, Xiang Y, Xin VW, et al. Dendritic cell biology and its role in tumor immunotherapy. J Hematol Oncol 2020;13:107. [Crossref] [PubMed]
30. Sun L, Su Y, Jiao A, et al. T cells in health and disease. Signal Transduct Target Ther 2023;8:235. [Crossref] [PubMed]
31. Niogret J, Berger H, Rebe C, et al. Follicular helper-T cells restore CD8(+)-dependent antitumor immunity and anti-PD-L1/PD-1 efficacy. J Immunother Cancer 2021;9:e002157. [Crossref] [PubMed]