Twenty-three pairs of frozen BLCA tissues and matched adjacent samples were obtained from Qilu Hospital (Qingdao) with all patients’ informed consent. All samples were used to detect the mRNA level of characteristic genes by qPCR, and six pairs of samples were used to detect the protein expression level of characteristic genes by IHC. This study was approved by the Ethics Committee of Qilu Hospital (Qingdao) of Shandong University. The clinical data of BLCA patients are shown in Supplementary Table S1 .

Normalized The Cancer Genome Atlas (TCGA)-BLCA gene expression data (TPM) and clinical data were downloaded from the UCSC Xena platform (http://xena.ucsc.edu/, accessed on 15 December 2022) and used as the training dataset for constructing the TEX model, which included 411 tumor samples and 19 paracancerous samples. The gene expression data (TPM) of normal bladder tissue in the GTEx database were also downloaded from the UCSC Xena platform, which included nine normal samples. The expression data and survival information of the GSE13507, GSE48276, and GSE19915 datasets were downloaded from the Gene Expression Omnibus (GEO) website (https://www.ncbi.nlm.nih.gov/geo/, accessed on 25 December 2022) as the validation datasets, where multiple probes corresponding to the same gene, the probe values were averaged.

Referring to the specific pathways related to TEX reported in previous studies (24), including tumor necrosis factor (TNF), interleukin (IL)-2, interferon (IFN)-γ, and T- cell cytotoxic pathways. The gene sets of the above pathways were downloaded from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb, accessed on 15 December 2022).

The single sample Gene Set Enrichment Analysis (ssGSEA) algorithm in the R package GSVA was used to calculate the enrichment scores of the gene sets of the above four pathways in the TCGA-BLCA samples. The R package survival performed univariate analysis to calculate the correlation between the enrichment scores and the prognosis of BLCA patients, including overall survival (OS) and progression-free interval (PFI). Wilcox test was used to calculate the correlation between the enrichment scores and the clinical features of BLCA patients, including age, gender, tumor grade, tumor stage, lymphatic invasion, and therapeutic effect. The mutation data of the genes in the above four pathways were downloaded from the TCGA database, and R package maftools drew the mutation map of the top 50 genes. The key modules related to the enrichment scores were screened using the weighted correlation network analysis (WGCNA). Genes in key modules were identified as candidate markers reflecting the TEX signature. Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed on these candidate markers through the R package clusterProfiler.

Based on the candidate genes screened above, markers whose expression was significantly associated with the prognosis of BLCA patients were further screened by univariate Cox regression analysis of R package survival and visualized by R package survminer. The prognosis-related markers were analyzed by lasso Cox regression using R package glmnet. The genes with non-zero regression coefficients were obtained as characteristic genes to construct a prognosis-related TEX model (Hereinafter referred to as the TEX model).

The risk score for each BLCA sample was calculated using the TEX model, the formula is A=∑coef gene i × expr gene i, where coef is the regression coefficient of the gene i in the model, and expr is the expression value of the gene i. The BLCA patients were divided into high- and low-risk groups based on the median risk score of all samples. The prognostic predictive efficacy of the constructed TEX model was analyzed in TCGA-BLCA datasets and external independent datasets. The R package maftools was used to visualize the gene mutations between the high- and low-risk groups. Differences in gene expression and clinical parameters between the two groups were also analyzed.

Based on the clinical features of BLCA patients in the TCGA-BLCA datasets, the Wilcox test of R package ggpubr was used to analyze the differences in the risk score among different groups, including age, gender, tumor grade, tumor stage, lymphatic invasion, etc. Univariate and multivariate Cox regression analyses of R package survival were used to observe whether the TEX model was an independent prognostic predictor of BLCA patients. In addition, a nomogram consisting of the TEX model with commonly used clinical features was constructed by R package rms to assess the model’s efficacy in predicting the progression of BLCA patients.

