The mRNA and lncRNA sequencing data of 394 BLCA samples and 19 controls were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). TCGA database was queried up to 24 January, 2024, to retrieve transcriptomic and clinical data (12). The dataset encompassed 413 patients with BLCA, including 394 tumor samples and 19 normal samples. For diverse data analyses, transcriptomic data in HTSeq-Counts and HTSeq-TPM formats were specifically chosen. Samples with incomplete clinical records or overall survival (OS) of <30 days were excluded.

Recent literature has identified several disulfidptosis-related genes (DRGs), including SLC7A11, SLC3A2, SLC2A1, NCKAP1, WASF2 and RAC1(4). A protein-protein interaction network for the six DRGs was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://cn.string-db.org/). Utilizing Pearson correlation analysis of these DRGs, a screening for DRlncRNAs (|R| >0.3, P<0.001) was conducted. Then ‘DESeq2’ package was performed for differential analysis (P<0.05, |log2-fold change| >1) (13).

Using a 7:3 ratio, the 364 BLCA samples were divided into training and testing cohorts. Initially, DRlncRNAs were selected through univariate Cox analysis of the training group (P<0.01). The selection was then refined using the least absolute shrinkage and selection operator (LASSO) algorithm to prevent overfitting of this signature. In the end, 8 DRlncRNAs were selected.

The predictive risk score for each patient with BLCA was calculated using the following formula: Predicted risk score=Σ[(coefficient of lncRNAn) x (expression of lncRNAn)]. Subsequently, patients were stratified into two cohorts based on the median risk score derived from the prognostic signature: High-risk and low-risk groups. To investigate OS differences between these groups in each study cohort, the ‘survival’ and ‘survminer’ software tools were employed to construct Kaplan-Meier (K-M) curves. The ‘timeROC’ package was employed to calculate the area under the receiver operating characteristic (ROC) curves (AUC) to assess predictive accuracy. Additionally, principal component analysis (PCA) was conducted to assess the separation of risk groups (14). Variables that independently impact BLCA survival were identified through univariate and multivariate regression analyses. Finally, K-M survival curves were generated for subgroups with distinct clinical characteristics, providing a comprehensive assessment of the prognostic signature's clinical applicability.

Utilizing the ‘RMS’ package, a prognostic nomogram was constructed for predicting the OS of individual patients with BLCA. This comprehensive model integrates both the risk score and relevant clinical data (15). To validate its precision, ROC curves were employed and Decision Curve Analysis (DCA) was conducted. Furthermore, the nomogram's performance was visually evaluated by plotting calibration curves.

