An overview of our research is shown in Figure 1A. First, we systematically constructed a landscape of the TME immune cell network that comprehensively demonstrated the interactions between immune cells (Supplementary Figure S1A). Then, CIBERSORT algorithms were performed to quantify the infiltration levels of immune cells in UC tissues (Supplementary Table S1). According to the immune cell infiltration data and clinical information of 348 patients (Supplementary Table S2), we performed unsupervised clustering to classify the UC patients into three distinct subtypes (Figure 1B), including 62 patients in TME-ClusterA, 137 patients in TME-ClusterB, and 149 patients in TME-ClusterC (Supplementary Figure S1B). And there were significant differences in prognosis outcomes among these clusters. The TME-ClusterB exhibited a prominent survival advantage, while the prognosis of patients in TME-ClusterA was the worst (log-rank test, p = 0.01, Figure 1C). And the distribution of immune cell infiltration in the IMvigor210 cohort was shown in Figure 1D. In addition, we also performed CIBERSORT analysis in The Cancer Genome Atlas (TCGA) cohort and used the same parameters for consistent clustering. We found that the TCGA cohort could also be divided into three categories, and also had a significant difference in prognosis among three categories (log-rank test, p = 0.00052, Figure 1E), which further indicates the rationality of stratifying urothelial cancer patients according to TME characteristics. Interestingly, we found that in the TCGA cohort, there was a partial overlap in the survival curves between ClusterB and ClusterC, which is most likely due to a batch effect between these two cohorts. And the main conclusion that ClusterB has the best prognosis and ClusterA has the worst prognosis obtained by clustering is not affected.

Immune-associated cells could reflect the characteristics of individual immune microenvironment to a certain extent, and the immune checkpoint is also considered to be an important factor in predicting the response to immunotherapy. The Kaplan-Meier analysis we performed also showed that patients with different levels of immune cell infiltration and immune checkpoint expression had significant difference in clinical prognosis (Supplementary Figure S2). In order to explore the characteristics of patients in different patterns, we carried out a detailed comparison of them. The expression levels of CD8⁺ effector T cells and immune checkpoints in patients of TME-ClusterB were higher than those in patients in the other clusters (p < 0.05) (Figures 2A,B). These results strongly imply that patients in TME-ClusterB may be more likely to benefit from immunotherapy, which is consistent with the favorable prognosis in TME-ClusterB patients. The high level of immune-associated cell infiltration level also indicated that this cluster may be associated with multiple immune-related responses or activities and could be identified as immune-enriched subtype. TME-ClusterA was associated with the activation of epithelial-mesenchymal transition (EMT), the transforming growth factor-β (TGF-β), and Wnt signaling pathways (Figures 2C–I). The expression of specific immune checkpoints was also lower in this cluster (Figure 2J). The patients in this cluster had the worst prognosis, and the infiltration levels of T regulatory cells and M0 and M2 macrophages in this cluster were significantly higher than those in other clusters. Based on these characteristics, this cluster could be identified as the stromal-activation subtype. Interestingly, we also observed abundant immune cell infiltration in TME-ClusterC, such as memory B cells, plasma cells, CD4⁺ memory resting T cells, monocytes, resting dendritic cells, activated dendritic cells, activated mast cells, and eosinophils, but the relative abundance of immune cells did not significantly change the prognosis of these patients, and their powerful antitumor effect was suppressed. So, we defined this group as the immune-suppressive subtype.

To clarify the unique biological role of each cluster in the TME, we performed a differential expression analysis with the limma package. Each cluster was compared with the other clusters in the cohort, and a total of 7,996 DEGs were identified (Supplementary Figure S3A and Supplementary Table S3). Based on these DEGs, the GSVA package was used to analyze the specific enrichment pathways of each cluster (Supplementary Table S4). We found that TME-ClusterA was significantly enriched in HALLMARK EPITHELIAL MESENCHYMAL TRANSITION, HALLMARK COAGULATION, and HALLMARK ANGIOGENESIS, which may be related to the poor prognosis outcome (Supplementary Figure S3B). We conducted a functional enrichment analysis by the clusterProfiler R package (Supplementary Table S5) and found enrichment mainly in neutrophil activation, neutrophil-mediated immunity, T cell activation, regulation of innate immune response, and other immune-related Gene Ontology (GO) terms (Figure 2K). This once again proved the close relationship between the DEGs and immune-related functions.

