Leveraging hypoxia-related genes signature for predicting the prognosis of bladder cancer

Introduction

Bladder cancer (BC), recognized as the tenth most prevalent malignancy globally, is responsible for an estimated 573,000 new cases and 213,000 deaths annually, representing 3.0% and 2.1% of the total cancer incidence and mortality rates, respectively (1). BC is classified into non-muscle-invasive BC (NMIBC) and muscle-invasive BC (MIBC) based on tumor infiltration depth, with 75% to 80% of cases initially presenting as NMIBC. Although new therapeutic methods such as minimally invasive techniques, immune checkpoint inhibitors (ICIs), and targeted therapies have been introduced, there is still no significant improvement on the 5-year survival rate, primarily due to its propensity for recurrence, infiltration, and metastasis (2). In addition, the accuracy of the existing clinical staging system, which is primarily based on pathologic findings, remains suboptimal in predicting the prognosis of BC patients and meeting clinical requirements (3). Current biomarkers, while offering some insights into tumor biology, often lack sufficient sensitivity and specificity to accurately reflect the biological heterogeneity of tumors or predict high-risk recurrence, metastasis, and treatment outcomes. These limitations highlight the urgent need for more robust molecular biomarkers to improve the precision of BC staging and prognostic predictions (4). Therefore, advancing research on biomarkers with greater predictive power remains a central focus in BC study.

Accumulating research indicates that hypoxia is associated with diseases such as cancer, diabetes, and inflammation (5). The hypoxia phenotype refers to a series of biological characteristics and physiological responses of cells in a low-oxygen environment. This adaptation enables cells to maintain their physiological functions by enhancing anaerobic metabolism, promoting angiogenesis, and regulating cell proliferation and apoptosis. Moreover, the hypoxia phenotype plays a crucial role in tumor growth, progression, and metastasis. The rapid growth and abnormal proliferation of tumor cells necessitate neoangiogenesis to secure additional oxygen and nutrients, fulfilling their growth requirements. However, abnormal neovascularization in the tumor facilitates the formation of highly permeable, dilated, and gyriform occluded areas in the tumor microenvironment (TME), leading to reduced blood perfusion, and consequently causing significant hypoxia in tumor cells (6). This alteration further promotes tumor plasticity and heterogeneity, leading to a more aggressive and metastatic phenotype. Hypoxia could also affect the immune system by modulating TME, thereby suppressing anti-tumor therapeutic efficacy. It has been reported that hypoxia could affect the genetic instability and malignant progression of MIBC and diminish its responsiveness to radiation therapy.

The development of a hypoxic phenotype is a multifaceted physiological process that involves the regulation of multiple signaling pathways and molecular mechanisms. Several studies have reported that hypoxia inducible factor (HIF), one of the best-known hypoxia-related genes (HRGs), serves as a crucial transcriptional regulator for over 100 genes that facilitates adaptation to hypoxic environments. Tumor cells could increase the transcriptional activity of HIF to adapt to hypoxia, influencing BC staging, grading, and overall survival (OS) negatively (7,8). Besides the HIF signaling pathway, additional signaling molecules like mTOR, PI3K/Akt, and MAPK are also involved in regulating the hypoxic phenotype. While the role of HRGs has already been described in various tumor types, research on HRGs in BC remains limited. Thus, a comprehensive investigation into the biological role of HRGs in BC development is imperative and carries significant clinical relevance and implications.

