Construction and verification of a prognostic model for bladder cancer based on disulfidptosis-related angiogenesis genes

Data gathering and processing

The study workflow is presented in Fig. 1A. An aggregate of 425 BLCA samples were procured from the TCGA repository (https://portal.gdc.cancer.gov/), comprising 406 tumor samples as well as 19 normal specimens. Among the 406 BLCA patients, 401 had complete survival data available for analysis (Table S1). Additionally, the GSE32548 cohort, comprising 131 BLCA patients, was obtained from the GEO repository (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32548) and used for external verification (Table S2). DRGs were sourced from a presently available publication (Xin et al., 2023). The 82 ARGs were retrieved from the GESA Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_ANGIOGENESIS) (Subramanian et al., 2005).

Unsupervised clustering analysis

Based on 11 DRGs, an unsupervised clustering analysis was executed on the 406 tumor specimens from the TCGA-BLCA cohort by applying the ConsensuClusterPlus software package (Wilkerson & Hayes, 2010). This analysis aimed to identify DRG-related subgroups. The key clinical information among the subgroups was organized and compared (Table S3). The “Survival” and “survminer” packages (version 4.22) were employed to generate Kaplan-Meier (KM) curves of OS, thereby determining the prognostic differences between the identified subgroups. P < 0.05 was indicative of statistical significance.

Differential analysis and enrichment analysis

To elucidate the mechanisms of DRGs, the “DESeq2” software package was applied to carry out differential expression analysis (Love, Huber & Anders, 2014). Genes exhibiting differential expression were identified in line with the criteria of P < 0.05, |log2FC|≥ 1. Following this step, the “clusterProfiler” package was leveraged to execute Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses (Yu et al., 2012), with identical filtering criteria. Furthermore, the “VennDiagram” software package was utilized to intersect disulfidptosis-related DEGs with the selected 82 ARGs.

Formulation and verification of a model based on disulfidptosis-related ARGs

Through LASSO regression analysis, two key genes were identified. The risk score was derived via the formula: Risk Score = Σ (Expi × Coefi), in which Coefi denotes the risk coefficient, while Expi signifies the expression level of each gene. BLCA samples were then segmented into low-risk and high-risk cohorts leveraging the median risk score. Prognostic analysis, along with ROC curve analysis, was carried out via “ggrisk” “survival” and “rms” software packages. In the test cohort, the robustness of the model was further affirmed.

Clinical relevance analysis

To ascertain the influence of BLCA-linked variables on the outcome of prognosis, both univariate and multivariate Cox regression analyses were implemented. Subsequently, we probed the linkage between the prognostic risk model and various clinical features, covering age, gender, clinical stage, as well as pathological TNM stage, as obtained from TCGA. To forecast BLCA patients’ survival probability, the “rms” and “survival” packages were applied to formulate a nomogram. For each patient, a total score was derived from the values assigned to the predictive factors, and this total score was thereafter used to project the likelihood of survival (Park, 2018). The efficiency of the constructed nomogram was evaluated by calibration curves.

Immunological analysis and tumor mutation burden

The “CIBERSORT” package was employed to quantify immune cell’s degree of infiltration in both high- and low-risk cohorts (Newman et al., 2015). Differences in immune checkpoint expression between the high- and low-risk cohorts were compared via “tinyarray” software package. Subsequently, the R package “estimate” was leveraged to estimate and compare Stromal scores, Immune scores, and ESTIMATE scores across the two cohorts at risk. Online tool TIDE (https://tide.nki.nl/) was then implemented to assess the potential for evasion of immune response. The TCGA database provided somatic mutation information for BLCA patients, while waterfall plots illustrating the gene mutation distribution between the high- and low-risk segments were generated via R package “maftools.” Survival curves integrating risk scores with tumor mutation burden (TMB) were constructed and analyzed to gauge the prognosis of BLCA patients.

Cell cultivation and qRT-PCR

The standard bladder epithelial cell strain (SV-HUC-1) and the bladder cancer cell strain (T24 and SW1710) were sourced from the Cell Bank of the Chinese Academy of Sciences (Shanghai). These cell strains were cultured in RPMI 1640 medium (Gibco, Waltham, MA, USA), enriched by 10% fetal bovine serum (Gibco, Waltham, MA, USA) and 1% penicillin/streptomycin (Gibco, Waltham, MA, USA), and incubated at 37 °C in 5% CO2. Total RNA was isolated using TRIzol reagent (Invitrogen, Waltham, MA, USA). Reverse transcription was performed using HiScript III SuperMix (Vazyme, Beijing, China). qPCR assays were carried out with a QuantStudio 3 instrument (Thermo, China). Briefly, Upon the denaturation at 95 °C for 1 min, these cell strains were amplified at 95 °C for 5 s and 60 °C for 15 s for 40 cycles. The mRNA expression levels of COL5A2 and SCG2 were quantified through the 2−ΔΔCt approach. The primers utilized in these reactions were presented here:

β-ACTIN forward primer: CGGGAAATCGTGCGTGAC; reverse primer: CAGGCAGCTCGTAGCTCTT

COL5A2 forward primer: GGTCTCAGTTCGCTTATGG; reverse primer: TGTAAGTGATGTTCTGGGAGG

SCG2 forward primer: ACCAGACCTCAGGTTGGAAAA; reverse primer: AAGTGGCTTTCATCGCCATTT

Statistical analysis

R software (version 4.4.1) and GraphPad Prism (version 9.2.0) were employed to execute all statistical analyses. The Spearman correlation coefficient was utilized to determine the correlation between two variables. Continuous variables were subject to the Student’s t-test or the Mann-Whitney U test, while the chi-square test was adopted to evaluate categorical variables. Differences were deemed to be statistically significant at P < 0.05.

Identification of disulfidptosis-related subgroups

In the present study, a comprehensive literature review was conducted to identify 24 genes linked to disulfidptosis. Among these 24 DRGs, 11 exhibited variations in the expression levels between BLCA tumor specimens and non-tumorous samples. Specifically, an elevation in the expression levels of PRN1, OXSM, SLC3A2, and CD2AP were observed in tumor tissues, while the expression levels of DSTN, FLNA, TLN1, IQGAP1, MYL6, NCKAP1, and NDUFS1 declined in tumor tissues (Fig. 2A). An unsupervised consensus clustering analysis was carried out on the 11 differentially expressed DRGs to ascertain the biological functions of DRGs in BLCA. Results from the analysis revealed that the TCGA-BLCA cohort could be split into two subgroups when k = 2 (Figs. 2B–2D). A notable disparity in the survival rates between the two subgroups was noted, as indicated by Kaplan-Meier survival analysis (Fig. 2E).

Differential analysis and functional enrichment of DEGs

Differential analysis was conducted between the two identified subgroups, and 2,166 DEGs were determined. Among these DEGs, 682 genes were downregulated, while 1,484 were upregulated (Fig. 3A). To elucidate the potential molecular mechanisms of these DEGs, GO and KEGG pathway analyses were implemented. As a result, the GO analysis for biological processes (BP) indicated that those DEGs were predominantly enriched in such processes as cell chemotaxis, arrangement of the extracellular matrix, organization of extracellular structures, and formation of external encapsulating structures. The GO analysis for cellular components (CC) indicated that these genes were predominantly enriched in the collagen-laden extracellular matrix, the outer side of the plasma membrane, and the lumen within the endoplasmic reticulum. With respect to molecular function (MF), the GO analysis demonstrated that the DEGs were principally involved in the glycosaminoglycan interaction, the construction of extracellular matrix structure, as well as endopeptidase activity (Fig. 3B). Prior research has demonstrated the pivotal role of such processes as the arrangement of the extracellular matrix, organization of extracellular structures, collagen-laden extracellular matrix, and the construction of extracellular matrix structure in regulating tumor angiogenesis (Mongiat et al., 2016; Pickup, Mouw & Weaver, 2014; Zhao et al., 2023b). KEGG pathway analysis uncovered substantial enrichment of DEGs within certain pathways, for instance, cytokine-cytokine receptor engagement, the PI3K-Akt signaling route, neuroactive ligand-receptor interplay, as well as the cytoskeleton of muscle cells (Fig. 3C). According to GO and KEGG analyses, 82 ARGs were selected from the molecular signature database and were subsequently intersected with the 2,166 DEGs. Ultimately, 11 common genes were determined, including COL3A1, COL5A2, EPGN, POSTN, SCG2, SERPINA5, SPHK1, SPP1, STAB1, TGFB2, SHH (Fig. 3D).

Formulation and verification of prognostic risk model

To ascertain the implications of the 11 disulfidptosis-related ARGs for prognosis, the LASSO algorithm was employed. Two optimal candidate genes (SCG2 and COL5A2) with the minimum lambda were determined and then were subsequently utilized to formulate the risk model (Figs. 4A and 4B). The formula was established as: Risk Score = COL5A2 * 0.12456 + SCG2 * 0.19130. Patients were stratified into high- and low-risk cohorts according to a specific cut-off of −0.04. Kaplan-Meier (KM) survival analysis indicated that high-risk individuals demonstrated a markedly less favorable prognosis when contrasted with low-risk individuals (Fig. 4C). The ROC curve analysis corroborated the prognostic validity of this model, with AUC readings standing at 0.613 for 3 years and 0.624 for 5 years (Fig. 4D). Scatter plots demonstrated that the high-risk population experienced a heightened mortality rate in comparison to the low-risk population, while heatmaps showed elevated expression levels of the risk-related candidate genes (COL5A2 and SCG2) in the high-risk group (Fig. 4E). To further affirm the model’s robustness, the GEO dataset (GSE32548) was selected for external validation. The results confirmed that patients who had higher risk scores had reduced overall survival (OS) and increased fatality (Fig. 4F). At the 3-year and 5-year time points, the AUC values were 0.803 and 0.762, respectively (Fig. 4G), thereby substantiating the stable predictive efficiency of the risk model.

Association between prognostic risk model and patient clinical features

To elucidate the interrelation between the risk model and the prognosis of BLCA patients, COX regression analysis and univariate and multivariate analyses were implemented on clinical features and risk scores. The results demonstrated that riskScore, Age, and Stage emerged as independent prognostic indicators for BLCA patients (Figs. 5A and 5B). Utilizing the aforementioned risk score along with clinical-pathological elements, specifically, age, gender, as well as disease progression stage, a nomogram was developed for the purpose of determining probabilities of survival at 1, 3, and 5 years (Fig. 5C). The nomogram’s predictive reliability for prognoses at 1, 3, and 5 years was substantiated by calibration graphs (Fig. S1). Furthermore, the linkage between this risk model and BLCA patients’ clinical features within the TCGA cohort was examined. Meanwhile, violin plots revealed that elevated risk scores were noted in older individuals (Fig. 5D), females (Fig. 5E), individuals with higher stages (Fig. 5F), and those with higher pathological TNM (Tumor Node Metastasis) stages (Figs. 5G–5I). In summary, the constructed prognostic risk model exhibited robust performance in projecting the prognosis of BLCA patients.

Association of risk prediction model with immune checkpoints, and tumor 3.5.1 microenvironment

The “CIBERSORT” and “ESTIMATE” algorithms were leveraged to probe into the interplay between this risk prediction model and the tumor microenvironment. It was noteworthy that substantial disparities in immune cell infiltration were detected between these two risk cohorts. The high-risk population manifested a notable decline in several anti-tumor immune cells (Figs. 6A and 6B). Specifically, the proportions of CD8+ T cells, Tregs, as well as activated dendritic cells decreased in the high-risk cohort. The evidence underpinned that tumors in the high-risk population likely possessed an immune-inhibiting microenvironment, impairing tumor-counteracting immune function and consequently promoting tumor progression. The correlation analysis demonstrated that dendritic cells activated, dendritic cells resting and T cells CD8 exhibited an inverse relationship with COL5A2, SCG2, and riskScore (Fig. 6C). This finding aligned with the previously mentioned results of the immune cell infiltration analysis. Between high- and low-risk populations, 16 immune checkpoint genes with significant differential expression were determined, including CD200, LAIR1, and HAVCR2. All 16 genes demonstrated elevated expression levels in the high-risk cohort (Fig. 6D). Moreover, our study indicated that significantly greater stromal, immune, and ESTIMATE scores were found in the high-risk population as opposed to the low-risk population (Fig. 6E). To ascertain the interplay between risk score and TME, TIDE was employed for the purpose of gauging BLCA immunological features in the high- and low-risk cohorts. Both the exclusion score and TIDE score were elevated in the high-risk population (Figs. 6F and 6G), suggesting that BLCA with high-risk scores exhibited enhanced immune escape and potential non-responsiveness to immune checkpoint inhibitors.

Linkage between prognostic risk model and TMB

In an attempt to delve into the somatic mutation profiles between the two risk cohorts, TMB was analyzed. The results demonstrated that the mutation rate among the 198 high-risk patients stood at 93.94% (Fig. 7A), while it was 96.97% for the 198 low-risk patients (Fig. 7B). The frequency of TMB appeared to be reduced in the high-risk population as opposed to the low-risk population. The interplay between TMB and OS among individuals with BLCA was examined. Our findings indicated that patients possessing H-TMB experienced a superior prognosis compared to those with L-TMB (Fig. 7C). Furthermore, patients were subdivided into four groups according to TMB and risk scores. The findings demonstrated that the L-TMB + high-risk group exhibited the poorest prognosis among the four groups (P < 0.0001) (Fig. 7D). This analysis offers insights into the subtle differences between the genetic profiles of BLCA and their associated prognosis.

In vitro validation of two gene expressions

To ascertain the credibility of the two genes incorporated in the risk prediction model, we conducted quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) with standard bladder epithelial cell strain (SV-HUC-1) and bladder cancer cell strains (T24 and SW1710). As illustrated in Figs. 8A and 8B, mRNA expression levels of COL5A2 and SCG2 were considerably elevated in the bladder cancer cell strains (T24 and SW1710) relative to standard bladder epithelial cell strain (SV-HUC-1). The evidence corroborated that the two risk genes could function as diagnostic biomarkers for BLCA.