Differentially expressed DRGs and their enrichment in biological signaling pathways were identified through Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. Subsequently, the Gene Set Enrichment Analysis (GSEA) software (version 3.0; Broad Institute website (https://www.gsea-msigdb.org/)) was leveraged to scrutinize GOBP, KEGG and WIKIPATHWAYS gene sets (16). GSEA was conducted on risk-stratified gene expression profiles, utilizing 1,000 resamples, an upper limit of 5,000 genes, and a lower limit of 5. The significance thresholds were set at q-values <0.25 and P-values <0.05.

The ‘TCGAbiolinks’ package was utilized to aggregate somatic mutation data profiles in the format of mutation annotation. Subsequently, a comparative analysis of mutation profiles and TMB scores was conducted between the high- and low-risk groups, using the ‘maftools’ program (17).

The ‘ESTIMATE’ software was utilized to analyze microenvironmental differences between the two subgroups. Additionally, seven algorithms (XCELL, CIBERSORT-ABS, TIMER, EPIC, QUANTISEQ, MCPCOUNTER and CIBERSORT) were used to investigate the correlation between risk scores and infiltrating immune cells (18). To evaluate infiltration scores in the BLCA microenvironment, the ‘GSVA’ tool was applied for single-sample GSEA (ssGSEA).

The tumor immune dysfunction and exclusion (TIDE) platform was utilized to predict the response to immunotherapy in BLCA (19). Furthermore, the ‘oncoPredict’ package was harnessed to assess the IC₅₀ of widely used chemotherapeutic agents and compare their sensitivity across different risk cohorts (20).

The BLCA cell lines (T24 and 5637) (IMMOCELL; (http://www.immocell.com/), where T24 is generally considered to have a higher malignancy level compared with 5637, and the SV-HUC-1 cell line, derived from human normal bladder epithelial cells, were cultivated in a controlled incubator at 37˚C with 5% CO₂, using RPMI-1640 medium (Thermo Fisher Scientific, Inc.). Total RNA was extracted from each sample using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific Inc.). cDNA was reverse-transcribed from the isolated RNA using the PrimeScript RT Reagent Kit (Takara Bio, Inc.), following the manufacturer's protocol. And the SYBR Green premixed qPCR kit (Hunan Accurate Bio-Medical Co., Ltd.) was used in a Roche LightCycler 480 II [Roche Diagnostics (Shanghai) Co., Ltd.]. Subsequently, RT-qPCR was performed. The thermocycling conditions were as follows: initial denaturation at 95˚C for 5 min, followed by 40 amplification cycles consisting of denaturation at 95˚C for 15 sec, annealing at 60˚C for 20 sec, and extension at 72˚C for 30 sec. A final extension step was performed at 72˚C for 5 min. The forward primer sequence for the reference gene GAPDH is 5'-GGAAGCTTGTCATCAATGGAAATC-3', and the reverse primer sequence is 5'-TGATGACCCTTTTGGCTCCC-3'. The relative expression of lncRNAs was quantified as a 2^-ΔΔCq value after assessing gene expression levels via RT-qPCR (21). Experiments were conducted in triplicate, with primer details provided in Table SI. Furthermore, the protein expression profiles of the identified DRGs were examined and compared between normal and BLCA tissues using the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/).

Expression levels of DRlncRNAs in cell lines were detected using RT-qPCR. Statistical significance was determined by unpaired Student's t-test and one-way analysis of variance (ANOVA) to compare the expression levels between different cell lines. Should the ANOVA indicate a significant difference among the groups, the comparison between two groups was conducted using LSD method as the post hoc test. Data are presented as the mean ± standard deviation (SD) from at least three independent experiments. All statistical analyses were conducted using SPSS version 26.0 (IBM Corp.) and RStudio version 4.2.1 [PBC (http://www.rstudio.com/)]. P<0.05 was considered to indicate a statistically significant difference.

The schematic diagram representing the study is depicted in Fig. 1. In the initial phase, following the identification of six DRGs, the STRING database (https://cn.string-db.org/) was utilized to construct a protein-protein interaction network (Fig. 2A). Subsequently, through Pearson correlation analysis, 1,156 lncRNAs correlated with these DRGs were identified. Further differential analysis revealed that 287 lncRNAs exhibited differential expression between cancerous and normal tissues (Fig. 2B).

The 364 patients were divided into two groups: A training group (n=255) and a test group (n=109). Detailed clinical information for each group is provided in Table I. Notably, 14 lncRNAs were identified in the training set through univariate analysis (Table SII). Subsequently, the LASSO algorithm refined this selection, resulting in the identification of 8 DRlncRNAs for the development of the signature (Fig. 2C and D; details in Table SIII). Additionally, the study also explored correlations among these final 8 DRlncRNAs (Fig. S1A) and their upregulation and downregulation patterns (Fig. S1B). Furthermore, the associations between these 8 lncRNAs and their associated genes were investigated (Fig. S1C).

Following the outlined steps, a prognostic signature was formulated for patients with BLCA, and the risk scores were established as specified: Predicted risk score = (-0.003727 x AL390719.2 expression) + (-0.126979 x ASMTL-AS1 expression) + (0.1324342 x AL031058.1 expression) + (0.3540896 x LINC02438 expression) + (0.1858756 x LINC01788 expression) + (0.1228996 x AC022613.2 expression) + (0.1901772 x RBMS3-AS3 expression) + (0.2711116 x AL122035.1 expression). Specific clinical data for both risk groups are provided in Table II.

The prognostic outcomes of the high- and low-risk groups in the training (P<0.0001), test (P=0.03) and entire (P<0.0001) cohorts exhibited significant disparities, as evidenced by K-M curves (Fig. 3A-C). ROC curves revealed anticipated AUC values at different time intervals (1-, 3- and 5-years) of 0.73, 0.70 and 0.70, respectively, in the training cohort (Fig. 3D); 0.80, 0.72 and 0.71, respectively, in the test cohort (Fig. 3E); and 0.75, 0.70 and 0.70, respectively, in the entire cohort (Fig. 3F).

These data affirm the outstanding predictive performance of the signature for patients with BLCA. Heatmaps depicting the expression of 8 DRlncRNAs, risk curves, risk survival status plots (Fig. S2A-C) and scatter dot plots (Fig. S2D-F) in each cohort vividly illustrate the unfavorable survival outcomes among high-risk patients.

Furthermore, PCA underscores the superior discriminatory accuracy of the risk signature when distinguishing between the two groups of patients with BLCA, surpassing the discriminatory power of individual genes, lncRNAs and DRlncRNAs (Fig. 4A-D). Notably, univariate analysis identified age, tumor stage and risk score (all P<0.001) as significant prognostic factors (Fig. 4E). In the multivariate analysis, age, tumor stage and risk score (all P<0.001) independently demonstrated predictive significance (Fig. 4F).

The heatmap visually portrays distinct expression patterns of the 8 lncRNAs alongside clinical features (tumor stage, age and sex) in both high- and low-risk patients (Fig. S3A). Scatter plots revealed that increasing risk scores associated with female sex (P=0.026) and mortality (P<0.001), while age and tumor stage do not exhibit significant associations with the risk score (Fig. S3B-E). Survival curves across diverse clinical subgroups underscore the robust predictive capacity of this signature, particularly among male patients across all tumor stages and age groups (Fig. 5A-F).

Following the aforementioned analyses, a prognostic nomogram that seamlessly integrates the risk score alongside other pertinent clinical features was constructed to predict the OS of patients with BLCA at 1, 3, and 5 years (Fig. 6A). Subsequent validation using data of patients with BLCA substantiated the effectiveness of this predictive tool (Fig. S4).

The AUC further attests to the accuracy of the risk score-based nomogram, with AUC values as follows: Risk score (0.693), age (0.602), sex (0.467) and stage (0.675). Notably, the nomogram incorporating the risk score (0.733) outperformed the nomogram without the risk score (0.703) and DCA confirmed the accuracy of the risk score-based nomogram (Fig. 6B and C). Additionally, calibration curves underscored the superior predictive capacity of the risk score-based nomogram (Fig. 6D and E).

The salient pathways identified through KEGG pathway enrichment analysis are illustrated in Fig. S5A and B. This analysis, conducted across three libraries and examined using GSEA, yielded consistent outcomes. Remarkably, the pathways enriched in high-risk patients were linked to the cell cycle, focal adhesion and WNT signaling (Fig. S6A-C), while those in low-risk patients were predominantly related to substance metabolism (Fig. S6D). Comprehensive details of all pathways are available in Table SIV.

As depicted in Fig. 7A, the low-risk subgroup exhibited an elevated TMB. The waterfall plots visually represented the 20 most frequently mutated genes in patients with BLCA (Fig. 7B and C). Notably, the top four mutated genes among high-risk patients were TP53 (56%), TTN (44%), ARID1A (27%) and KMT2D (25%). Meanwhile, among low-risk patients, the prominent mutated genes were TTN (42%), TP53 (40%), KDM6A (29%) and KMT2D (26%).

In the cohort of high-risk patients, elevated stromal, immune and ESTIMATE scores were observed, alongside diminished tumor purity scores (Fig. 8A-D). These findings suggest a potential link between the unfavorable prognosis in high-risk patients and an immunosuppressive microenvironment that facilitates tumor immune evasion. The accompanying bubble chart delineates the associations of risk scores with immune cell populations (Fig. 8E).

Increased immune cell infiltration was revealed in high-risk patients, as discerned through ssGSEA analysis (Fig. 8F). Notably, immunosuppressive cell subsets [including myeloid-derived suppressor cells (MDSCs), Th2 cells (type 2 T helper cells), and regulatory T cells (Tregs; T follicular helper cells)] were significantly upregulated in the high-risk group. CD56dim natural killer cells and monocytes were more highly infiltrated in the low-risk group. For a comprehensive understanding of the intricate interplay, the detailed immune cell associations are presented in Fig. S7A-D. Overall, the high-risk cohort consistently exhibited augmented immune activity, as evidenced by the ssGSEA analysis of immune-related functions (Fig. 8G). Except for Type II IFN response, which was expressed more strongly in the low-risk group, the other 12 immune-related functions (APC co inhibition, APC co-stimulation, CCR, check-point, cytolytic activity, HLA, inflammation-promoting, MHC class I, para-inflammation, T cell co-inhibition, T cell co-stimulation and type I IFN response) were expressed more strongly in the high-risk group.

In the exploration of therapeutic strategies for BLCA, it has been discerned that low-risk patients exhibit significantly lower TIDE scores, coupled with elevated Microsatellite instability scores, while T cell dysfunction scores remain relatively stable (Fig. S8A-D). These empirical findings posit that low-risk patients with BLCA may manifest reduced vulnerability to tumor immune subversion, thereby potentially augmenting the efficacy of immunotherapy.

Additionally, pharmaceutical agents to which high-risk patients evince heightened sensitivity are listed in Table SV, while drugs that elicit augmented responsiveness in low-risk patients are presented in Table SVI. The pharmacological agents whose sensitivity remains largely unaltered between the two patient cohorts are meticulously catalogued in Table SVII. A summary of antineoplastic drug target pathways is provided in Table SVIII. Noteworthy, compounds targeting the PI3K/mTOR signaling pathway appear to exert a more pronounced impact on high-risk populations, whereas those directed at apoptosis regulation exhibit heightened efficacy against low-risk populations.

HPA, the valuable resource for understanding protein expression patterns, was meticulously explored to discern the protein expression disparities between BLCA and normal samples (Fig. S9). Analyzing DRGs protein expression may reveal their role in BLCA, and abnormal expression could serve as prognostic biomarkers or therapeutic targets. It was found that the proteins encoded by NCKAP1, RAC1, SLC2A1 and SLC3A2 have higher expression levels in tumors, while the protein encoded by WASF2 has higher expression levels in normal tissues. Survival analysis from the HPA indicated that high expression of proteins encoded by NCKAP1, RAC1, SLC2A1 and SLC3A2 is associated with a poor 5-year survival rate in BLCA, whereas WASF2 expression was not statistically significant. Detailed results are provided in Table SIX.

The expression levels of DRlncRNAs in the SV-HUC-1 cell line, which represents a normal bladder epithelial cell line, and the BLCA cell lines (T24 and 5637) were assessed using RT-qPCR (Fig. S10). The results revealed that ASMTL-AS1 expression is higher in the 5637 cell line compared with the T24 cell line, while the expression levels of the other seven genes do not show any statistically significant differences (Fig. S10). Additionally, as demonstrated in Fig. S11, AL390719.2, ASMTL-AS1, AL031058.1 and LINC02438 were upregulated in BLCA cell lines (T24 and 5637), whereas AC022613.2, RBMS3-AS3 and AL122035.1 were downregulated. No significant difference in expression was observed for LINC01788. In addition, to explore the potential roles of DRGs in the two cell lines, bar charts of the expression levels of DRGs in the two cell lines were obtained from HPA. The results indicated that RAC1, SLC7A11 and WASF2 were expressed higher in 5637 cells, while SLC2A1, NCKAP1 and SLC3A2 were expressed higher in T24 cells (Fig. S12).