
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)13669-9
10.1016/j.heliyon.2024.e37638
e37638
Research Article
Disulfidptosis and ferroptosis related genes define the immune microenvironment and NUBPL serves as a potential biomarker for predicting prognosis and immunotherapy response in bladder cancer
Zhang Xuezhou zhangxuezhou1128@163.com
a1
Hong Baoan a1
Li Hongwei b1
Sun Zhipeng a
Zhao Jiahui a
Li Mingchuan a
Wei Dechao a
Wang Yongxing a
Zhang Ning Azurology_zn@126.com
a⁎
a Department of Urology, Beijing Anzhen Hospital, Capital Medical University, Beijing, PR China
b Department of Pathology, Key Laboratory of Carcinogenesis and Translational Research (Ministry of Education/Beijing), Peking University Cancer Hospital & Institute, Beijing, PR China
⁎ Corresponding author. Department of Urology, Beijing Anzhen Hospital, Capital Medical University, 2 Anzhen Road, Chaoyang District, Beijing, 100029, PR China. Azurology_zn@126.com
1 The authors contributed equally to this work.

07 9 2024
15 9 2024
07 9 2024
10 17 e3763828 5 2024
5 8 2024
6 9 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Background

Ferroptosis and disulfidptosis are regulatory forms of cell death that play an important role in tumorigenesis and progression. However, few biomarkers about disulfidptosis and ferroptosis related genes (DFRGs) have been developed to predict the prognosis of bladder cancer (BC).

Methods

We conducted a bioinformatics analysis using public BC datasets to examine the prognostic significance of differentially expressed DFRGs. A Lasso regression was employed to create a prognostic prediction model from these DFRGs. Hub DFRGs that play a role in immunotherapy response and immunoregulation were pinpointed. Immunohistochemistry (IHC) experiment was performed to assess NUBPL and c-MYC expression in BC patients who underwent surgery or received immune checkpoint inhibitor (ICI) immunotherapy at Peking University Cancer Hospital.

Results

We constructed a valid model to predict the prognosis of BC based on DFRGs and performed relevant validation, the results demonstrated that the model was an independent prognostic factor for BC. Further analysis indicated that the model score, combined with the expression of various immune factors and tumor mutation burden (TMB), could predict the prognosis for BC. In addition, we also found that NUBPL was strongly associated with prognosis and response to ICI treatment, and NUBPL may influence BC malignant progression through the c-MYC pathway.

Conclusions

Our research findings highlight the satisfactory predictive value of DFRGs in the immune microenvironment and suggest that NUBPL may be a highly promising biomarker for predicting the prognosis and efficacy of ICI treatment in BC patients.

Graphical abstract

Image 1

Keywords

Bladder cancer
Disulfidptosis
Ferroptosis
NUBPL
Immunotherapy
==== Body
pmc1 Introduction

According to the latest data from the American Cancer Society (ACS) in 2023, the incidence of BC is the 4th highest among male patients, with approximately 62,420 new cases along with 12,160 new deaths [1]. The high incidence and mortality of BC directly threaten human health and survival. In the last 30 years, immunotherapy, represented by immune checkpoint inhibitor (ICI), has prolonged the survival and brought new hope to patients with advanced BC. However, immunotherapy is still facing many difficulties, such as low overall response rate, lack of more accurate predictive indicators, and inability to accurately screen the potential beneficiary population, which limits the wide application of immunotherapy [2]. Therefore, further exploration of effective biomarkers for predicting survival prognosis and response to immunotherapy for clinical application is a very urgent task in the field of BC treatment.

Regulated cell death (RCD) plays an important role in the growth, invasion and apoptosis of tumor cells and also provides new strategies for tumor treatment [3,4]. Among the various forms of RCD, ferroptosis stands out as a unique type characterized by its dependence on iron and reactive oxygen species (ROS), which distinguishes it biochemically, genetically, and morphologically from other types of cell death, including necroptosis, pyroptosis, autophagy, and apoptosis [5,6]. In a noteworthy advancement, recent research has revealed that cells with high SLC7A11 expression undergoing glucose deprivation can accumulate intracellular disulfides, leading to a novel and distinct form of cell death. This process, which differs from both ferroptosis and apoptosis, has been termed disulfidptosis by the investigators [7]. However, disulfidptosis and ferroptosis related genes (DFRGs) and the therapeutic value of their interconnectivity had never been explored in BC patients.

Therefore, in this research we focused on discovering the DFRGs and its prognostic value in BC and established a signature based on DFRGs. This study stands out for its innovative approach in bladder cancer modeling by integrating the processes of ferroptosis and disulfidptosis to construct a novel model based on DFRGs, which exhibits notable predictive capabilities. In addition, we collected clinical data from BC patients in the Peking University Cancer Hospital (PUCH) and confirmed for the first time that NUBPL was strongly associated with the prognosis and the response to ICI immunotherapy in BC. Furthermore, NUBPL may influence BC malignant progression through the c-MYC pathway.

2 Materials and methods

2.1 Data collection and acquisition of DFRGs

We obtained RNA-seq expression profiles and clinical data of BC patients from two sources: The Cancer Genome Atlas (TCGA) BC cohort (TCGA-BLCA) and the Gene Expression Omnibus (GEO) BC cohort (GSE13507), accessible via their respective websites. The TCGA dataset provided comprehensive patient details, including gender, age, cancer stage, TNM classification, overall survival (OS), and current survival status. Our study included 600 BC samples (TCGA:412; GSE13507:188) and 87 non-tumorous bladder tissues from BC patients (TCGA:19; GSE13507:68).

