Transcription profiles and clinical information of BLCA patients stored in TCGA were downloaded from the TCGA-GDC website (https://portal.gdc.cancer.gov/). The reference file “Homo_sapiens.GRCh38.gtf obtained from Ensembl website (http://asia.ensembl.org/index.html) was used for gene annotation. GEO database (https://www.ncbi.nlm.nih.gov/gds/?term=) was comprehensively retrieved with the following terms: 1. searching terms were “(bladder) AND (cancer OR carcinoma)”; 2. the entry type was “series”; 3. the organism was “Homo sapiens”; and 4. the expression profile array was sequenced in tissues and the number of samples was larger than 50. Finally, seven individual GEO datasets were identified, including GSE120736 (platform: GPL10558), GSE128959 (platform: GPL6244), GSE13507 (platform: GPL6102), GSE31684 (platform: GPL570), GSE32548 (platform: GPL6947), GSE48075 (platform: GPL6947), and GSE83586 (platform: GPL16570). TCGA sequence data with FPKM form were converted into TPM form to better match microarray data and the code was obtained from https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/. Sva package in R environment was used to perform data combining and magnitude harmonization, which resulted in a range between 1 and 20 of the expression value of each cohort (Leek et al., 2012). In detail, the ComBat function in sva package was used to remove batch effects and improve reproducibility. IMvigor210 cohort, a urothelial carcinoma cohort treated with the anti-PD-L1 antibody atezolizumab, was used for the prediction of patient response to immunotherapy, which was downloaded under Creative Commons 3.0 license. Differentially expressed genes (DEGs) were identified using limma package with a threshold of |logFC| > 1 and adj.p < 0.05 (Ritchie et al., 2015).

The ClusterProfiler package was used to conduct gene set enrichment analysis (GSEA), which provided a gene classification method for pathway enrichment (Yu et al., 2012). The Hallmark gene set (MSigDB) was selected as a reference set to explore the difference in the oncogenetic pathways. Single sample GSEA (ssGSEA) was performed using the package GSVA in an R environment to evaluate the enrichment scores of 53 immune terms. The gene set used to quantify 53 immune terms has been uploaded in Figshare (https://figshare.com/s/d94792fd1f413eb68ca3) (Ren et al., 2021).

In the combined cohort that has been quantified into 53 immune terms, univariate Cox analysis was first performed to identify prognosis-related terms with the threshold of p-value <0.05. Random forest algorithm was then used for dimension reduction (ntree = 1,000), which was performed in R software using the randomForestSRC package (Alaa et al., 2019). The top 10 important terms were then identified for multivariate cox analysis and IBS calculation using the following formula: riskscore = coef1*term1 + coef2*term2 + coef3*term3 + … + coefN*termN. The patients were divided into high- and low-risk groups based on the median riskscore. ROC curve and Kaplan–Meier survival curve were used to evaluate the prognosis value of riskscore.

TCGA-BLCA tumor mutation data were downloaded from the cBioPortal website (http://www.cbioportal.org/datasets - Bladder Urothelial Carcinoma, TCGA, PanCancer Atlas). R software was used to take out the same sample as the expression data, including “Tumor_Sample_Barcode”, “Hugo_Symbol”, “Chromosome”, “Start_Position”, “End_Position”, “Variant_Classification”, “Variant_Type”, “Reference_Allele”, “Tumor_Seq_Allele1”, and “Tumor_Seq_Allele2”. The mutation types were distinguished into nonsynonymous and synonymous mutations. The tumor immune dysfunction and exclusion (TIDE) score was calculated based on the gene expression data, which could partly predict the immunotherapy response rate of enrolled patients (http://tide.dfci.harvard.edu/login/) (Fu et al., 2020). In addition, the result of comprehensive immunogenomic analyses was obtained from The Cancer Immunome Database (TCIA) website (https://tcia.at/home).

The connectivity map (CMAP) is a gene expression profile database based on intervention gene expression developed by the Broad Institute. It is mainly used to reveal the functional relationship between small molecular compounds, genes, and disease status. The top 500 upregulated and downregulated genes between high- and low-risk patients were identified as input files of CMAP analysis. The small bioactive molecules that were potential targets for BLCA immune signature were then screened based on these DEGs. ComplexHeatmap package was used for result visualization (Gu et al., 2016).

All the analyses were performed using R software (R Core Team, version 4.0.0) and GraphPad Prism version 7.0. All statistical tests were two-sided, and p-value <0.05 was considered statistically significant. Continuous variables with normal distribution in two groups were compared using an independent t-test. Continuous variables with skewed distribution in two groups were compared with Wilcoxon rank-sum test. Kruskal–Wallis test was used to compare the data in more than two groups. The log-rank test was used to compare the statistical difference in the Kaplan–Meier survival curve. The ggplot2 package was used to perform most graphic drawings (https://ggplot2.tidyverse.org/). The R code and primary data in this study were uploaded in Figshare (https://figshare.com/s/d94792fd1f413eb68ca3).

The flowchart of the whole study is shown in Figure 1. After comprehensively retrieving the GEO and TCGA database, eight BLCA cohorts met our criteria, including GSE120736, GSE128959, GSE13507, GSE31684, GSE32548, GSE48075, GSE83586, and TCGA database (Figure 2A). The detailed clinicopathological characteristics of patients in the TCGA cohort and each GEO cohort have been summarized in Supplementary Table S1–8. Meanwhile, we found a lot of inter-assay variances in these datasets (Figure 2A; Comp 1: 78.7% variance, Comp 2: 6.6% variance). Therefore, we used sva package to eliminate potential inter-assay differences as much as possible and to generate a combined BLCA cohort (Figure 2B; 1,614 samples). Furthermore, we used the ssGSEA package to quantify patients’ gene expression profiles into 53 immune terms for further analysis (Figure 2C).

Patients with complete survival information, including survival time and status, were selected for further analysis to explore the effect of identified immune terms on patients’ prognosis. First, the patients in combined cohorts with complete prognostic information were randomly divided into training and validation groups at a ratio of 1:1. Then, we performed univariate Cox analysis to identify prognosis-related terms. The result showed that 10 terms were protective factors, and six were risk factors for the BLCA patients (Figure 3A). Co-expression relationships of these prognosis-related terms are shown in Figure 3B. Next, we performed a random forest algorithm to calculate the relative importance of these prognosis-related terms (Supplementary Figure S1, oob error rate = 0.373 when ntree = 1,000). The top 10 important terms (CD8_T_cells, core serum response up, CSR_activated, cytotoxic cells, lymph vessels, lymphs_PCA, macrophages, NK cells, T_cell_receptors_score, and Tfh cells) were selected for further multivariate cox analyses (Figure 3C). Then, a prognosis model based on these immune terms was established with the formula of “riskscore = CD8_T_cells * -0.1168 + Tfh cells * -0.1346 + lymph vessels * -0.0163 + macrophages * 0.3199 + CSR_activated * 0.0280 + cytotoxic cells * -0.2001 + T_cell_receptors score * -0.1205 + core serum response up * 0.1471 + NK cells * 0.0745 + Lymphs_PCA * -0.1213. The risk plot of patients in the training cohort is shown in Figure 3D. The Kaplan–Meier survival curve indicated that the patients in the high-risk group might have worse overall survival than patients in the low-risk group (Figure 3E; p < 0.0001). In parallel, the ROC curve showed a satisfactory predicted efficacy on patients’ prognosis in the training group (Figure 3F; 1-year AUC = 0.773, 3-year AUC = 0.767, and 5-year AUC = 0.756). Also, our prognosis model showed an effectively predicted efficacy in the internal validation group (Figures 3G–I; 1-year AUC = 0.763, 3-year AUC = 0.700, and 5-year AUC = 0.704) and IMvigor210 external cohort (Figure 3J–L; 1-year AUC = 0.638, 3-year AUC = 0.682, and 5-year AUC = 0.658). Univariate and multivariate Cox regression analysis showed that our model has satisfactory prognostic efficiency independent of clinical factors (Supplementary Figure S2).

Clinical and biological relevance were explored to probe into the possible mechanism of our model on patients’ prognosis. We found that high-risk patients might have adverse baseline characteristics (Figure 4A). Moreover, patients with high grade tend to have higher riskscore than patients with low grade (Figure 4B). Interestingly, we found that riskscore was not associated with disease relapse but with disease progression (Figure 4C). In addition, we explored the riskscore difference between non-muscle invasive and muscle invasive BLCA and found that muscle invasive patients had higher riskscore values (Figure 4D). Estimate package quantified tumor microenvironment into immunescore, stromalscore, and estimatescore, standing for the content of immune, stromal, and non-tumor cells, respectively (Figures 4E–G). Our result showed that riskscore was significantly negatively correlated with immunescore and estimatescore (Figures 4F,G; immunescore, r = -0.340, p < 0.001; estimatescore, r = -0.140, p < 0.001).

With the Hallmark gene set as a reference set, we performed GSEA to explore the underlying oncogenetic difference between high- and low-risk patients. The result showed that the epithelial–mesenchymal translational (EMT), mTORC1 signaling, mitotic spindle, glycolysis, E2F target, and G2M checkpoint pathways were aberrantly activated in high-risk patients (Figure 5A). No remarkable difference was found in non-synonymous mutation and synonymous mutation between high- and low-risk patients (Figures 5B,C). However, we observed more TP53 mutation in high-risk patients (low-risk: 78; high-risk: 120) and FGFR3 mutation in low-risk patients (low-risk: 40; high-risk: 17) (Figure 5D).

Nowadays, immunotherapy has become an important therapeutic strategy for BLCA patients. Therefore, we assessed riskscore differences in multiple checkpoint molecules (Figure 5E). To be noted, we found statistically significant differences of critical immune checkpoint molecules PD-L1 and PD-L2 between high- and low-risk patients, indicating that the immunotherapy response rate might differ in these two groups (Figure 5E). Furthermore, we performed TIDE analysis of our combined cohort (Figure 6A). TIDE score >0 was defined as nonresponders, while they were defined as responders otherwise (Figure 6A). More responders were observed in the low-risk group than in the high-risk group (Figure 6B). TCIA results also reached the same conclusion. The result demonstrated that the low-risk group had higher TCIA scores (CTLA4_negative_PD1_positive, CTLA4_positive_PD1_negative, CTLA4_positive_PD1_positive) than the high-risk group, indicating that low-risk might have a better immunotherapy efficacy of anti-CTLA4 and anti-PD-1 therapy (Figures 6C–E). Meanwhile, we validated this conclusion in a urothelial carcinoma cohort treated with the anti-PD-L1 antibody atezolizumab, IMvigor210. The result showed that riskscore was lower in immunotherapy complete/partial response patients (Figure 6F). Moreover, the ROC curve showed that the riskscore could effectively predict the immunotherapy response of patients in IMvigor210 (Figure 6G, AUC = 0.659; SD/PD and CR/PR).

Our model could effectively indicate patients’ prognosis and immunotherapy response rate. Specific marker genes contributed to the quick genotyping in clinical features. For the model immune terms, Lymphs_PCA, Tfh cells, CD8_T_cells, cytotoxic cells, and T_cell_receptors score was obviously activated in high-risk patients, whereas core serum response up, NK cells, lymph vessels, CSR_activated, and macrophages were downregulated (Figure 7A). At the same time, 10 characteristic genes were identified between high- and low-risk patients, including KCNH4, UGT1A1, TPO, SHANK1, PITX3, MYH1, MYH13, KRT3, DEC1, and OBP2A (Figure 7B). The clinical correlation and KM survival curves of these feature genes are shown in Supplementary Figure S3. The result showed that TPO, SHANK1, PITX3, MYH13, KRT3, and DEC1 were associated with worse clinical features and prognosis, yet the KCNH4 and UGT1A1 were the contrary. Based on the logistic regression, the expression level of these 10 genes could effectively predict the risk group of bladder cancer patients (Supplementary Figure S4A, AUC = 0.856). Meanwhile, these 10 genes also have some predictive power on the immunotherapy response rate in the IMvigor210 cohort (Supplementary Figure S4B, AUC = 0.620). Based on the CMAP analysis, potential compounds targeting the riskscore were screened, for example, amantadine, carteolol, and helveticoside (Figure 7C). The mechanism of these compounds was mainly focused on ATPase inhibitor, bacterial cell wall synthesis inhibitor, cyclooxygenase inhibitor, progesterone receptor agonist, and protein synthesis inhibitor.