Based on the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/, accessed on 20 January 2023), the Wilcox test of R package ggpubr was used to analyze the difference in chemotherapeutic drug sensitivity between high- and low-risk groups, and the differential expression of characteristic genes in the TEX model. R package Corrplot was used to analyze the correlation between the expression of characteristic genes and chemotherapeutic drug sensitivity.

The TEX model’s efficacy in predicting immunotherapy’s effect in BLCA patients was evaluated by three visualization forms: KM curves drawn by R package survival, scatter plots drawn by R package ggpubr, and ROC curves drawn by R package pROC. Referring to previous reports (25), the TIP tool (http://biocc.hrbmu.edu.cn/TIP/, accessed on 30 January 2023) was used to analyze the anti-cancer immune process activity between the high- and low-risk groups of TCGA-BLCA patients, and the Wilcox test of R package ggpubr analyzed the difference.

Total RNA was extracted from the tumor and matched adjacent tissue samples using TRIzol LS reagent (Thermo Fisher Scientific, USA). The cDNA was synthesized by Reverse Transcriptase M-MLV (Takara, Japan) following the manufacturer’s instructions. In addition, qRT-PCR was performed with SYBR^® Premix Ex Taq™ II (Takara, Japan). GAPDH was used as an endogenous control for qPCR. All primers of twelve characteristic genes of the TEX model were listed in Supplementary Table S2 . The relative expression was compared using 2−^ΔCt between tumor and adjacent samples, with ΔCt = Ct_{characteristic genes} – Ct_GAPDH.

For IHC, 5-micron frozen sections were fixed in pre-cooled acetone. After washing and antigen retrieval, sections were placed in 3% hydrogen peroxide to block endogenous peroxidase. After antigen blocking, primary and secondary antibodies were added for incubation. The specific information on the antibodies used is shown in Table S3 . After DAB color development and cell nucleus restraining, dehydration and mounting were performed. Protein expression levels were analyzed by Image-Pro Plus software (Media Cybernetics, USA).

Analysis of the correlation between the enrichment scores of four specific pathways related to TEX in each TCGA-BLCA sample and patients’OS revealed that the high enrichment of IL-2 signaling and IFN-γ signaling was significantly positively correlated with the OS of BLCA patients (P<0.05, Figure 1A ). BLCA patients were classified into high and low groups based on the enrichment of each pathway using the optimal threshold. The results showed that the enrichment scores of T-cytotoxic pathway, IFN-γ signaling, and IL-2 signaling significantly affected the PFI of BLCA patients. Similarly, the patients with lower enrichment scores had shorter PFI time (P<0.01, Figure 1B ). Interestingly, the difference analysis of the enrichment scores of the four pathways among different clinical features showed that IFN-γ signaling and T-cytotoxic pathway were significantly highly enriched in samples with higher tumor grades (P<0.05), the other two pathways also showed this trend ( Figure 1C ), suggesting that the specific pathways related to TEX may become more abnormal in patients with a higher malignant degree. In addition, the enrichment score of TNF signaling was negatively correlated with lymph node metastasis in BLCA patients (P<0.05, Figure 1D ). There was also a specific correlation between the enrichment scores of the four pathways and the therapeutic effect. The enrichment scores of patients with complete response (CR) and progressive disease (PD) were significantly higher than those of patients with stable disease (SD), suggesting that the activity of the pathways involved has an important impact on the development of the disease. There may be more significant heterogeneity among BLCA patients ( Figure 1E ).

The gene set of the four pathways includes a total of 160 genes, of which 147 genes have obtained mutation data through the TCGA database. The mutation map of the top 50 genes with the highest mutation rate is shown in Figure 2A . The results showed that missense mutations dominated the mutation type.

Key modules associated with enrichment scores were screened using WGCNA analysis. According to the soft threshold and the scale-free fitting index, it was found that there was no intersection point between power and the horizontal line (correlation 0.8) ( Figure 2B ). Referring to the previously reported method (26), the estimated value of powerEstimate equal to 5 is selected as the β value. It can be seen that k is negatively correlated with p(k) (correlation coefficient 0.89>0.8) ( Figure 2C ), indicating that the selected β value is suitable for the establishment of a gene scale-free network. Six modules were obtained ( Figure 2D ), and two modules, Meyellow with a significant positive correlation and Meturquoise with a significant negative correlation, were selected according to the correlation between modules and enrichment scores. A total of 4571 genes were included as candidate markers for TEX ( Figure 2E ). The R package clusterProfiler was used to perform GO and KEGG analysis on the obtained 4571 genes, and the top 20 pathways were selected for display ( Figure 2F ).

Based on the above 4571 candidate marker for TEX, univariate Cox was used to screen OS-related genes, and a total of 147 genes were found with P<0.001 as a threshold. Based on the 147 prognostic-related TEX markers, according to the previous method (27), 28 characteristic genes with non-zero coefficients were obtained for constructing the prognosis-related TEX model using LASSO Cox analysis of R package glmnet ( Figures 3A–C ).

Based on the median risk score of the TEX model, the TCGA-BLCA samples were divided into two groups, namely TEX^high (n=201) and TEX^low (n=202) ( Table S4 ). Survival analysis showed that BLCA patients in the TEX^high group had a shorter OS (P<0.001, Figures 3D, E ). The ROC curve results showed that the TEX model’s AUCs for predicting patients’ survival at 1, 3, and 5 years were 0.767, 0.773, and 0.773, respectively, which were all greater than 0.75. The TEX model has high accuracy in predicting the prognosis of BLCA patients ( Figure 3F ).

The effect of the TEX model on the prognosis of BLCA patients was validated in multiple independent external datasets. Based on the best cut-off values of TEX model risk scores for all samples, BLCA patients in GSE13507 dataset were divided into TEX^high and TEX^low groups. The results showed that patients in the TEX^high group had shorter OS (P = 0.026, Figure S1A ). Similar results were obtained in the GSE48276 and GSE19915 datasets (P = 0.049 and P = 0.012, Figures S1B, C ).

The patients in the training dataset TCGA_BLCA were stratified according to different clinical characteristics. The results showed that the risk scores of the TEX model were closely related to age, tumor stage, tumor grade, lymph node metastasis, and treatment effect. BLCA patients with older age, higher tumor stage or grade, positive lymph nodes, and SD or PD after treatment had higher risk scores ( Figures 4A ; S2 ).

Univariate and multivariate Cox regressions were used to analyze the prognostic correlation between the TEX model and the clinical features of the training dataset TCGA_BLCA. The results of both univariate analysis (P < 0.001, HR =3.44, 95%CI = 2.48-4.77) and multivariate analysis (P < 0.001, HR=3.15, 95%CI = 2.24-4.43) indicated that TEX model could be used as an independent risk factor for BLCA patients ( Figure 4B ). The univariate and multivariate Cox regression analysis in the external validation dataset GSE13507 also showed that the TEX model could be an independent risk factor for BLCA patients ( Figure S3 ). Using the TEX model and multiple clinical parameters, including age, stage, grade, gender, etc., to construct a nomogram to predict the prognosis of BLCA patients in 1, 3, and 5 years ( Figure 4C ). Survival probability prediction plots indicate that the constructed nomogram exhibited a strong prognostic predictive ability ( Figure 4D ).

Wilcox test was used to analyze the difference in sensitivity of chemotherapeutic drugs between the TEX^high and TEX^low groups. The results showed that cisplatin, gemcitabine, gefitinib, camptothecin, and afatinib had higher estimated IC50 values in the TEX^high group ( Figure 5A ). When the 25 characteristic genes used to construct the TEX model were divided into high- and low-expression groups according to their median expression values, it was found that the expression of 21 genes significantly impacted the sensitivity of chemotherapy drugs ( Figure S4 ). The correlation analysis between the expression levels of 25 characteristic genes and the sensitivity of chemotherapeutic drugs intuitively showed that the high expression of epidermal growth factor receptor (EGFR), epithelial membrane protein 1 (EMP1), and EH domain binding protein 1 (EHBP1) was positively correlated with higher IC50 value. In contrast, the increased expression of transcription elongation factor A N-Terminal and central domain containing (TCEANC), gasdermin B (GSDMB), cation channel sperm associated 4 (CATSPER4) and CLK4 associating serine/arginine rich protein (CLASRP) was related to the lower IC50 value ( Figure 5B ).

Immunotherapy is one of the important methods for the clinical treatment of BLCA patients. The analysis results of the prognosis-related data of BLCA patients receiving immunotherapy showed that the SD/PD group had a higher risk score than the CR/PR group ( Figure 5C ). The AUC value of the ROC curve was 0.7 ( Figure 5D ), so the TEX model can better predict the effect of immunotherapy in BLCA patients. The difference analysis of anti-cancer immune process activity between high- and low-risk groups showed that the recruitment process of CD8+ T cells (no significant difference), macrophages, neutrophils, NK cells, and Myeloid-derived suppressor cells (MDSCs) were more active in TEX^high group. In contrast, Th2 and regulatory T (Treg) cells were more involved in the TEX^low group ( Figure 5E ).

We analyzed the 20 most frequently mutated genes in the high- and low-risk groups to explore the underlying molecular mutations behind the differences in immune activity processes. The results showed that the mutation rate of tumor suppressor gene was higher in the TEX^high group, such as tumor protein P53 (TP53) (55% vs 43%), AT-rich interaction domain 1A (ARID1A) (28% vs 22%) and lysine methyltransferase 2D (KMT2D) (27% vs 26%). It is particularly striking that the mutation rate of RB transcriptional corepressor 1 (RB1) was as high as 22% in the TEX^high group while less than 14% in the TEX^low group ( Figure 5F ). RB1 is an immune-related prognostic biomarker and a promising target for immunotherapy in several cancers, including BLCA (28, 29).

In view of the importance of the 28 characteristic genes in model construction, we analyzed their expression differences between the TEX^high and TEX^low groups, and the results showed that nine genes were highly expressed in the TEX^high group, while 19 were highly expressed in the TEX^low group ( Figure S5 ). To facilitate the use of the TEX model in clinical practice, we alternatively performed the differential expression analysis of 28 characteristic genes in TCGA and GTEx datasets, including 407 BLCA samples and 28 paracancerous or normal bladder tissue samples. The results showed that Charged Multivesicular Body Protein 4C (CHMP4C), Signal Transducing Adaptor Family Member 2 (STAP2), Prickle Planar Cell Polarity Protein 3 (PRICKLE3), Zinc Finger Protein 165 (ZNF165), Gasdermin B (GSDMB) and SH2 Domain Containing 2A (SH2D2A) were highly expressed in bladder cancer. In contrast, Amylase Alpha 2B (AMY2B), EH Domain Binding Protein 1 (EHBP1), Epithelial Membrane Protein 1 (EMP1), Maestro (MRO), Fibronectin Leucine Rich Transmembrane Protein 2 (FLRT2) and Procollagen C-Endopeptidase Enhancer 2 (PCOLCE2) were low expressed ( Figure S6A ). We then performed differential expression analysis on 19 paired samples (BLCA tissues and paired paracancerous tissues) in the TCGA dataset. The results showed a similar expression trend ( Figure S6B ). To further verify the results, we detected the above 12 characteristic genes by qPCR in 23 pairs of BLCA and paracarcinoma samples, and the results were consistent with the database analysis ( Figure 6A ). Six characteristic genes highly expressed in BLCA comparing to paracarcinoma were verified by IHC ( Figures 6B ; S7 ).