This study was conducted in a combination of bioinformatics, machine learning, and single-cell RNA sequencing (scRNA-seq) analysis to comprehensively identify the potential functions of differentially expressed HRGs, aiming to improve the understanding of the pathogenesis of BC. Leveraging the GeneCards database, which integrates extensive gene annotations and diverse biological and clinical datasets, we identified HRGs with strong relevance to hypoxic processes. This comprehensive resource provided valuable support for analyzing the relationship between genes and biological processes related to hypoxia. We selected the hub prognostic HRGs to construct a hypoxia-related risk model for predicting prognosis, informing pharmaceutical strategies, and elucidating the immune microenvironment landscape for BC patients from the perspective of the hypoxic phenotype. We present this article in accordance with the TRIPOD reporting checklist (available at https://tau.amegroups.com/article/view/10.21037/tau-2025-118/rc).

Methods

Data download

The bladder urothelial carcinoma (BLCA) dataset was downloaded from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) via the TCGAbiolinks package (9) and analyzed as a test set. After excluding samples lacking clinical information, the count sequencing data of 410 BC samples with prognostic OS clinical information were obtained. The count sequencing data of BC were also standardized to fragments per kilobase per million (FPKM) format. The clinical data corresponding to the samples were obtained from the UCSC Xena database (10) (http://genome.ucsc.edu). Two BC-related datasets, GSE13507 (11) and GSE48075 (12), respectively containing 165 and 73 BC samples with prognostic OS clinical information, were downloaded from the Gene Expression Omnibus (GEO) database via the R package GEOquery (13). The single-cell dataset GSE135337 (14) was also downloaded which contains seven BC samples and one control sample. Besides, the HRGs were collected through the GeneCards database (15) (https://www.genecards.org/) which provides comprehensive information on human genes. Using the word “Hypoxia” as the search keyword, we screened genes with correlation value greater than 5, and a total of 126 HRGs were obtained. Meanwhile, the microsatellite instability (MSI) data were downloaded from the cBioPortal database (16) (https://www.cbioportal.org/). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Hypoxia score definition

Single-sample gene set enrichment analysis (ssGSEA) algorithm could quantify the relative abundance of each gene expression in single sample. We calculated the hypoxia score for each sample in the TCGA-BLCA dataset based on the expression matrix of HRGs using the gene set variation analysis (GSVA) package (17) via the ssGSEA algorithm. After plotting Kaplan-Meier (KM) survival curves with the TCGA-BLCA dataset, we characterized the relationship between HRGs features and OS for each sample to determine the optimal hypoxia score. Then the BC patients were divided into high-scoring and low-scoring groups based on hypoxia score, and P value <0.05 was considered to be statistically significant.

Prognostic model construction

The differentially expressed genes (DEGs) between different groups were identified using the R package “limma” (18). The screening criteria for DEGs were |log fold change (FC)| ≥0.5 and adjusted (adj.) P<0.05. In order to construct the prognostic model of BC, univariate Cox regression was performed to obtain HRGs associated with prognosis, and then the genes with P<0.05 were screened for the least absolute shrinkage and selection operator (LASSO) regression and multivariate Cox regression. LASSO is a regularized regression approach commonly used for feature selection, using 10-fold cross-validation based on the R package “glmnet” (12). Subsequently, the coefficients of each gene in the prognostic model were extracted to calculate the risk score for each sample. Additionally, the function “surv_cutpoint” of the survminer package was used to define the optimal cut-off expression for the best grouping. The risk score was calculated using the following formula:

Riskscore=∑iCoefficient(hubgenei)∗mRNAexpression(hubgenei)[1]

GSEA

GSEA was used to assess gene distribution trends from a predefined gene set in the gene table to determine their contribution to the phenotype (19). GSEA was performed using the R package clusterprofiler. The minimum gene set size was set at 10 genes, and the maximum size at 500 genes. The P value correction was performed using Benjamini-Hochberg (BH) method. Enrichment was considered significant for P<0.05 and false discovery rate (FDR) value (q value) <0.25.

Somatic mutation (SM) analysis

To delineate the SM landscape in BC, we downloaded “Masked Somatic Mutation” subtype of SM data from the TCGA-BLCA dataset, processed the data using VarScan software, and subsequently analyzed and visualized it with the R package “Maftools” (20). Furthermore, the “OncogenicPathways” and “PlotOncogenicPathways” functions were utilized for the biological analysis to identify the oncogenic pathways affected by the mutations. Additionally, the “drugInteractions” function was applied to explore the interactions between drugs and mutations, aiming to identify mutations associated with specific drug sensitivities and resistances.

Immune infiltration analysis

The ssGSEA algorithm was used to quantify the relative abundance of each immunocytes infiltrate. Each infiltrating immune cell type was labeled. Enrichment scores calculated by ssGSEA were used to represent the relative abundance of immunocytes in each sample. The ggplot2 package was used to visualize data graphs.

Drug sensitivity analysis potential drug prediction

The potentially druggable genes for samples in high- and low-risk groups were filtered for druggability using the Drug Gene Interaction database (DGIdb; http://dgidb.genome.wustl.edu/). Besides, the drug sensitivity analysis for hub genes was performed using the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) (21), the Cancer Cell Line Encyclopedia (CCLE) database (https://portals.broadinstitute.org/ccle/about) (22), and the CellMiner database (https://discover.nci.nih.gov/cellminer/home.do) (22). Additionally, connectivity map (Cmap; https://clue.io/) was used to identify small molecule compounds targeting hub genes to explore potential drugs. The GDSC database serves as a comprehensive repository for drug sensitivity and molecular response markers in diverse cancer cell types. The CCLE database determines drug sensitivity using pharmacological profiles of 24 anticancer drugs across 479 cell lines. Moreover, the CellMiner database offers raw data linking experimental outcomes to hub genes for comprehensive gene analysis.

Molecular docking

The molecular structures of candidate drugs and herbal monomers, as solasonine and rhein, were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The crystal structures of hub genes were acquired from the Protein Data Bank (PDB) database (https://www.rcsb.org). The protein structures were then prepared by adding hydrogen atoms and merging nonpolar hydrogen atoms using AutoDock Tools (ADT; version 1.5.6) (23) and saved as PDBQT files. The docking pocket of protein was defined using the grid box module. CB-Dock2 (https://cadd.labshare.cn/cb-dock2/php/index.php) and AutoDock Vina (version 1.1.2) (24) were utilized to perform ligand-receptor docking. After evaluating the binding energy (affinity), the conformation with the best affinity was visualized in PyMOL (version 2.5.2) software. The 2D diagram of the protein-ligand was visualized using the Discovery Studio server.

Cell culture

The human ureteral epithelial immortalized cell line SV-HUC-1, and the human BC cell lines T24 and RT-112 were purchased from the Shanghai Institutes for Biological Sciences (Shanghai, China). Ten percent fetal bovine serum was used in Ham’s F-12K medium for SV-HUC-1, and in RPMI-1640 medium for T24 and RT-112. A total of three cell lines were supplemented with 100 U/mL penicillin and 100 µg/mL streptomycin, and cultured in a cell culture incubator at 37 °C with 5% CO₂.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR) analysis

The qRT-PCR was used to detect the messenger RNA (mRNA) expression levels of candidate genes. Total cellular RNA was first extracted using Trizol (Invitrogen, Carlsbad, CA, USA), and then reverse-transcribed to cDNA according to the instructions of the cDNA First Strand Synthesis Kit (TransGen Biotech, Beijing, China; AT341). Finally, the cellular mRNA expression levels of candidate genes were detected by qRT-PCR using AceQ^TM Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China; Q511-03), with β-actin as an internal reference. The primers sequences are as follows: JUN forward, 5'-TCCAAGTGCCGAAAAAGGAAG-3', JUN reverse, 3'-CGAGTTCTGAGCTTTCAAGGT-5', MYC forward, 5'-GGCTCCTGGCAAAAGGTCA-3', MYC reverse, 3'-CTGCGTAGTTGTGCTGATGT-5', EGFR forward, 5'-AGGCACGAGTAACAAGCTCAC-3', EGFR reverse, 3'-ATGAGGACATAACCAGCCACC-5', SLC2A1 forward, 5'-TCTGGCATCAACGCTGTCTTC-3', SLC2A1 reverse, 3'-CGATACCGGAGCCAATGGT-5', β-actin forward, 5'-CATGTACGTTGCTATCCAGGC-3', β-actin reverse, 3'-CTCCTTAATGTCACGCACGAT-5'.

Quality control and analysis for single-cell sequencing datasets

The “CreateSeuratObject” function in the R package Seurat v4.0 (25) was used to import the count matrices for BC and control samples in the scRNA-seq dataset GSE135337 and created as Seurat objects. Parameters were set to filter for genes expressed in at least three cells and to ensure that each cell expresses at least 200 genes. As for a higher proportion of mitochondrial genes within the cellular genome often suggests cellular stress rather than a steady state, cells with mitochondrial gene content exceeding 10% were excluded.

All samples in the GSE135337 dataset underwent normalization for sequencing depth using the “SCTransform” function. Then, the principal component analysis (PCA) was applied to identify the principal component (PC), and the “Elbowplot” function was used to visualize the P value distribution. The first 15 PCs were selected for dimensionality reduction using the t-distributed stochastic neighbor embedding (tSNE) algorithm. The “FindNeighbors” function, utilizing the default parameters and the 15 PCs dimensions, was used to construct the k-nearest neighbors of Euclidean distances in PCA space. The “FindClusters” and “clustree” functions were used at a resolution of 0.3 to classify the cells into different clusters, and the “RunTsne” function was used to perform visualization. Finally, leveraging the reference dataset BlueprintEncodeData, we annotated the scRNA-seq dataset GSE135337 with cell types via “SingleR” function in the R package SingleR (26) for cell type identification.

Pseudo-time analysis

The R package “Monocle” (27) was utilized for pseudo-time analysis (28), which enables the arrangement of each cell along its corresponding developmental trajectory, determined by the temporal gene expression of each cell. Based on the gene expression status, the sample could be classified into multiple cell populations under differentiated states, leading to generate an intuitive lineage development dendrography for cell differentiation and development trajectory prediction.

Statistical analysis

Data processing and analyses in this study were conducted using R software (version 4.2.0). Continuous variables were presented as mean ± standard deviation (SD). The Wilcoxon Rank Sum Test was used for comparison between the two groups. Correlation coefficients were calculated via a Spearman correlation analysis if not specified. A P value <0.05 was considered as the criterion for a significant statistical difference.

Results

Flow chart and related data information

The workflow of this study design is illustrated in Figure 1. The specific information of the datasets involved is shown in Table 1, and the information of HRGs involved is shown in Table S1.

Figure 1 The workflow of this study design. BLCA, bladder urothelial carcinoma; DEGs, differentially expressed genes; GSEA, gene set enrichment analysis; HRGs, hypoxia-related genes; KM, Kaplan-Meier; LASSO, least absolute shrinkage and selection operator; MSI, microsatellite instability; PCA, principal component analysis; ROC, receiver operating characteristic; SNP, single-nucleotide polymorphism; TCGA, The Cancer Genome Atlas.

Hypoxia score construction

The diagnostic performance and prognostic value of hypoxia score for BC were evaluated through receiver operating characteristic (ROC) and KM survival analyses, respectively, revealing a better diagnostic capability [area under the curve (AUC) =0.731, Figure 2A] and significant prognostic value (P<0.05, Figure 2B). Following differential expression analysis of TCGA-BLCA dataset, a total of 1094 DEGs were identified between the high-scoring and low-scoring groups, including 764 upregulated and 330 downregulated genes (Figure 2C). The expression of the top 10 positively and negatively regulated genes was displayed in the heatmap (Figure 2D). The hypoxia score clearly differentiated between high-scoring and low-scoring groups as evidenced by the PCA analysis (Figure 2E). Furthermore, significant differences in MSI between the groups were noted, with the low-scoring group exhibiting a higher MSI score, indicative of increased microsatellite mutations in the low-scoring group (Figure 2F).

Figure 2 Hypoxia score construction and analysis. (A) ROC curve of hypoxia scores in the TCGA-BLCA dataset. (B) KM curve of hypoxia scores in the TCGA-BLCA dataset. (C) Volcano plots of DEGs between the high-scoring group and low-scoring group in the TCGA-BLCA dataset. (D) Heatmap of the top 10 positively and negatively regulated DEGs. (E) The PCA plots of the high-scoring group in blue and low-scoring group in red. (F) Comparison plots of MSI scores in the high-scoring group and low-scoring group of the TCGA-BLCA dataset. ***, P<0.001. AUC, area under the curve; BLCA, bladder urothelial carcinoma; DEGs, differentially expressed genes; HRG, hypoxia-related gene; KM, Kaplan-Meier; MSI, microsatellite instability; PC, principal component; PCA, principal component analysis; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Evaluation prognosis prediction ability of the hypoxia-related risk model

Through LASSO Cox regression analysis on the TCGA-BLCA dataset, prognostic genes were selected (Figure 3A), and 14 genes were further identified based on the minimum partial likelihood deviation (Figure 3B). By intersecting with DEGs that identified according to the hypoxia score, four hub genes were finally obtained, including JUN, MYC, EGFR, and SLC2A1, to construct the prognostic risk model (Figure 3C). The predictive ability of this model was further evaluated using time-dependent ROC curves, suggesting good predictive performance for OS of BC patients, with AUC values of 0.612 at 1-year, 0.627 at 3-year, and 0.632 at 5-year (Figure 3D). Prognostic survival KM curves were also plotted, revealing statistically significant predictive ability for prognostic survival in BC patients (P<0.01, Figure 3E). Additionally, univariate Cox regression analysis was conducted to assess whether the risk model and clinical variables served as independent prognostic factors for BC using the TCGA-BLCA cohort. The analysis revealed that the hazard ratio (HR) of the risk score was 2.57 [95% confidence interval (CI): 1.68–3.93] (P<0.01, Figure 3F). Moreover, a nomogram plot was generated to predict 1-, 3-, and 5-year survival in BC patients by using age, gender, stage and risk score (Figure 3G), with calibration curves demonstrating good calibration and discrimination of the nomogram (Figure 3H).

Figure 3 Evaluation prognosis prediction ability of the hypoxia-related risk model. (A) LASSO regression analysis and (B) LASSO variable trajectory diagram of the TCGA-BLCA dataset. (C) Venn diagram between genes identified using Lasso algorithm and DEGs identified using hypoxia score. (D) The time-dependent ROC analysis for risk score in the TCGA-BLCA cohort. (E) KM curves for OS of risk score in the TCGA-BLCA cohort. (F) Forest plots of univariate Cox regression analysis involving risk score and clinical factors. (G) The nomogram of risk score and clinical factors. (H) The calibration curve of nomogram. AUC, area under the curve; BLCA, bladder urothelial carcinoma; CI, confidence interval; DEGs, differentially expressed genes; HR, hazard ratio; KM, Kaplan-Meier; LASSO, least absolute shrinkage and selection operator; OS, overall survival; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Validation of the prognostic risk model

We further validated the predictive performance of the prognostic risk model using the GSE13507 and GSE48075 datasets. The hub genes were all upregulated in the high-risk group in both datasets (Figure S1). The time-dependent ROC curves of the risk score showed that the 1-, 3-, and 5-year OS AUCs were 0.675, 0.609, and 0.579 for the GSE13507 dataset (Figure 4A), and 0.674, 0.651, and 0.6 for the GSE48075 dataset (Figure 4B), indicating an excellent predictive ability. Next, prognostic survival KM curves suggest that the risk model was statistically significant in predicting prognostic survival of BC patients (P<0.01, Figure 4C,4D). Furthermore, the nomogram based on the risk score and risk factors also demonstrated good predictive effect (Figure 4E,4F), with good calibration (Figure 4G,4H).

Figure 4 Validation of the prognostic risk model. Time-dependent ROC analysis for risk score in the (A) GSE13507 and (B) GSE48075 datasets. KM curves for OS of risk score in the (C) GSE13507 and (D) GSE48075 datasets. The nomogram of risk score and clinical factors for the (E) GSE13507 and (F) GSE48075 datasets. The calibration curve of nomogram for the (G) GSE13507 and (H) GSE48075 datasets. AUC, area under the curve; KM, Kaplan-Meier; OS, overall survival; ROC, receiver operating characteristic.

Analysis of clinical correlation and GSEA for risk score

The relationship between risk score and clinical characteristics in BC patients was analyzed. The results showed that the proportions of older patients and female patients were higher in the high-risk group than in the low-risk group (Figure 5A,5B). Besides, risk score differed between tumor stages, with stage 4 generally having a higher risk score and stage 1 having a lower risk score (Figure 5C). Moreover, the results of GSEA revealed the significantly enriched P53 downstream pathway, PI3K/AKT signaling pathway, signaling by EGFR, and WNT signaling pathway associated with the risk score (Figure 5D, table available at https://cdn.amegroups.cn/static/public/tau-2025-118-1.pdf).

Figure 5 Clinical correlation and GSEA in the high- and low-risk groups. (A) The relationship of risk score with age and (B) gender in BC patients. (C) The boxplot of risk score in BC patients with different tumor stages. (D) GSEA of risk score in BC patients. *, P<0.05; **, P<0.01; ***, P<0.001; ns, no statistical significance. BC, bladder cancer; GSEA, gene set enrichment analysis; NES, normalized enrichment score.

Analysis of tumor mutation characteristics for risk score

The tumor mutation burden characteristics in high- and low-risk groups showed that TP53, TTN, KMT2D, KDM6A, ARID1A, and MUC16 had a high mutation frequency in BC patients (Figure 6A,6B). TP53 has the highest mutation frequency, with 68% in the high-risk group, and 45% in the low-risk group. Biological function alterations caused by mutations in the high-risk group were mainly focused on signaling pathways such as RTK-RAS and Hippo (Figure 6C), while those caused by mutations in the low-risk group were mainly focused on signaling pathways such as RTK-RAS and NOTCH (Figure 6D). The potentially druggable genes were predicted using the DGIdb database, revealing that the genes enriched in the druggable genome in the high-risk group were EP300, ERBB2, FAT4, HMCN1, and MUC16 (Figure 6E), and in the low-risk group were ATM, BIRC6, FAT4, HMCN1, and MUC16 (Figure 6F).

Figure 6 Mutation characteristics in high- and low-risk groups. Landscape of mutated genes in the (A) high-risk group samples and (B) low-risk group samples. Biological functions affected by mutations in the (C) high-risk group samples and (D) low-risk group samples. Classification of potentially druggable genes in the (E) high-risk group samples and (F) low-risk group samples. The top 5 genes are shown in parentheses after each classification, and all genes are shown if there are fewer than 5. The X-axis represents the number of genes in the druggable gene classification. TMB, tumor mutational burden.

Immune infiltration assessment

The infiltration of 28 immune cells between high- and low-risk group samples in the TCGA-BLCA dataset was estimated by ssGSEA algorithm (Figure 7A). The immune infiltration score of different immune cells was calculated, which showed that the high-risk group exhibited increased infiltration of multiple immune cells, including regulatory T cells (Tregs), myeloid-derived suppressor cells (MDSCs), neutrophils, and others (P<0.05, Figure 7B).

Figure 7 Immune infiltration assessment. (A) Heatmap of 28 immune cells infiltration abundance in TCGA-BLCA samples in the high- and low-risk groups. (B) Comparison of immune cell infiltration between the high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001; ns, no statistical significance. BLCA, bladder urothelial carcinoma; TCGA, The Cancer Genome Atlas.

Hub genes drug sensitivity analysis

To investigate the relationship between the drug sensitivity and the mRNA expression of the four hub genes in the risk model, a ridge regression model was constructed to predict the half-maximal drug inhibitory concentration (IC₅₀) of the commonly used antitumor drugs by pRRophetic algorithm, utilizing gene expression profiles in TCGA and cell line expression profiles in CellMiner, GDSC, and CCLE databases. The association between hub genes expression and drug sensitivity was visualized, revealing 718, 1,279, and 69 drugs for interaction with hub genes through the CellMiner (Figure 8A), GDSC (Figure 8B), and CCLE database (Figure 8C), respectively.

Figure 8 Hub genes drug sensitivity analysis. Evaluation of the antitumor drug sensitivity for hub genes via (A) CellMiner, (B) GDSC, and (C) CCLE databases. The interaction pairs with higher correlation were presented. CCLE, Cancer Cell Line Encyclopedia; FDR, false discovery rate; GDSC, Genomics of Drug Sensitivity in Cancer.

Prediction of potential drugs for hub genes

To identify potential drugs, we also explored the interaction between hub genes and small molecule compounds using the Cmap database. The drug-gene interaction network with the highest scores was presented, revealing that EGFR had high sensitivity with AG-490 (score =98.13), while MYC (score =37.22) and JUN (score =37.22) had high sensitivity with TWS-119 (Figure 9A). Both AG-490 and TWS-119 were associated with tumor progression. Subsequently, the corresponding molecular docking with potential drugs was visualized via CB-Dock2 (Figure 9B-9D).

Figure 9 Prediction of potential drugs for hub genes. (A) The drug-gene interaction network with the highest score for hub genes interacting with small molecule compounds. (B) Molecular docking maps of EGFR with AG-490, (C) MYC with TWS-119, and (D) JUN with TWS-119.

Molecular docking analysis of solasonine and rhein towards hub genes

The crystal structures of hub genes were acquired, including JUN (PDB code: 1JUN), EGFR (PDB code: 1M14), MYC (PDB code: 1A93), and SLC2A1 (PDB code: 6THA). Binding affinity was measured in kcal/mol, with a smaller affinity value indicating stronger binding. Generally, an affinity of <−7.00 kcal/mol represents a satisfactory binding strength between the ligand and receptor. The molecular docking analysis showed that both solasonine and rhein exhibited significant binding strengths with EGFR (−9.4 and −8.0 kcal/mol, Figure 10A,10B) and SLC2A1 (−10.5 and −9.5 kcal/mol, Figure 10C,10D), leading to stable complex formation. While both solasonine and rhein showed weak affinities to JUN (−6.2 and −5.4 kcal/mol) and MYC (−6.5 and −5.2 kcal/mol) (Figure S2).

Figure 10 Molecular docking of herbal monomers with hub genes. Molecular docking of herbal monomer solasonine with (A) EGFR and (B) SLC2A1. Molecular docking of herbal monomer rhein with (C) EGFR and (D) SLC2A1.

Expression validation of hub genes in vitro

The expression of JUN, EGFR, MYC, and SLC2A1 was identified in various urinary tract cancer cells using the CCLE database (Figure 11A). Among these, T24 and RT-112 cells exhibited relatively high expression levels of the four genes, making them suitable models for subsequent gene expression validation. We then compared their expression between normal and BC cells by performing qRT-PCR in SV-HUC-1, T24, and RT-112 cell lines. As result, JUN, EGFR, MYC, and SLC2A1 were significantly highly expressed in BC cells compared with normal urothelial cells (P<0.01, Figure 11B), indicating their positive correlation with BC development.

Figure 11 Expression validation of hub genes in vitro. (A) The mRNA expression of JUN, EGFR, MYC, SLC2A1 in 26 urinary tract cancer cells extracted from CCLE database. (B) The relative expression of hub genes in SV-HUC-1, T24, and RT-112 cell lines. **, P<0.01; ***, P<0.001. CCLE, Cancer Cell Line Encyclopedia; mRNA, messenger RNA.

Clustering and annotation of single cells

We processed the BC single-cell RNA-seq dataset GSE135337 for dimensionality reduction clustering using tSNE (Figure 12A). At a resolution of 0.3, the cells were classified into 12 independent clusters (Figure 12B). We further identified the cell clusters as a total of four cell types by cell annotation with R package SingleR, including epithelial cells, macrophages, fibroblasts, and endothelial cells (Figure 12C). The expression of cell marker genes EPCAM, KRT8, CD14, CSF1R, AIF1, DCN, PDPN, TAGLN, PECAM1, VWF, and CLDN5 indicated a relatively accurate cellular annotation (Figure 12D).

Figure 12 Clustering and annotation of single cells. (A) Clustering visualization of GSE135337 dataset by tSNE. (B) Classifying cells into 12 independent clusters by tSNE at a resolution of 0.3. (C) The clusters were annotated into four main cell types. (D) Dot plots of cellular marker genes for each cluster. tSNE, t-distributed stochastic neighbor embedding.

Pseudo-time analysis for malignant epithelial cells

We subsequently selected malignant epithelial cells for pseudo-time analysis to explore the developmental trajectories of different cell states (Figure 13A,13B). The findings showed that the developmental starting point was state2, and the developmental endpoint was state1. Dot plots and heat maps depicted the temporal expression of hub genes in different cell states. As shown in Figure 13C,13D, the expression of JUN varied greatly during the development of state2 to state1, and the expression of MYC was higher in state1.

Figure 13 Pseudo-time analysis for malignant epithelial cell. (A) Diagram of the temporal sequence malignant epithelial cell differentiation and development. The colors transition from dark to light to represent the proposed chronological order. (B) Developmental trajectory of malignant epithelial cell differentiation. (C) Dot plot and (D) heat map showing the temporal expression of hub genes.

Discussion

BC exhibits highly heterogeneous morphologically and clinically. Present diagnostic and therapeutic approaches face numerous challenges, including the limited accuracy for non-invasive diagnostic techniques, the inherent risk and complications associated with invasive diagnostic methods, and individual variability in treatment protocols. The complexity of treatment partly arises from the alteration of tumor cell biology within the hypoxic TME. Hypoxia, as a microenvironmental pressure prevalent in solid tumors, promotes tumor growth by activating specific hypoxia-responsive pathways and also facilitates immune evasion by modulating the tumor immune microenvironment (29). Thus, elucidating the biological role of hypoxic phenotype in BC development comprehensively is urgently necessary.

According to the ROC analysis and KM survival analysis, hypoxia score system showed a favorable performance in diagnostic and prognostic prediction for BC patients. By comparing the DEGs expressions between the high- and low-scoring groups, CXCL8 (also known as interleukin-8) was demonstrated to be the most highly expressed gene in the high-scoring group, and UPK1A was the highest in the low-scoring group. CXCL8 has been shown to activate endothelial cells (30), facilitating tumor cells to obtain more nutrients and oxygen for survival and proliferation, thus inducing tumor progression and chemotherapy resistance (31,32). Immunohistochemical analysis of UPK1A is a potential specific biomarker for distinguishing BC from other tumors. However, its sensitivity in MIBC decreased significantly, falling below 50% (33), which is consisted with our findings. Besides, S100A7 and S100A8, members of the S100 protein family, also exhibited significant upregulation in the high-scoring group. This upregulation is correlated with deeper muscle invasion in BC (34), as well as tumorigenic processes and immune evasion across various tumor types (35). The detrimental effect of hypoxia on immune cell function is considered essential for immune evasion mechanisms (36). Consistently, our results showed a decreased MSI in the high-scoring group. This suggested that functional hypoxia induced immune escape, which makes tumor cells more resistant to immune attack, consequently resulting in unsatisfactory effectiveness and prognosis.

By the intersection of the genes obtained from LASSO and previously obtained DEGs using hypoxia-scoring system, four hub genes, JUN, MYC, EGFR, and SLC2A1, were filtered. These hub genes were all involved in the high-hypoxia-scoring group, indicating an activation of them in response to hypoxia conditions. JUN and MYC, as oncogenic transcription factors, are crucial targets of the WNT and NOTCH pathways and also have binding sites in the promoter sequence of EGFR (37). SLC2A1 is the most common glucose transporter responsible for mediating glucose metabolism in tumor cells (38). Four hub genes have been found to modulate tumor metabolism, promote tumorigenesis, correlate with tumor immune invasion and immunotherapy efficacy, and might serve as rational oncotargets in cancer therapy (38-41). The prognostic risk model based on hub genes exhibited an excellent predictive ability in stratifying BC risk groups in the TCGA and GEO datasets, indicating the valuable prognostic implication of the genic cluster in BC. Our model offers a simple, efficient, and flexible method for prognostic prediction in BC patients, which might address the limitations of traditional classification methods in distinguishing patients with identical tumor-node-metastasis (TNM)/American Joint Cancer Committee (AJCC) stages and aid in devising personalized therapeutic strategies.

A higher proportion of elderly and female patients was observed in the high-risk group compared to the low-risk group. Additionally, the risk score varied across tumor stages, with patients at higher tumor stages generally exhibiting higher risk scores. The pathogenic role of hypoxia in senescence induction has recently been demonstrated. Hypoxia could drive senescence through multiple pathways, and the senescence could in tune exacerbate local tissue hypoxia (42). Estrogen exhibits oncogenic properties and facilitates neovascularization in various tumors under hypoxic conditions. Hypoxia and estrogen are interchangeable, both similarly regulating epithelial-endothelial cell interactions, inducing epithelial-mesenchymal transition (EMT), and then promoting tumor cell metastasis (43). GSEA showed a clear enrichment of P53 downstream pathway, PI3K/AKT signaling pathway, signaling by EGFR, and WNT signaling pathway, which were all involved in tumor progression and prognosis.

Tumor hypoxia could promote gene mutation (44). The analysis of mutation characteristics in the high- and low-risk groups revealed that TP53, TTN, KMT2D, KDM6A, ARID1A, and MUC16 had a high mutation frequency in BC. These frequently mutated genes might be critical genetic events that drive tumor progression or modulate responses to therapy. Among them, TP53 had the highest mutation frequency in both groups. Increasing evidence suggests that co-mutation of TP53 and RB1 within the hypoxic microenvironment correlates with a poor prognosis in BC patients and might render tumors resistant to ICIs (23,24). Additionally, TP53 also serves as a predictive biomarker for clinical outcomes and therapeutic response in MIBC patients (45). We also found that in the high-risk group mutations primarily affected RTK-RAS and Hippo pathways, while in the low-risk group, mainly focused on RTK-RAS and NOTCH pathways. Further druggable genes predicted in the high-risk group, such as EP300, ERBB2, and MUC16, were closely associated with the RTK-RAS signaling pathway, while FAT4 was involved in the Hippo signaling pathway. As for HMCN1, it might participate in multiple signaling pathways, but its specific role requires further investigation.

We also compared the tumor immune microenvironment in different risk groups, showing that the high-risk group exhibited increased infiltration of various immune cells, including Tregs, MDSCs, and neutrophils. Tregs, traditionally viewed as unilateral suppressors of antitumor immune responses, are associated with unfavorable outcomes in cancer patients. Activated Tregs could express immune checkpoint proteins, such as PD1, PDL1, and CTLA4, suppressing tumor-specific T-lymphocyte functions and establishing an immunosuppressive microenvironment that enhances BC aggressiveness (46,47). Therefore, high-risk patients might benefit more from ICIs therapy than low-risk patients. Meanwhile, MDSCs mediate tumor immune escape by inhibiting T-cell responses, thereby facilitating tumor growth and metastasis (48). Neutrophils are critical players at the forefront of immune responses, which could be regulated by hypoxia and generally exhibit a dual role in tumor development because of their plasticity (49). They have been implicated in promoting tumor proliferation by suppressing the immune system (50). Furthermore, higher neutrophil infiltration was reported to be observed in MIBC tissues compared to NMIBC, underscoring their association with worse outcomes (51).

To explore the potential of hub genes as therapeutic targets for BC, we assessed their interactions with common anticancer drugs. Multiple drugs were found to influence the stability or expression of these hub genes. By further identifing the highest scoring drug-gene interaction networks, we found the associations of EGFR with AG-490, as well as JUN and MYC with TWS-119. It has been demonstrated that AG-490, a type of EGFR inhibitor, could effectively inhibit tumor cell proliferation by limiting the expression of cell cycle protein D1 (52). TWS119, as a p-GSK-3β inhibitor and WNT pathway activator, not only mediates the proliferation and cytokine production of tumor-infiltrating lymphocytes (TILs), but also enhanced the proliferation and cytolytic activity of human γδT cells against cancer (53,54). However, their therapeutic efficacy in BC remains to be investigated. In addition, the capacity of both solasonine and rhein to interact with hub genes was evaluated. Solasonine and rhein have shown promising results as potential therapeutic options for BC in previous studies (55,56). As a result, the selective and significant interactions of these herbal monomers with targeted proteins were revealed, particularly EGFR and SLC2A1, suggesting their potential as promising candidates for targeted therapeutics in BC.

Single-cell sequencing, in contrast to conventional sequencing, not only offers more nuanced insights into tumor cellular heterogeneity, but also generates a comprehensive map of the TME at the single-cell level, including genomics, transcriptomics, epigenetics, and multi-omics (57). Our single-cell sequencing analysis showed that malignant epithelial cells constituted the predominant cell type in BC tissue. Consistent with our findings, 75% of BC patients were reported to test positive for EPCAM, which is thought to be significantly associated with OS in BC patients (58). We subsequently selected malignant epithelial cells for pseudo-time analysis to explore the developmental trajectories of different cell states. The results revealed significant dynamic changes in JUN expression along the developmental trajectory, while MYC expression was markedly upregulated at the terminal stage. The dynamic variation of JUN suggests its critical regulatory role in cell differentiation and tumor progression, reflecting the complexity of transcriptional regulatory networks and the responses of tumor cells to microenvironmental stimuli at different developmental stages. Meanwhile, the elevated MYC expression at the terminal stage highlights its role in maintaining the malignant phenotype and promoting late-stage tumor cell proliferation, further underscoring its potential as a therapeutic target. These findings indicate that JUN and MYC are pivotal in BC progression, offering insights into its etiology and hypoxia-related mechanisms. Further research into these hub genes might elucidate the etiology of BC and its relationship with hypoxia, providing new perspectives for future therapeutic strategies.

However, there are several limitations in this study. First, the prognostic risk model was primarily developed and validated using publicly available datasets, lacking validation with large-scale, multi-center clinical samples. Second, while drug sensitivity predictions identified potential therapeutic targets, these results require further validation through in vitro and in vivo experiments to assess their clinical applicability. Addressing these limitations in future studies will strengthen the reliability and translational potential of this hypoxia-related prognostic model for BC.

Conclusions

In summary, we developed a hypoxia-related prognostic risk model based on four hub HRGs, including JUN, MYC, EGFR, and SLC2A1. This model was validated as an independent predictor for BC patient prognosis, which offers insights into immune microenvironment and carcinogenesis mechanisms of BC. Analyzing mutation patterns, drug sensitivities, and the developmental trajectories of genes within this model would be helpful for understanding of carcinogenesis and refining therapeutic strategies.

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-118/rc

Data Sharing Statement: Available at https://tau.amegroups.com/article/view/10.21037/tau-2025-118/dss

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

Funding: This research was supported by the National Natural Science Foundation of China (to Y.D.) (No. 82204866), the Natural Science Foundation of Jiangsu Province (to Y.D.) (No. BK20231163), the Pengcheng Talent-Medical Young Reserve Talent Training Project (to Y.D.) (No. XWRCHT20220012), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (to Y.A.C.) (No. SJCX23_1399), and the Jiangsu Provincial TCM Science and Technology Development Plan Project (to Z.D.S.) (No. MS2023081).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tau.amegroups.com/article/view/10.21037/tau-2025-118/coif). Y.D. reports funding from the National Natural Science Foundation of China (No. 82204866), the Natural Science Foundation of Jiangsu Province (No. BK20231163), and the Pengcheng Talent-Medical Young Reserve Talent Training Project (No. XWRCHT20220012). Y.A.C. reports funding from the Postgraduate Research & Practice Innovation Program of Jiangsu Province (No. SJCX23_1399). Z.D.S. reports funding from the Jiangsu Provincial TCM Science and Technology Development Plan Project (No. MS2023081). 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. The 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. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Lenis AT, Lec PM, Chamie K, et al. Bladder Cancer: A Review. JAMA 2020;324:1980-91. [Crossref] [PubMed]
3. Sun Z, Jing C, Xiao C, et al. An autophagy-related long non-coding RNA prognostic signature accurately predicts survival outcomes in bladder urothelial carcinoma patients. Aging (Albany NY) 2020;12:15624-37. [Crossref] [PubMed]
4. Witjes JA, Bruins HM, Cathomas R, et al. European Association of Urology Guidelines on Muscle-invasive and Metastatic Bladder Cancer: Summary of the 2020 Guidelines. Eur Urol 2021;79:82-104. [Crossref] [PubMed]
5. Guo M, Lin J, Cao X, et al. Genetic variants in hypoxia-inducible factor pathway are associated with colorectal cancer risk and immune infiltration. J Cell Mol Med 2024;28:e18019. [Crossref] [PubMed]
6. Li Y, Zhao L, Li XF. Hypoxia and the Tumor Microenvironment. Technol Cancer Res Treat 2021;20:15330338211036304. [Crossref] [PubMed]
7. Zhou Q, Huang W, Xiong J, et al. CDCA8 promotes bladder cancer survival by stabilizing HIF1α expression under hypoxia. Cell Death Dis 2023;14:658. [Crossref] [PubMed]
8. Gilkes DM, Semenza GL, Wirtz D. Hypoxia and the extracellular matrix: drivers of tumour metastasis. Nat Rev Cancer 2014;14:430-9. [Crossref] [PubMed]
9. Colaprico A, Silva TC, Olsen C, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016;44:e71. [Crossref] [PubMed]
10. Goldman MJ, Craft B, Hastie M, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol 2020;38:675-8. [Crossref] [PubMed]
11. Lee JS, Leem SH, Lee SY, et al. Expression signature of E2F1 and its associated genes predict superficial to invasive progression of bladder tumors. J Clin Oncol 2010;28:2660-7. [Crossref] [PubMed]
12. Choi W, Porten S, Kim S, et al. Identification of distinct basal and luminal subtypes of muscle-invasive bladder cancer with different sensitivities to frontline chemotherapy. Cancer Cell 2014;25:152-65. [Crossref] [PubMed]
13. Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics 2007;23:1846-7. [Crossref] [PubMed]
14. Lai H, Cheng X, Liu Q, et al. Single-cell RNA sequencing reveals the epithelial cell heterogeneity and invasive subpopulation in human bladder cancer. Int J Cancer 2021;149:2099-115. [Crossref] [PubMed]
15. Stelzer G, Rosen N, Plaschkes I, et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr Protoc Bioinformatics 2016;54:1.30.1-1.30.33.
16. Gao J, Aksoy BA, Dogrusoz U, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal 2013;6:pl1. [Crossref] [PubMed]
17. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
18. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
19. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
20. Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28:1747-56. [Crossref] [PubMed]
21. Yang W, Soares J, Greninger P, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res 2013;41:D955-61. [Crossref] [PubMed]
22. Barretina J, Caponigro G, Stransky N, et al. The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature 2012;483:603-7. [Crossref] [PubMed]
23. Guzmán L, Villalón K, Marchant MJ, et al. In vitro evaluation and molecular docking of QS-21 and quillaic acid from Quillaja saponaria Molina as gastric cancer agents. Sci Rep 2020;10:10534. [Crossref] [PubMed]
24. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem 2010;31:455-61. [Crossref] [PubMed]
25. Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell 2021;184:3573-3587.e29. [Crossref] [PubMed]
26. Aran D, Looney AP, Liu L, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 2019;20:163-72. [Crossref] [PubMed]
27. Trapnell C, Cacchiarelli D, Grimsby J, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol 2014;32:381-6. [Crossref] [PubMed]
28. Haghverdi L, Büttner M, Wolf FA, et al. Diffusion pseudotime robustly reconstructs lineage branching. Nat Methods 2016;13:845-8. [Crossref] [PubMed]
29. Wicks EE, Semenza GL. Hypoxia-inducible factors: cancer progression and clinical translation. J Clin Invest 2022;132:e159839. [Crossref] [PubMed]
30. Keane MP, Belperio JA, Xue YY, et al. Depletion of CXCR2 inhibits tumor growth and angiogenesis in a murine model of lung cancer. J Immunol 2004;172:2853-60. [Crossref] [PubMed]
31. Akiyama K, Ohga N, Hida Y, et al. Tumor endothelial cells acquire drug resistance by MDR1 up-regulation via VEGF signaling in tumor microenvironment. Am J Pathol 2012;180:1283-93. [Crossref] [PubMed]
32. Kikuchi H, Maishi N, Annan DA, et al. Chemotherapy-Induced IL8 Upregulates MDR1/ABCB1 in Tumor Blood Vessels and Results in Unfavorable Outcome. Cancer Res 2020;80:2996-3008. [Crossref] [PubMed]
33. Reiswich V, Könemann S, Lennartz M, et al. Large-scale human tissue analysis identifies Uroplakin 1a as a putative diagnostic marker for urothelial cancer. Pathol Res Pract 2022;237:154028. [Crossref] [PubMed]
34. Dong Y, Zhu GY, Hao L, et al. High S100A7 expression is associated with early muscle invasion and poor survival in bladder carcinoma. Ann Diagn Pathol 2022;56:151847. [Crossref] [PubMed]
35. Bresnick AR, Weber DJ, Zimmer DB. S100 proteins in cancer. Nat Rev Cancer 2015;15:96-109. [Crossref] [PubMed]
36. Barsoum IB, Koti M, Siemens DR, et al. Mechanisms of hypoxia-mediated immune escape in cancer. Cancer Res 2014;74:7185-90. [Crossref] [PubMed]
37. Abbaszadegan MR, Riahi A, Forghanifard MM, et al. WNT and NOTCH signaling pathways as activators for epidermal growth factor receptor in esophageal squamous cell carcinoma. Cell Mol Biol Lett 2018;23:42. [Crossref] [PubMed]
38. Wang J, Ye C, Chen C, et al. Glucose transporter GLUT1 expression and clinical outcome in solid tumors: a systematic review and meta-analysis. Oncotarget 2017;8:16875-86. [Crossref] [PubMed]
39. Vasilevskaya I, O’Dwyer PJ. Role of Jun and Jun kinase in resistance of cancer cells to therapy. Drug Resist Updat 2003;6:147-56. [Crossref] [PubMed]
40. Venkatraman S, Balasubramanian B, Thuwajit C, et al. Targeting MYC at the intersection between cancer metabolism and oncoimmunology. Front Immunol 2024;15:1324045. [Crossref] [PubMed]
41. Alharbi KS, Javed Shaikh MA, Afzal O, et al. An overview of epithelial growth factor receptor (EGFR) inhibitors in cancer therapy. Chem Biol Interact 2022;366:110108. [Crossref] [PubMed]
42. Gao H, Nepovimova E, Heger Z, et al. Role of hypoxia in cellular senescence. Pharmacol Res 2023;194:106841. [Crossref] [PubMed]
43. George AL, Rajoria S, Suriano R, et al. Hypoxia and estrogen are functionally equivalent in breast cancer-endothelial cell interdependence. Mol Cancer 2012;11:80. [Crossref] [PubMed]
44. Huang J, Yao X, Zhang J, et al. Hypoxia-induced downregulation of miR-30c promotes epithelial-mesenchymal transition in human renal cell carcinoma. Cancer Sci 2013;104:1609-17. [Crossref] [PubMed]
45. Wu X, Lv D, Cai C, et al. A TP53-Associated Immune Prognostic Signature for the Prediction of Overall Survival and Therapeutic Responses in Muscle-Invasive Bladder Cancer. Front Immunol 2020;11:590618. [Crossref] [PubMed]
46. Li H, Lu H, Cui W, et al. A TP53-based immune prognostic model for muscle-invasive bladder cancer. Aging (Albany NY) 2020;13:1929-46. [Crossref] [PubMed]
47. Winerdal ME, Krantz D, Hartana CA, et al. Urinary Bladder Cancer Tregs Suppress MMP2 and Potentially Regulate Invasiveness. Cancer Immunol Res 2018;6:528-38. [Crossref] [PubMed]
48. Wesolowski R, Markowitz J, Carson WE 3rd. Myeloid derived suppressor cells - a new therapeutic target in the treatment of cancer. J Immunother Cancer 2013;1:10. [Crossref] [PubMed]
49. Chen S, Zhang Q, Lu L, et al. Heterogeneity of neutrophils in cancer: one size does not fit all. Cancer Biol Med 2022;19:1629-48. [Crossref] [PubMed]
50. Ocana A, Nieto-Jiménez C, Pandiella A, et al. Neutrophils in cancer: prognostic role and therapeutic strategies. Mol Cancer 2017;16:137. [Crossref] [PubMed]
51. Xiao LY, Su YL, Huang SY, et al. Chitinase 3-like-1 Expression in the Microenvironment Is Associated with Neutrophil Infiltration in Bladder Cancer. Int J Mol Sci 2023;24:15990. [Crossref] [PubMed]
52. Zhang J, Liu H, Zhang W, et al. Identification of lncRNA-mRNA Regulatory Module to Explore the Pathogenesis and Prognosis of Melanoma. Front Cell Dev Biol 2020;8:615671. [Crossref] [PubMed]
53. Tang YY, Sheng SY, Lu CG, et al. Effects of Glycogen Synthase Kinase-3β Inhibitor TWS119 on Proliferation and Cytokine Production of TILs From Human Lung Cancer. J Immunother 2018;41:319-28. [Crossref] [PubMed]
54. Chen YQ, Zheng L, Aldarouish M, et al. Wnt pathway activator TWS119 enhances the proliferation and cytolytic activity of human γδT cells against colon cancer. Exp Cell Res 2018;362:63-71. [Crossref] [PubMed]
55. Dong Y, Hao L, Shi ZD, et al. Solasonine Induces Apoptosis and Inhibits Proliferation of Bladder Cancer Cells by Suppressing NRP1 Expression. J Oncol 2022;2022:7261486. [Crossref] [PubMed]
56. Ma L, Wei HL, Wang KJ, et al. Rhein promotes TRAIL-induced apoptosis in bladder cancer cells by up-regulating DR5 expression. Aging (Albany NY) 2022;14:6642-55. [Crossref] [PubMed]
57. Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity. Nat Rev Immunol 2018;18:35-45. [Crossref] [PubMed]
58. Woo HJ, Kim SH, Kang HJ, et al. Continuous centrifugal microfluidics (CCM) isolates heterogeneous circulating tumor cells via full automation. Theranostics 2022;12:3676-89. [Crossref] [PubMed]