To further explore the association between the DEGs and phenotypes, we conducted another unsupervised clustering analysis (Supplementary Figure S3C) and found that the cohort could also be divided into three cohorts with significant prognosis differences (log-rank p = 0.045) (Figure 2L). We referred to these as Gene-Clusters A, B, and C. Patients in Gene-ClusterB had the best prognostic outcomes, while patients in Gene-ClusterA had the worst. Using a chi-square analysis to compare the gene clusters and TME patterns, we found good consistency between these two grouping methods (χ² contingency tests, p < 2.2e-16). The distribution of patients by TME patterns and TME gene clusters are shown and specific information can be found in Supplementary Table S6. There were also significant differences in the level of immune checkpoint expression (Figure 2M) and immune cell infiltration (Figure 2N) among these gene clusters, indicating that these gene clusters could also represent the TME characteristics of the patients.

For the TME-related DEGs, we first matched expression data with clinical information, then reduced the dimension by using the univariate Cox regression model and used the more stringent p < 0.01 as the screening criterion to select 318 prognosis-related genes for further analysis (Supplementary Table S7). Next, we divided the IMvigor210 cohort into a training set and testing set at a ratio of 7:3. In the training set (n = 244), we performed modified Lasso regression analysis to construct the TMSig. In the process of cyclic calculation, we found that the maximum AUC value of the TMSig at 2 years was 0.906 (Supplementary Figure S4A). We defined the gene signature present at this time as the best candidate model. The prognosis of the low-score group was significantly better than that of the high-score group in the training set (log-rank p < 0.0001) (Figure 3A). Similar to the results obtained with the training set, the low-score group had a better prognosis in the internal testing set (log-rank p < 0.0001) (Figure 3B) and in the entire IMvigor210 cohort (log-rank p < 0.0001) (Figure 3C). The ROC curves proved the robust predictive ability of the TMSig, and the AUC at 1 year was 0.88 and 0.906 at 2 years in the training set. The AUC at 1 year was 0.747, and that at 2 years was 0.805 in the testing cohort. In the whole cohort, the AUC was 0.840 at 1 year and 0.876 at 2 years (Supplementary Figures S4B–D). Then we performed the univariate Cox regression algorithm to analyze TMSig together with other clinical characteristics of the patients in the training set. And further included them in the multivariate Cox regression algorithm after screening out the features with p < 0.05. And the same analysis was performed not only in the train cohort, but also in the test cohort and the entire IMvigor210 cohort. The p-value of TMSig < 0.05 in each time of analysis, proving that it can be served as an independent prognostic factor for patients (Figure 3D). And the TMSig also showed better prognostic predictive power in three independent cohorts (Figures 3E–G).

The TMSig score had a significantly different distribution in BOR, immune phenotype, IC level and other subgroups (Kruskal−Wallis, p < 0.05) and had a difference in TC level, but it was not significant (Supplementary Figures S4E–H). To further verify the reliability of the TMSig for predicting the prognostic outcomes, we performed a stratified analysis based on the clinical information of the IMvigor210 cohort. Through the Kaplan-Meier analysis, we found the TMSig has a great performance in several clinical subgroups (Immune phenotype: immune desert type, immune excluded type, immune inflamed type, IC level: IC0, IC1, IC2, Sex: Male, Female, BOR: SD/PD, CR/PR, TC level: TC0, TC1, TC2, Tobacco history: NEVER, PREVIOUS OR CURRENT, Supplementary Figures S4I–W). We also conducted an external verification of the prognostic value of the TMSig in the independent TCGA-BLCA cohort and found that it was of great significance for predicting both the overall survival (OS) and disease-specific survival (DSS) rates of these patients (Supplementary Figures S5A,B). In addition, stratified analysis of the TCGA cohort showed that the TMSig had significant prognostic significance in patients with higher disease stages (Supplementary Figures S5C,D), which inspired us to conclude that the TMSig may play a unique role in predicting the prognosis of patients with advanced neoplasia.