We identified 484 genes associated with ferroptosis by exploring the FerrDb database (http://www.zhounan.org/ferrdb/) and collated a list of 10 genes linked to disulfidptosis from the existing scientific literature [7]. Additionally, we sourced immune cell composition data from the Tumor Immune Estimation Resource (TIMER) platform.

The IMvigor210 trial, a phase II clinical study, evaluated the efficacy of the programmed death-ligand 1 (PD-L1) inhibitor atezolizumab in patients with locally advanced or metastatic urothelial carcinoma who had not responded to prior platinum-based chemotherapy [8]. The inclusion criteria for this cohort of patients were: (1) patients with pathologically confirmed BC; (2) patients with complete immunotherapy response information. Exclusion criteria: (1) patients with pathology of non-BC; (2) imperfect information on immunotherapy efficacy and gene expression. Eventually our study included 195 BC samples in IMvigor210 cohort.

2.2 Identification of differentially expressed DFRGs

We employed the "limma" package in R to identify DFRGs that were differentially expressed between BC and normal tissue in the TCGA-BLCA cohort, using an FDR <0.05 and |log2FC| > 1 as thresholds. Heatmaps illustrated the most significantly altered genes, while volcano plots visualized the overall differential expression landscape.

2.3 Identification of prognosis-related DFRGs

To validate the association between the expression of DFRGs and BC patient outcomes, we analyzed two cohorts: the TCGA BC (TCGA-BLCA) and the GEO BC (GSE13507). Using the "survival" package in R, univariate Cox regression analysis determined the DFRGs with prognostic significance. A network plot visualized the connections among the expression levels of these key prognostic DFRGs. For TCGA BC, gene-level copy number datasets were sourced from the Genomic Data Commons (GDC) through UCSC XENA (https://xena.ucsc.edu/). We used the "RCircos" package in R to examine the frequency of copy number variations (CNV.frequency%) and the genomic alterations in prognosis-related DFRGs.

2.4 Consensus clustering analysis of DFRGs

BC patients were stratified into distinct groups through consensus clustering analysis performed with the "ConsensusClusterPlus" package in R. Kaplan-Meier analysis then compared survival outcomes across these different clusters. To verify the clustering's robustness, we applied principal component analysis (PCA), uniform manifold approximation and projection (UMAP), and t-distributed stochastic neighbor embedding (tSNE).

We analyzed expression patterns of prognosis-associated DFRGs and their clinical correlations, presenting the findings in heatmaps and boxplots. We also explored immune infiltration differences across BC subtypes. R packages "GSVA" and "GSEABase" facilitated Gene Set Variation Analysis (GSVA) and Gene Set Enrichment Analysis (GSEA), respectively, highlighting KEGG pathway differences among clusters [9]. We also compared the expression of immune checkpoint genes programmed death-ligand 1 (PD-L1), cytotoxic T-lymphocyte associated protein 4 (CTLA4) and programmed death 1 (PD-1) across clusters.

2.5 Construction of DFRGs-based prognostic prediction model

R software package "survival" was used to perform univariate Cox regression analysis to determine the prognosis-related DFRGs. Incorporating these genes into a Lasso-Cox model, we conducted cross-validation to establish a DFRGs-based prognostic index for BC prognosis using the "glmnet" package in R. The formula for the DFRG-based signature risk score was Ʃ(βi × Expi), where i represents the DFRGs count. Lasso regression, ideal for genomic feature-rich datasets, shrinks some coefficients to improve model generalizability.

The cohort was split into training and testing groups, further divided into high- and low-risk categories based on median risk scores. Kaplan-Meier analysis compared overall survival (OS) in both high- and low-risk groups within the training and testing sets. We quantified the predictive accuracy of our prognostic signature using ROC curve AUC values. Multivariate Cox regression determined the model's significance, and a heatmap correlated risk scores with signature genes.

2.6 Identification and evaluation of a predictive nomogram

We constructed a nomogram incorporating clinicopathological features and risk scores to predict outcomes in BC patients. Calibration curves were generated to compare the predicted versus actual survival rates, and decision curve analysis (DCA) was utilized to determine the clinical utility of the nomogram. DCA was widely used in clinical medicine as a tool for evaluating the effectiveness of medical decision-making. With the horizontal axis as the threshold and the vertical axis as the net benefit, the DCA curves can visualize the performance of different strategies in the balance of avoiding over-treatment and under-diagnosis [10].

2.7 Relationship between risk score and immune cell infiltration

To quantify the composition of immune cells within the tumor, CIBERSORT was conducted [11]. CIBERSORT estimated the relative proportions of different immune cell types, ensuring the sum of proportions in each sample equaled one. Spearman rank correlation was used to investigate the association between risk scores and the levels of immune infiltration. Additionally, we calculated the estimate score for the expression profiles of high- and low-risk groups in R.

2.8 Single cell data and drug sensitivity analysis

The Tumor Immune Single-Cell Hub (TISCH) database, a large-scale online database of single-cell RNA-seq, was accessed to examine the heterogeneity of the tumor microenvironment (TME) across various cell types and datasets [12]. Drug sensitivity, a critical factor for predicting treatment response, was predicted for BC samples in the high- and low-risk groups using the "oncoPredict" and "ggplot2″ R packages, referencing the Genomics of Drug Sensitivity in Cancer (GDSC) database [13].

2.9 Relationship between risk score and immunotherapy response

Kaplan-Meier analysis was used to compare the survival outcomes difference between high- and low-risk score groups in IMvigor210-BC cohort. Subsequently, we analyzed clinicopathological features, immune factors and their relationships with the risk score. We further analyzed whether the risk score integrated with various immune factors including PD-L1, CTLA4 and TMB expression could predict the prognosis in BC patients. We further investigated the relationship between DFRGs and immunotherapy responses in the IMvigor210-BC cohort.

2.10 RNA Isolation and reverse Transcription-Quantitative Polymerase chain Reaction (RT-qPCR)

To further validate the expression levels of DFRGs that were utilized to develop the signature in BC tissues and normal tissues, RT-qPCR was performed. The Total RNA Extraction Kit (Cat:R1200, Solarbio Life Sciences, Beijing, China) was used to extract total RNA from tissue specimens based on the manufacturer's instructions. For the detection of mRNA levels, the total RNA (500 ng) was transcribed into cDNA using a PrimeScript™ RT reagent kit (Perfect Real Time). Primer sequences were as follows: NUBPL, forward primer 5′-GCCTGGCCCGTTGTTATGAA-3′ and reverse primer 5′-GCCTCATTAGGTTGCTCTGTGA-3’; GAPDH, forward primer 5′- GAAAGCCTGCCGGTGACTAA-3′ and reverse primer 5′- GCCCAATACGACCAAATCAGAG-3’. The relative expression levels of NUBPL were determined by the comparative threshold (Ct) cycle method using Formula 2−(ΔΔCt).

2.11 Immunohistochemistry (IHC)

NUBPL is a disulfidptosis related gene and previous study has shown that NUBPL plays an important role in colorectal cancer invasion and migration [14], while the role of NUBPL in BC has not been reported. In addition, c-MYC is the most common human oncogene. Therefore, we would like to preliminarily explore the expression of NUBPL and its relationship with c-MYC by IHC in BC. Tissue samples were collected from 78 BC patients who underwent surgery from October 2014 to December 2021 in the Peking University Cancer Hospital (PUCH). In the PUCH-BC cohort, no patients received radiotherapy, immunotherapy or adjuvant chemotherapy before surgery and follow-up was performed every 3 months. Informed consent was acquired from all patients and the study was approved by the ethics committee (Institutional Review Board approval number: 2020KT143-GZ01) and conducted in accordance with Good Clinical Practice and the Declaration of Helsinki and its latest amendments. All clinical characteristics were compiled in Table 1. We further validated the association between NUBPL expression and BC prognosis using the TCGA-BLCA cohort, GSE13507 cohort, and PUCH-BC cohort.Table 1 Correlation between clinicopathological features and NUBPL expression in PUCH-BC cohort.

Table 1Clinical features		Total:78(%)	NUBPL expression	P-valuea	
negative	positive	
Age	<60	19(24.36 %)	5	14		
≥60	59(75.64 %)	13	46	0.7575	
Gender	Male	50(64.10 %)	13	37		
Female	28(35.90 %)	5	23	0.5768	
T	T1-T2	49(62.82 %)	10	39		
T3-T4	29(37.18 %)	8	21	0.5799(T1-2 vs. T3-4)	
N	N0-N1	59(75.64 %)	12	47		
N2-N3	19(24.36 %)	6	13	0.3542(N0-1 vs. N2-3)	
M	M0	69(88.46 %)	17	52		
M1	9(11.54 %)	1	8	0.6759	
c-MYC	Low	19(24.36 %)	8	11		
High	59(75.64 %)	10	49	0.0319	
a Fisher's exact test. PUCH: Peking University Cancer Hospital, BC: bladder cancer.

In a subset of 30 patients undergoing immunotherapy at PUCH from November 2019 to August 2021, treatments were administered every three weeks, either as monotherapy or in combination with platinum-based chemotherapy. Treatment efficacy was evaluated using the immune-modified Response Evaluation Criteria in Solid Tumors (iRECIST) version 1.1, with outcomes including complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). All clinical information and efficacy evaluation was summarized in Table 2.Table 2 Clinical information, efficacy evaluation and NUBPL expression of the 30 BC patients treated with ICI in PUCH cohort.

Table 2Clinical features	Total:30(%)	NUBPL expression	P-valuea	
low	high	
Age	<60	8(26.67 %)	3	5		
≥60	22(73.33 %)	6	16	0.6662	
Gender	Male	17(56.67 %)	3	14		
Female	13(43.33 %)	6	7	0.1232	
T	T1	8(26.67 %)	4	4		
T2	5(16.67 %)	1	4		
T3	12(40.00 %)	2	10		
T4	3(10.00 %)	1	2		
Tx	2(6.67 %)	1	1	0.4097(T1-2 vs. T3-4)	
N	N0	7(23.33 %)	1	6		
N1	10(33.33 %)	5	5		
N2	10(33.33 %)	2	8		
N3	3(10.00 %)	1	2	0.6908(N0-1 vs. N2-3)	
M	M0	14(46.67 %)	5	9		
M1	15(50.00 %)	4	11		
Mx	1(3.33 %)	0	1	0.6999(M0 vs. M1)	
Efficacy	CR	2(6.67 %)	2	0		
PR	4(13.33 %)	2	2		
SD	14(46.67 %)	3	11		
PD	10(33.33 %)	2	8	0.0492(CR + PR vs. SD + PD)	
a Fisher's exact test. BC: bladder cancer, ICI: immune checkpoint inhibitor, PUCH: Peking University Cancer Hospital.

For histological analysis, consecutive tissue sections were prepared and subjected to immunohistochemical staining. After deparaffinization and antigen retrieval, sections were incubated with NUBPL polyclonal (No:17393-1-AP, Proteintech Company) and c-MYC monoclonal antibodies (No:67447-1-Ig, Proteintech Company). Staining was executed using the avidin-biotin peroxidase method and counterstained with hematoxylin. Positive control is human colorectal cancer tissue. Immunohistochemical staining of the samples was scored and analyzed by a senior pathologist. Five randomly selected fields of view were categorized according to the percentage of staining (0: 0 %, 1: 1–25 %, 2: 26–50 %, 3: 51–75 %, 4: 76–100 %) and the degree of staining (0: negative staining, 1: yellowish, 2: light brown, 3: dark brown). The final score was the multiplication of the two scores and score ≥8 was considered positive.

2.12 Statistical analysis

Statistical analyses were performed using R software (version 4.0.2), with p < 0.05 denoting significance. Student's t-test was used to compare two experimental groups and one-way analysis of variance (ANOVA) was used to compare several treatments groups. To ensure the integrity and dependability of our statistical findings, we applied multiple test corrections for differential expression analysis and batch effects were considered and removed when combining TCGA and GEO data.

3 Results

3.1 Identification of differentially expressed DFRGs related to prognosis in BC

Fig. 1 outlines our study's workflow, where we identified 492 DFRGs, excluding shared genes (Supplementary Table 1). To validate the association between the expression of DFRGs and BC patient outcomes, we then pinpointed 75 differentially expressed DFRGs (Fig. 2A) using the "limma" package with our predefined criteria (Supplementary Table 2). The differential gene volcano plot is depicted in Fig. 2B.Fig. 1 Flowchart of the study: Fig. 1 shows the flowchart of our study.

Fig. 1

Fig. 2 Identification of prognosis-related DFRGs: (A) Seventy-five differentially expressed DFRGs between BC and normal samples (FDR< 0.05 and |log2FC| > 1); (B) The volcano map of all the differential DFRGs; (C) The network plot showed the relationship between the expression levels of the 25 prognosis-related DFRGs of rank; (D,E) The copy number variation frequency (CNV.frequency%) and the chromosome region and alteration of the above 25 prognosis-related DFRGs were screened.

Fig. 2

We conducted survival analysis on the 75 genes, and univariate Cox regression revealed 25 DFRGs with significant prognostic value in BC (p < 0.05) (Supplementary Fig. 1). A network plot illustrated the interactions among these 25 DFRGs (Fig. 2C). Fig. 2D provides a detailed display of the frequency of copy number variations (CNV.frequency%) for the aforementioned genes. Genes such as TFAP2A and SLC3A2 exhibit a notably higher frequency of "gain" compared to "loss". Conversely, genes like NT5DC2 and CDC25A show a significantly higher frequency of "loss" compared to "gain". Fig. 2E further illustrates the specific chromosomal locations of these variations for the genes mentioned. In this section, we identified differentially expressed DFRGs associated with BC prognosis.

3.2 BC molecular subgroups by using Consistent clustering

BC patients were stratified into distinct groups through consensus clustering analysis and Kaplan-Meier analysis then compared survival outcomes across these different clusters. We applied consensus clustering to the 25 prognostic DFRGs using the "Consensus Cluster Plus" package in R to better understand their impact on BC. The optimal number of clusters was determined to be two (k = 2), as depicted in Fig. 3A, which resulted in distinct prognostic outcomes as demonstrated by the OS analysis (p = 0.010), with cluster B exhibiting poorer survival (Fig. 3B). Multidimensional analyses including PCA, tSNE, and UMAP confirmed the distinct separation of the two subtypes (Fig. 3C–E). The expression of DFRGs across the subtypes was visualized using boxplots (Fig. 3F), and the immune cell infiltration patterns were outlined for each cluster (Fig. 3G). A heatmap detailed the expression of the 25 DFRGs and associated clinicopathological features for the two subtypes (Fig. 4A).Fig. 3 BC molecular subgroups by using Consistent clustering: (A) The above 25 prognosis-related DFRGs were used for Consensus Clustering, when k = 2, the cohort could be well classified into two subtypes; (B) Overall survival analysis showed a significant difference in prognosis between the two subtypes (p = 0.010) and DFRGcluster B has a worse overall survival; (C–E) PCA, tSNE and UMAP analyses were used to test the accuracy of this clustering and the results showed that the two clustering subtypes could be well identified; (F) The boxplot of the expression patterns with DFRGs in the two subtypes (∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001); (G) Immune infiltration patterns in the two subtype clusters (∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001).

Fig. 3

Fig. 4 Construction of DFRGs-based prognostic model: (A) Heat map of the above 25 DFRGs expression and corresponding clinicopathological features of two subtypes; (B) The GSVA analysis focus on the differential enrichment of KEGG pathways showed Cluster B, the poor prognostic group, is mainly involved in “pathways in cancer”, “ECM receptor interaction” pathways; (C) The GSEA enrichment analysis and the top 5 most significant pathways with Cluster B were displayed; (D–F) Immune-checkpoint gene PD-L1, CTLA4 and PD-1 were all highly expressed in the clusterB (∗∗∗p < 0.001); K-M survival curve showed that the DFRGs-based model was closely related to the OS of BC patients in (G) all- (p < 0.001), (H) test- (p < 0.001) and (I) train- (p < 0.001) groups; (J) The heat map indicated the expression patterns of the 7 hub DFRGs; (K) The risk score was significantly different among the two previous subtypes (p = 1.6e-09); (L) Alluvial diagram showed the changes of DFRGs-cluster, DFRGs-based model risk and living status.

Fig. 4

The distribution of the DFRGs in the clusters was assessed using the "GSVA" package in R, which identified differential KEGG pathway enrichment between clusters A and B (Fig. 4B). The top five pathways significantly associated with Cluster B were highlighted through GSEA (Fig. 4C), with full details provided in Supplementary Table 3. Expression levels of immune-checkpoint genes PD-L1, CTLA4, and PD-1 were significantly higher in cluster B (p < 0.001), which suggested that cluster B had a high immunoregulatory (Fig. 4D–F). In this section, we found that Cluster B, which is associated with poor prognosis, was characterized by enrichment in cancer-related pathways and extracellular matrix (ECM) receptor interactions.

3.3 Construction and validation of DFRGs-based prognostic model

We used COX regression analysis as well as Lasso-Cox methods in order to construct a prognostic model. Patients were divided into training and testing groups for prognostic model development and validation. Lasso regression (Supplementary Figs. 2A and B) yielded a model comprising 7 DFRGs (CAV1, MEG3, MIOX, NT5DC2, PVT1, SLC3A2, and SREBF1), detailed in Supplementary Table 4. Kaplan-Meier curves illustrated the model's strong association with OS in all patient groups (p < 0.001) (Fig. 4G–I), with high-risk scores correlating with reduced 5-year survival. ROC analysis confirmed the model's predictive accuracy (Supplementary Figs. 2C–E), and multivariate Cox regression verified its independence as a prognostic factor for BC (p < 0.001) (Supplementary Fig. 2F).

The heatmap in Fig. 4J displayed the expression of the 7 hub DFRGs, while Fig. 4K showed the risk score differences between the two subtypes. Alluvial diagrams are commonly utilized to categorize populations based on characteristics such as gender, age, and family history, illustrating the patterns of disease progression across different groups. Fig. 4L provides a detailed representation of the model risk scores and prognostic status assessments among patients with different clusters. Overall, cluster B is associated with poorer clinical outcomes compared to cluster A. In this section, we constructed and validated a prognostic model based on DFRGs.

3.4 Establishment of a prognostic nomogram for BC patients

A nomogram integrating the DFRG-based model and clinicopathological data was constructed (Fig. 5A), with calibration plots confirming its predictive reliability (Fig. 5B). The cumulative distribution curve indicated higher hazards for the nomoRisk high group (Fig. 5C). DCA demonstrated the nomogram's utility in predicting 1-, 3-, and 5-year survival in BC patients(Fig. 5D–F).Fig. 5 Establishment of a prognostic nomogram: (A) The nomogram was constructed combined the DFRGs-based model with clinicopathological information; (B) The calibration plot showed the nomogram had a good predictive performance; (C) Cumulative distribution curve showed nomoRisk high group had a higher cumulative hazard; (D–F) Decision curve analysis (DCA) showed the above nomogram exhibited as a good method for predicting 1-, 3-, and 5- years survival in BC patients.

Fig. 5

3.5 Correlation analysis of risk score, DFRGs and tumor immune activity

The immune landscape's role in tumorigenesis was explored through the analysis of immune cell content in the low- and high-risk groups (Fig. 6A). Differences in immune cell populations were significant, notably in plasma cells (p = 0.025), regulatory T cells (p = 0.026), resting NK cells (p = 0.021) and M0 Macrophages (p = 0.022) (Fig. 6B). Correlations between risk scores, DFRGs, and immune cells were illustrated in Fig. 6C, with regulatory T cells (p = 0.014), plasma cells (p = 0.015), naïve B cells (p = 0.016) negatively correlated, and resting NK cells (p = 0.0066), activated mast cells (p = 0.046), M0 Macrophages (p = 0.017) positively correlated with risk scores (Fig. 6D–I).Fig. 6 Correlation analysis of risk score, DFRGs and tumor immune activity: (A) The relative immune cell content of each sample from low- and high-groups were showed; (B) The content of immune cells between the high- and low-risk group and the results showed significant differences in plasma cells (p = 0.025), regulatory T cells (p = 0.026), resting NK cells (p = 0.021) and M0 Macrophages (p = 0.022); (C) Correlation between riskScore, DFRGs and immune cells were showed (∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001); (D–I) There was a negative correlation between regulatory T cells (p = 0.014), plasma cells (p = 0.015), naïve B cells (p = 0.016) and riskScore, a positive correlation between resting NK cells (p = 0.0066), activated mast cells (p = 0.046), M0 Macrophages (p = 0.017) and riskScore.

Fig. 6

3.6 Single-cell data analysis and drug sensitivity of DFRGs

Using single-cell (sc) RNA-seq data from the TISCH database (BLCA_GSE145281), we analyzed DFRG expression within the tumor microenvironment (TME). This data covered 5 medium cell types and 17 cell clusters (Supplementary Fig. 3A), with expression levels of 6 DFRGs evaluated (Supplementary Figs. 3B and C). The expression data for the MIOX gene are not included in the Results section, as they were unavailable in the single-cell dataset utilized for this study. We also investigated drug sensitivity in the two groups (Supplementary Fig. 4).

3.7 DFRGs score combined with immune factors could predict prognosis in BC

To explore the relationship between risk score and immunotherapy response, we performed survival analysis and the results showed that the high-risk group had a significant worse prognosis in BC patients received immunotherapy (p = 0.008) (Fig. 7A). ROC curve analysis for 1-, 3-, and 5-year overall survival (OS) demonstrated good predictive performance with the AUC values presented in Fig. 7B. Moreover, this score was associated with age (p = 0.0071), the grade (p < 0.001), the stage (p < 0.001), T stage (p = 0.0082), and metastasis status (p = 0.036) of BC, but not in gender (p = 0.15), N stage (p = 0.052), and immunescore(p = 0.25) (Supplementary Figs. 5A–H). Further analyses suggested that the DFRGs-based model risk score integrated with various immune factors including PD-L1(p < 0.001), CTLA4(p < 0.001) and TMB(p < 0.001) expression, could predict prognosis for BC patients (Fig. 7C–E).Fig. 7 DFRGs score combined with immune factors could predict prognosis in BC: (A) Survival analysis showed that the high-risk group had a significant worse prognosis in BC patients received immunotherapy (p = 0.008); (B) ROC curve analysis for OS at 1-, 3-, and 5- years showed the AUC values which indicated a good predictive performance (AUC at 1 years: 0.564; AUC at 3 years: 0.575; AUC at 5 years: 0.550); (C–E) Further analyses suggested that the DFRGs score integrated with various immune factors including PD-L1(p < 0.001), CTLA4(p < 0.001) and TMB(p < 0.001) expression could predict prognosis for BC patients; (F–H) The immune-checkpoint gene PD-L1, PD-1 and CTLA4 were all highly expressed in the high-risk group (p < 0.001); (I–L) High expression of EGR1(p = 0.026) and WWTR1(p = 0.049) had a significant worse response, while high expression of CDCA3(p = 0.015) and RRM2(p = 0.031) had a significant better response for BC patients who received immunotherapy(Non-response:PD/SD, Response:CR/PR, CR: complete response, PR: partial response, PD: progressive disease, SD: stable disease).

Fig. 7

The correlation between the DFRGs-based model's risk score and immune factors in BC was examined, revealing higher expression of immune-checkpoint genes in the high-risk group (p < 0.001), suggesting the high-risk group maybe had a high immunoregulatory (Fig. 7F–H). Prognostic analysis indicated significant survival advantages associated with high expression of PD-L1(p = 0.042), CTLA4(p = 0.003) and TMB(p < 0.001) (Supplementary Figs. 5I–K). Among the 25 DFRGs, certain genes like EGR1(p = 0.026) and WWTR1(p = 0.049) were associated with poorer responses, while CDCA3(p = 0.015) and RRM2(p = 0.031) correlated with better responses in patients receiving immunotherapy (Fig. 7I–L), potentially serving as predictive markers for ICI response in BC.

3.8 NUBPL was correlated with prognosis and immunotherapy response in BC

We would like to preliminarily explore the expression of NUBPL and its relationship with c-MYC by RT-qPCR and IHC in BC. NUBPL expression and its prognostic correlation were validated in TCGA-BLCA and GSE13507 cohorts, with higher expression levels in BC tissues (p = 2.507e-16) (Fig. 8A).Fig. 8 NUBPL was correlated with prognosis and immunotherapy response in BC: (A) NUBPL was highly expressed in BC tissue compared with normal bladder tissue in TCGA-BLCA and GSE13507 cohort (p = 2.507e-16); (B) NUBPL was upregulated in BC by RT-qPCR; (C) BC patients with high expressions of NUBPL had a significantly worse OS rate in TCGA-BLCA and GSE13507 cohort (p = 0.002); (D) BC patients with high expressions of NUBPL had a worse OS rate compared to patients with low expressions of NUBPL in PUCH-BC cohort (p = 0.027); (E) The protein expression levels of NUBPL and c-MYC in PUCH-BC and normal tissues (NUBPL, No:17393-1-AP, Proteintech Company; c-MYC, No:67447-1-Ig, Proteintech Company); (F) The c-MYC expression levels were generally increased in BC patients with high NUBPL expression.

Fig. 8

RT-qPCR were performed to validate the expression levels of NUBPL in BC tissues and normal tissues. The results demonstrated significant differences in the expression levels of NUBPL between BC and normal tissues. Compared to normal tissues, NUBPL was upregulated in BC (Fig. 8B). High expression levels of NUBPL were associated with poorer prognosis in TCGA-BLCA and GSE13507 cohorts (p = 0.002) (Fig. 8C). Similar findings were observed in the PUCH-BC cohort (p = 0.027) (Fig. 8D). IHC analysis in the PUCH-BC cohort showed high NUBPL expression in a majority of BC tissues (76.92 %, 60/78) (Fig. 8E). The predictive value of NUBPL for ICI response was assessed in 30 advanced or metastatic BC patients, with a significant difference observed between response groups (p = 0.0492), indicating that high NUBPL expression could predict poor ICI response (Table 2).

Then we preliminarily explored the expression of NUBPL and its relationship with c-MYC by IHC in BC (Fig. 8F). IHC results showed that c-MYC expression levels were generally increased in patients with high NUBPL expression. Fisher's exact test showed that there was a positive correlation between the expression levels of NUBPL and c-MYC (p = 0.0319) (Table 1), suggesting that NUBPL may influence BC malignant progression through the c-MYC pathway.

4 Discussion

Bladder cancer (BC) is the most common malignant tumor of the urinary system with the 2nd highest incidence after prostate cancer, and also is the major disease that directly threaten human health and survival [15,16]. Recent years, immunotherapy based on PD-L1/PD-1 immune checkpoints had brought light to the treatment of BC patients [17,18]. However, the overall response rate to PD-L1/PD-1 ICIs in BC is currently low and there is a paucity of related biomarkers [19], the above problems are the current difficulties in BC treatment. Identifying new, reliable biomarkers for predicting disease prognosis and immunotherapy response is crucial for tailoring individual treatments in BC.

In recent years, studies have shown that ferroptosis is a form of RCD that plays an important role in tumor biology and therapy [20]. Ferroptosis, an iron-dependent form of regulated cell death (RCD) marked by lipid peroxidation, is increasingly recognized for its role in cancer biology and treatment. This process is intricately linked to the cellular metabolism of lipids, iron, and amino acids [21]. Genes related to ferroptosis have been proposed as prognostic markers in various cancers, including those of the prostate, colon, and breast [[21], [22], [23], [24]]. Disulfidptosis, distinct from both apoptosis and ferroptosis, is characterized by an abnormal build-up of intracellular disulfides in cells with high SLC7A11 expression under glucose deprivation [7]. In our study, we innovatively integrated disulfidptosis with ferroptosis to create a prognostic model for BC, utilizing differentially expressed DFRGs. Our model demonstrated the significant prognostic potential of these biomarkers in BC.

Two DFRGs-clusters were identified, with DFRGs-cluster B showing poorer overall survival and associations with critical pathways for tumor cell invasion and metastasis, such as “pathways in cancer” and “ECM receptor interaction” [25,26]. Furthermore, we observed enhanced immune infiltration of CD4 and CD8 T cells, as well as increased expression of immune checkpoints PD-L1, CTLA4, and PD-1 in cluster B, suggesting cluster B heightened immune regulatory activity. These results suggested that DFRGs may influence tumorigenesis, disease progression, and the immune landscape in BC.

To account for patient heterogeneity, we developed a DFRGs-based model assessing disulfidptosis and ferroptosis patterns in BC, focusing on genes CAV1, MEG3, MIOX, NT5DC2, PVT1, SLC3A2, and SREBF1. The model's risk score was significantly associated with patient prognosis, expression of immune checkpoints, and the degree of immune cell infiltration. Prior research has established PD-L1 and TMB as key predictors for the efficacy of ICI therapy in BC [27,28]. Our findings suggest a strong link between DFRGs, ICI effectiveness, and the immune microenvironment in BC.

An intriguing discovery from our research is that high expression levels of EGR1 and WWTR1 correlate with poorer ICI responses, while increased expression of CDCA3 and RRM2 is associated with better outcomes. It's been shown that CD8+ T cells activated by immunotherapy can promote ferroptosis-driven lipid peroxidation in cancer cells, enhancing the anti-tumor effects of ICI therapy [29]. Furthermore, recent studies have connected ferroptosis with the immune microenvironment in cancers such as liver cancer, head and neck squamous cell carcinoma, and lung adenocarcinoma [[30], [31], [32]]. Thus, our research underscores the potential of DFRGs biomarkers in predicting responses to immunotherapy in BC patients.

NUBPL is involved in the assembly of the respiratory chain dehydrogenase complex in the inner membrane of assembled mitochondria, which is associated with mitochondrial metabolism but whose role in cancer has not been systematically studied [33]. Previous study has shown that NUBPL expression levels were increased in colorectal cancer and high expression of NUBPL could promote epithelial-mesenchymal transition (EMT) and metastatic phenotype in colorectal cancer cells [14]. Here in this research, we confirmed that NUBPL was highly expressed in BC and it was closely related to the prognosis and PD-L1 immunotherapy treatment of BC. Meanwhile, IHC results showed that c-MYC expression levels were generally increased in patients with high NUBPL expression, suggesting that NUBPL may promote the progression of BC through the c-MYC pathway. However, further studies are needed to elucidate the specific role and mechanism of NUBPL in BC progression.

Our preliminary study does have some limitations. Firstly, the initial sample size was relatively small, which constrains the statistical robustness and the broader applicability of our findings. To mitigate this issue, we are in the process of enlarging the sample size and aim to incorporate additional patient data in subsequent studies to enhance the study's reliability and representativeness. We plan to gather a more extensive cohort of BC patients through partnerships with other medical centers across the country. To guarantee the representativeness of the sample, we will employ random sampling techniques and ensure that the sample encompasses diverse ages, genders, and geographic locations. Furthermore, the diversity and scope of the sample will be augmented through multi-center collaboration and we will verify the correlation between NUBPL expression and prognosis and immunotherapy response in a larger patient population. Additionally, while bioinformatics analyses of public databases provide valuable insights, experimental studies are essential to confirm the role of NUBPL in cancer progression and its interaction with PD-L1 targeted treatments in BC. In future studies, we plan to use cellular and animal models to validate NUBPL on relevant signaling pathways and their role in tumor development and drug response.

5 Conclusion

Despite these constraints, our research underscores the promising role of DFRGs within the immune landscape of BC and proposes NUBPL as a prospective marker for prognostication and immunotherapy response prediction in this disease context.

Fundings

This study was supported by Beijing Hospitals Authority Youth Programme (QML20231105 ), Xisike-leader Oncology Research Foundation (Y-2019AZZD-0369 ) and Beijing 10.13039/100016966 Bethune Charitable Foundation (mnzl202003 ).

Ethical approval

Informed consent was acquired from all patients and the study was approved by the ethics committee (Institutional Review Board approval number, 2020KT143-GZ01) and conducted in accordance with Good Clinical Practice and the Declaration of Helsinki and its latest amendments.

Data availability statement

The TCGA-BLCA cohort dataset is available at The Cancer Genome Atlas (TCGA) website (https://xena.ucsc.edu/). The GSE13507 cohort dataset for this study is available at Gene Expression Omnibus (GEO) with accession number: GSE13507 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE13507).

CRediT authorship contribution statement

Xuezhou Zhang: Writing – original draft, Validation, Software, Data curation. Baoan Hong: Writing – review & editing, Resources, Funding acquisition, Data curation, Conceptualization. Hongwei Li: Validation, Resources, Data curation. Zhipeng Sun: Software. Jiahui Zhao: Validation, Resources. Mingchuan Li: Visualization, Resources. Dechao Wei: Resources. Yongxing Wang: Supervision, Resources. Ning Zhang: Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization.

Declaration of competing interest

The authors declare that they have no conflict of interest.

Appendix A Supplementary data

The following are the Supplementary data to this article.Multimedia component 1

Multimedia component 1

Multimedia component 2

Multimedia component 2

Multimedia component 3

Multimedia component 3

Multimedia component 4

Multimedia component 4

Fig. S1 Univariate Cox regression analysis: Univariate Cox regression analysis showed that 25 differentially expressed DFRGs had significant prognostic value in BC (p < 0.05).

Fig. S1

Fig. S2 Lasso regression analysis and validation of DFRGs-based prognostic model: (A,B) Lasso regression analysis were used to construct a prognostic model; (C–E) ROC curve analysis for OS at 1-, 3-, and 5- years showed the AUC values which indicated a good predictive performance of the DFRGs-based model in all-(AUC at 1 years: 0.697; AUC at 3 years: 0.649; AUC at 5 years: 0.648), test-(AUC at 1 years: 0.643; AUC at 3 years: 0.630; AUC at 5 years: 0.647) and train-(AUC at 1 years: 0.743; AUC at 3 years: 0.680; AUC at 5 years: 0.660) groups; (F) Multivariate Cox regression analysis confirmed that the DFRGs-based model was an independent prognostic factor for BC (p < 0.001).

Fig. S2

Fig. S3 Single-cell data analysis of DFRGs: (A) In the BLCA_GSE145281 dataset, there are 17 cell clusters and 5 medium cell types, and the distribution and number of various cell types are shown; (B,C) Expressions and percentages of the 6 DFRGs showed that SLC3A2 is uniformly expressed in the 5 cell types while NT5DC2 is mainly expressed in B cells.

Fig. S3

Fig. S4 Drug sensitivity of DFRGs: The potential sensitivity of high- and low-risk group to clinical drugs.

Fig. S4

Fig. S5 Correlations between model scores and clinical features, correlations between immune checkpoint genes and prognosis: (A-H) The score was associated with age (p = 0.0071), the grade (p < 0.001), the stage (p < 0.001), T stage (p = 0.0082), and metastasis status (p = 0.036) of BC, but not in gender (p = 0.15), N stage (p = 0.052), and immunescore(p = 0.25); (I–K) Prognostic analyses showed that OS have significant advantages in patients with high expression of PD-L1(p = 0.042), CTLA4(p = 0.003) and TMB(p < 0.001).

Fig. S5

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e37638.
==== Refs
References

1 Siegel R.L. Miller K.D. Wagle N.S. Cancer statistics, 2023 CA Cancer J Clin 73 1 2023 Jan 17 48 36633525
2 Patel V.G. Oh W.K. Galsky M.D. Treatment of muscle-invasive and advanced bladder cancer in 2020 CA Cancer J Clin 70 5 2020 Sep 404 423 32767764
3 Koren E. Fuchs Y. Modes of regulated cell death in cancer Cancer Discov. 11 2 2021 Feb 245 265 33462123
4 Peng F. Liao M. Qin R. Regulated cell death (RCD) in cancer: key pathways and targeted therapies Signal Transduct Target Ther 7 1 2022 Aug 13 286 35963853
5 Jiang X. Stockwell B.R. Conrad M. Ferroptosis: mechanisms, biology and role in disease Nat. Rev. Mol. Cell Biol. 22 4 2021 Apr 266 282 33495651
6 Lei G. Zhuang L. Gan B. Targeting ferroptosis as a vulnerability in cancer Nat. Rev. Cancer 22 7 2022 Jul 381 396 35338310
7 Liu X. Nie L. Zhang Y. Actin cytoskeleton vulnerability to disulfide stress mediates disulfidptosis Nat. Cell Biol. 25 3 2023 Mar 404 414 36747082
8 Rosenberg J.E. Hoffman-Censits J. Powles T. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial Lancet 387 10031 2016 May 7 1909 1920 26952546
9 Hnzelmann S. Castelo R. Guinney J. GSVA: gene set variation analysis for microarray and RNA-Seq data BMC Bioinf. 14 1 2013 7-7
10 Vickers A.J. Cronin A.M. Elkin E.B. Extensions to decision curve analysis, a novel method for evaluating diagnostic tests, prediction models and molecular markers BMC Med Inform Decis Mak 8 2008 Nov 26 53 19036144
11 Newman A.M. Liu C.L. Green M.R. Robust enumeration of cell subsets from tissue expression profiles Nat. Methods 12 5 2015 May 453 457 25822800
12 Sun D. Wang J. Han Y. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment Nucleic Acids Res. 49 D1 2021 Jan 8 D1420 D1430 33179754
13 Yang W. Soares J. Greninger P. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells Nucleic Acids Res. 41 Database issue 2013 Jan D955 D961 23180760
14 Wang Y. Wu N. Sun D. NUBPL, a novel metastasis-related gene, promotes colorectal carcinoma cell motility by inducing epithelial-mesenchymal transition Cancer Sci. 108 6 2017 Jun 1169 1176 28346728
15 van Hoogstraten L.M.C. Vrieling A. van der Heijden A.G. Global trends in the epidemiology of bladder cancer: challenges for public health and clinical practice Nat. Rev. Clin. Oncol. 2023 Mar 13
16 Tran L. Xiao J.F. Agarwal N. Advances in bladder cancer biology and therapy Nat. Rev. Cancer 21 2 2021 Feb 104 121 33268841
17 Balar A.V. Kamat A.M. Kulkarni G.S. Pembrolizumab monotherapy for the treatment of high-risk non-muscle-invasive bladder cancer unresponsive to BCG (KEYNOTE-057): an open-label, single-arm, multicentre, phase 2 study Lancet Oncol. 22 7 2021 Jul 919 930 34051177
18 Iacovelli R. Ciccarese C. Brunelli M. First-line avelumab for patients with PD-L1-positive metastatic or locally advanced urothelial cancer who are unfit for cisplatin Ann. Oncol. 33 11 2022 Nov 1179 1185 35926813
19 McGrail D.J. Pilié P.G. Rashid N.U. High tumor mutation burden fails to predict immune checkpoint blockade response across all cancer types Ann. Oncol. 32 5 2021 May 661 672 33736924
20 Zhang C. Liu X. Jin S. Ferroptosis in cancer therapy: a novel approach to reversing drug resistance Mol. Cancer 21 1 2022 Feb 12 47 35151318
21 Stockwell B.R. Friedmann Angeli J.P. Bayir H. Ferroptosis: a regulated cell death nexus linking metabolism, redox biology, and disease Cell 171 2 2017 Oct 5 273 285 28985560
22 Ghoochani A. Hsu E.C. Aslan M. Ferroptosis inducers are a novel therapeutic approach for advanced prostate cancer Cancer Res. 81 6 2021 Mar 15 1583 1594 33483372
23 Yan H. Talty R. Johnson C.H. Targeting ferroptosis to treat colorectal cancer Trends Cell Biol. 33 3 2023 Mar 185 188 36473802
24 Li P. Yang B. Xiu B. Development and validation of a robust ferroptosis-related gene panel for breast cancer disease-specific survival Front. Cell Dev. Biol. 9 2021 Nov 25 709180
25 Debnath J. Gammoh N. Ryan K.M. Autophagy and autophagy-related pathways in cancer Nat. Rev. Mol. Cell Biol. 2 2023 Mar 1 16
26 Huang J. Zhang L. Wan D. Extracellular matrix and its therapeutic potential for cancer treatment Signal Transduct Target Ther 6 1 2021 Apr 23 153 33888679
27 Inman B.A. Longo T.A. Ramalingam S. Atezolizumab: a PD-L1-blocking antibody for bladder cancer Clin. Cancer Res. 23 8 2017 Apr 15 1886 1890 27903674
28 Chan T.A. Yarchoan M. Jaffee E. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic Ann. Oncol. 30 1 2019 Jan 1 44 56 30395155
29 Wang W. Green M. Choi J.E. CD8+ T cells regulate tumour ferroptosis during cancer immunotherapy Nature 569 7755 2019 May 270 274 31043744
30 Huang Y. Wang S. Ke A. Ferroptosis and its interaction with tumor immune microenvironment in liver cancer Biochim. Biophys. Acta Rev. Canc 1878 1 2023 Jan 188848
31 Yu J. Chen Y. Pan X. Relationships of ferroptosis and pyroptosis-related genes with clinical prognosis and tumor immune microenvironment in head and neck squamous cell carcinoma Oxid. Med. Cell. Longev. 2022 2022 Oct 5 3713929
32 Deng B. Xiang J. Liang Z. Identification and validation of a ferroptosis-related gene to predict survival outcomes and the immune microenvironment in lung adenocarcinoma Cancer Cell Int. 22 1 2022 Sep 24 292 36153508
33 Sheftel A.D. Stehling O. Pierik A.J. Human ind1, an iron-sulfur cluster assembly factor for respiratory complex I Mol. Cell Biol. 29 22 2009 Nov 6059 6073 19752196