To improve the clinical application of TMSig, we constructed dynamic nomogram (TMSigDynNomapp: https://the-nomogram.shinyapps.io/TMSigDynNomapp/, Figure 3H), while calibration plots showed that comprehensive signature has accurate predictive power at different time points (Figure 3I). Decision curve analysis also showed that comprehensive signature can provide better clinical benefit to patients compared to applying gender, IC level, and other indicators for prediction (Figure 3J). Compared with other previously reported bladder cancer-related signatures, TMSig also has more robust predictive power (Figure 3K).

We used ROC curves to evaluate the ability of the TMSig score to predict the efficacy of immunotherapy among patients in the IMvigor210 cohort and compared the score with known effective predictors such as TMB (Samstein et al., 2019), TNB (Wolf et al., 2019), and M1 macrophages (Zeng et al., 2020). It was found that the accuracy of the TMSig in effectively predicting the response to immunotherapy was not inferior to that of other biomarkers. (TMSig score AUC: 0.826, TMB AUC: 0.728, TNB AUC: 0.767, M1 macrophage AUC: 0.702) (Supplementary Figure S5E). To fully demonstrate the robustness of the TMSig for predicting immunotherapy response, we included two independent data sets for external validation. The AUC predicted by the TMSig score in the data from Miao et al. was 0.75 (Supplementary Figure S5F), and the AUC predicted by the TMSig score in the GSE35640 dataset was 0.687 (Supplementary Figure S5G). All these results indicate the great potential of the TMSig for discriminating patients who may benefit from immunotherapy.

It is well known that patients with different infiltration levels of immune cells have different prognostic outcomes or treatment responses. Therefore, we performed a Spearman correlation analysis to explore the relationship between the TMSig score and the infiltration level of various immune cells and found that there was a significant positive correlation between the TMSig score and M0 macrophages, resting mast cells, neutrophils, and eosinophils (p < 0.05, cor > 0) and a significant negative correlation with the infiltration level of follicular helper T cells, activated NK cells, gamma delta T cells, memory B cells, CD4⁺ memory-activated T cells, and M1 macrophages (p < 0.05, cor < 0) (Figure 4A). The strongest positive correlation was between the TMSig score and M0 macrophages (Figure 4B), and the strongest negative correlation was between the TMSig score and M1 macrophages (Figure 4C). There was a moderate but significant negative correlation between the TMSig score and TNB (Kruskal-Wallis, p = 0.00049) (Figure 4D), and the correlation between the TMSig score and TMB also demonstrated the same trend (Figure 4E). Intriguingly, combining the TMSig score with TMB or TNB contributed to the survival assessment (Kaplan-Meier analysis, TMSig score +TNB binary: p < 0.0001; TMSig score +TMB binary: p < 0.0001) (Figures 4F,G). We should clear that the correlation between TMSig and M1 macrophages is strong and deserves focused attention. But its correlation with TMB or TNB is moderate, which could provide direction for our study, but the exact relationship needs to be verified by further studies.

The high-score group had a higher TMSig score (Wilcoxon, p < 2.2e−16) (Figure 5A), indicating that the two groups have unique features not only in prognosis but also in immune-related characteristics. TIGIT, CD274, CTLA4, PDCD1, and LAG3 presented higher expression levels in the low score group (p < 0.05) (Figure 5B), which suggests that people with low TMSig scores might have a better response to immunotherapies targeting immune checkpoints. The immune infiltration cell analysis also showed that these two groups had significantly different marker immune cells (Figure 5C). Then, we used the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to evaluate each patient’s potential response to immunotherapy and observed that the responsiveness of immunotherapy in the low-score group was higher than that in the high-score group (p = 0.02) (Figure 5D). Moreover, subclass mapping was performed with another group of 47 melanoma patients who responded to immunotherapy (Roh et al., 2017). We were encouraged by the observation that a low score indicated potential patient response to PD-1 treatment. (Bonferroni corrected p = 0.008) (Figure 5E). These results reconfirmed the application value of the TMSig. To further explore the significantly enriched pathways of the DEGs between the two groups, we carried out GSEA using the clusterProfiler and fgsea R packages. It was found that HALLMARK ANGIOGENESIS, HALLMARK TGF-ß_ SIGNALING, HALLMARK APOPTOSIS, HALLMARK HYPOXIA, and HALLMARK P53 PATHWAY were significantly enriched in the upregulated genes (Supplementary Table S8), which may be related to poor prognosis.

In addition, we found that the high-score samples of the Baylor subtype and UNC subtype were mainly concentrated in the basal type (Wilcoxon, Baylor subtype, p = 2.1e-04; UNC subtype, p = 2.8e-04) (Figures 5F,G), and the TMSig score in the basal MDA subtype was also significantly higher than that in the luminal type and p53-like type (Kruskal−Wallis, p = 3e−06) (Figure 5H). In the TCGA subtypes, high scores were mainly distributed in clusters III and IV (Figure 5I), and it is generally believed that clusters I/II and III/IV in the molecular subtypes officially obtained by the TCGA are similar to differentiated (or luminal) and basal tumors, respectively (Mo et al., 2018). We also verified the distribution of TMSig scores among the four classification methods in the independent TCGA-BLCA dataset, and the results were consistent with those of the IMvigor210 cohort (Figures 5J–M). Therefore, we inferred that the TMSig is closely related to molecular subtypes of bladder cancer, and different scores may indicate different molecular characteristics, which is of great significance for further understanding the characteristics of these two groups.

The CTRP and PRISM datasets contain gene expression profiles and drug sensitivity profiles for hundreds of CCLs and can be used to construct predictive models of drug response. After removing duplicate drugs, these two datasets share 168 compounds, for a total of 1752 compounds (Figure 6A). We removed drugs with deletion values greater than 20% and cell lines derived from haematopoietic and lymphoid tissue. Finally, 680 CCLs for 354 compounds in the CTRP dataset and 480 CCLs for 1285 compounds in the PRISM dataset were used for subsequent analyses. The specific screening process is shown in Figure 6B. Before proceeding further, we first demonstrated that the results of drug response estimation are reliable. Cisplatin is a common therapeutic agent for bladder cancer patients, and a recent study showed that high GULP1 expression enhanced the sensitivity of patients to cisplatin (Teramoto et al., 2021). We divided the patients into high and low expression groups according to the expression level of GULP1. The Wilcoxon rank sum test was used to compare the difference in AUC estimates of cisplatin between the two groups, and the results showed that the AUC estimates were significantly higher (p = 0.003) in patients with high GULP1 expression (Figure 6C), consistent with the clinical presentation of cisplatin. After verifying the reliability of the calculation method, we adopted a similar analysis method to Yang et al. (2021). First, differential drug response analysis was performed between the group with high TMSig score (upper decile) and the group with low TMSig score (lower decile) to identify the group with high TMSig score (log2FC > 0.10) with low estimated AUC values. Then, by Spearman correlation analysis between AUC value and TMSig score, compounds with negative correlation coefficients (Spearman’s r < −0.30 for CTRP or <0.45 for PRISM). These analyses yielded one CTRP-derived compound (PD318088) and two Prism-derived compounds (Levocarnitine, YM−976) (Figures 6D,E). Secondly, the fold-change difference of the expression level of candidate drug target genes between tumor tissues and normal tissues (including paired analysis and unpaired analysis) was calculated. A higher fold change value indicated a greater potential of candidate agent for UC treatment (PD318088: MAP2K1, MAP2K2; YM−976: PDE4B, PDE4D) (Figures 6F,G). Finally, we searched at PubMed (https://www.ncbi.nlm.nih.gov/PubMed/) to find evidence of candidate compounds for UC treatment. Overall, PD318088 and YM−976, with relatively sufficient evidence, are considered to be the most promising potential treatment drugs for people with high TMSig score.