The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) is a landmark cancer genomics program that molecularly describes over 20,000 primary cancer, and matches normal samples spanning 33 cancer types. This joint effort between National Cancer Institute (NCI) and the National Human Genome Research Institute began in 2006, and has produced over 2.5 petabytes of genomic, epigenomic, transcriptomic, and proteomic data. The data, which has already led to improvements in our ability to diagnose, treat, and prevent cancer, will remain publicly available for anyone in the research community to use. We downloaded FPKM standardized RNA-seq data, clinical information and tumor mutation burden (TMB) information from the TCGA-BLCA cohort in TCGA database.

ImmPort (https://www.immport.org/) is funded by the National Institute of Health (NIH) and National Institute of Allergy and Infectious Diseases (NIAID) in support of the NIH mission to share data with the public. We clicked the “Resources” button on the Immport database homepage, then clicked the “Gene Lists” button on the “Resources” page, and finally clicked the “Gene Summary” to download immune-related genes.

Gene Expression Omnibus (GEO) is a public functional genomics data repository supporting MIAME-compliant data submissions. Array- and sequence-based data are accepted. Tools are provided to help users query and download experiments and curated gene expression profiles. We downloaded two data sets (GSE31684 and GSE39281) recording bladder cancer transcriptome genes (RNA-seq) and clinical information in the GEO database. After processing the data with Perl, we obtained two gene expression matrices. Then, we used the “sva” package in the R language (version 4.0.5) to merges the two expression matrices and eliminate batch effects.

(A) TMB analysis: we used BLCA mutation data in the TCGA database and Perl language to calculate the number of base mutations in each BLCA sample. (B) Single-sample gene set enrichment analysis (ssGSEA) and hierarchical cluster analysis: we used R packages (GSVA, GSEABase and limma) to perform ssGSEA to calculate the immune score of each BLCA sample according to 29 immune gene sets composed of different types of immune cells with different functions, pathways and checkpoints (Alhamdoosh et al., 2017). Firstly, the rank of gene expression values in a given BLCA sample was normalized, and then the enrichment score (ES) was calculated using the empirical cumulative distribution function. Each ssGSEA score XI was converted to XI′ by bias normalization to obtain the scores of different immune cells and immune-related functions in each sample. Then, we used the hierarchical clustering method of Euclidean distance and Ward linkage to do the immune stratification of BLCA patients. Meanwhile, we also made use of the T-distribution stochastic neighbor embedding (tSNE) algorithm to determine the immune stratification of BLCA patients through RtSEN package (Gardner et al., 2021). (C) Evaluation of tumor immune microenvironment: based on ESTIMATE algorithm, BLCA transcriptome data was utilized to predict stromal cell score, immune cell score and tumor purity, and then the content of these two types of cells was predicted, from which StromalScore, ImmuneScore and EstimateScore were determined (Yoshihara et al., 2013). (D) Tumor-infiltrating immune cells analysis: CIBERSORT, an R tool, was used for the deconvolution of the expression matrix of human immune cell subtypes according to linear support vector regression. This method is based on a known reference set and provides a set of gene expression characteristics of 22 immune cell subtypes. Therefore, we used the CIBERSORT method to do the calculation for the abundance of infiltrating immune cells in BLCA samples (Newman et al., 2015). (E) Immune differential genes determining the immune stratification: the limma package was utilized to select differentially expressed genes (DEGs) among people with different immune stratification (| log2 fold change | > 1.50 and FDR <0.05), and then we obtained immune-related genes from ImmPort (Bhattacharya et al., 2014). DEIGs were obtained through the intersection of immune genes and DEGs. (E) Prognostic markers: the survival package was utilized to do the univariate Cox regression analysis (p < 0.05) to identify the markers of significant prognosis-related immunity genes (PIMGs).

Gene set enrichment analysis was performed via the GSEA software (version 4.1.0) to analyze TCGA-BLCA transcriptomes for the identification of the key signaling pathways involved in DEGs.

The major Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways involved in the up-regulation of the Immunity_H and Immunity_L subgroups (p < 0.05, FDR <0.01) were selected. R (version 4.0.5) was used to perform further analysis, and visualize the results. Then, we obtained transcription factors associated with the occurrence and development of bladder cancer from the CISTROME project, extracted differentially expressed transcription factors (DETFs) from the total DEGs, and used Pearson correlation coefficient analysis to construct the regulatory network of PIMGs and DETFs (R > 0.3 and FDR <0.01) (Mei et al., 2017). Finally, the protein-protein interaction (PPI) network analysis was performed using STRING (String-db.org/).

We used the LASSO Cox regression model in R package (Dalal et al., 2012) “glmnet” to find genes significantly associated with the prognosis to construct the prognostic model of BLCA (PMB). The risk score was calculated as the following formula: r i s k S c o r e = ∑ i = 1 9 β i ∗ L P I M G i , where LPIMG_i represented the i-th LPIMG (Lasso-prognosis-related immunity genes), and β _i represented the expression coefficient of LPIMG_i obtained from Lasso regression analysis. All cases were classified into a low-risk group and a high-risk group based on the median risk score, and we performed the Kaplan-Meier survival analysis to compare the survival status between the high-risk group and the low-risk group. In order to verify the predictive power of PMB, the receiver operator characteristic (ROC) curve was drawn to calculate the area under the curve (AUC) of 1-, 3-, and 5-year survival. We conducted Kaplan-meier, logarithmic rank, ROC curve and calibration analysis using “timeROC,” “rms,” “survival,” and “survminer” software packages in R language. Based on the risk score calculated by PMB, Pearson correlation coefficient, Spearman correlation coefficient and corrplot package were used to evaluate the correlation between the risk score and overall survival, immune cell infiltration, immune checkpoint molecules and TMB. p < 0.05 of the critical value for the significant correlation was set. Eventually, univariate and multivariate Cox regression analysis of the risk scores of the constructed PMB and patients’ clinical characteristics (age, sex, stage) was performed to verify the accuracy of the independence of PMB-based risk characteristics. Based on the above factors, we created a nomogram using the R packages of “rms”, “nomogramEX” and “regplot.” Finally, the ROC and calibration chart were drew to determine the suitability of our established nomogram for potential clinical applications.

In order to fully evaluate the immunological characteristics of BLCA, we used the ssGSEA to analyze 414 tumor samples from the TCGA-BLCA cohort. According to the ssGSEA scores and hierarchical clustering method, BLCA cases were divided into two clusters. The average score of the immune microenvironment of the first cluster was 0.62, and the average score of the immune microenvironment of the second cluster was 0.49. Thus, the first cluster was set as the Immunity_H (high) group, and the second cluster as the Immunity_L (low) group (Figures 1A,B). The tSNE was further used to analyze the immune levels for different BLCA patients and the same classification was obtained (Supplementary Figure S1A). The results of ESTIMATE analysis indicated that EstimateScore (419.27 ± 1649.47), ImmuneScore (750.39 ± 886.17), and StromalScore (-331.12 ± 910.28) in the Immunity_H group were significantly higher than those which were (−2283.37 ± 727.70), (−620.62 ± 352.59), and (−1662.75 ± 487.21), respectively, in the Immunity_L group (Wilcox test, p < 0.001) (Figure 1C). CIBERSORT was used to detect the degree of immune cell infiltration in the tumor, which found that the differences between the Immunity_H group and the Immunity_L group in T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, NK cells resting, NK cells activated, Macrophages M1 and Mast cells activated were significant (Supplementary Figure S1B). The expression of human leukocyte antigen (HLA) genes in the both groups was examined, which suggested that most of HLA genes significantly increased in the Immunity_H group and significantly decreased in the Immunity_L group (Wilcox test, p < 0.05) (Figure 1D). Based on our results, we believed that immune response might play important roles in the development of BLCA.

Identification of immune-related genes associated with bladder cancer and their correlation with prognosis.

We further studied the expression of differential genes of immune stratification in BLCA patients. The FDR values and log2 fold change multiples of the immune differential genes in the Immunity_H group and the Immunity_L group were showed in Figure 2A. After primarily screening, we totally identify 994 DEGs, of which 812 genes were up-regulated and 82 genes were down-regulated (Figure 2B). Subsequently, 308 DEGs were selected as DEIG using the ImmPort database (Figure 2C). Univariate Cox regression analysis indicated that 13 PIMGs had significant association with the survival of BLCA patients in DEIGs (p < 0.01), of which seven genes, including NPR2, TGFB3, PDGFRB, PDGFRA, VIM, RBP1, RBP1 and TNC, increased the risk of prognosis, while the rest, including CD3D, CIITA, GNLY, LCK, PDCD1 and ZAP70, were conducive to survival (Figure 2D).

LASSO Cox regression analysis was performed on 13 selected PIMGs (Figures 3A,B). Finally, 9 LPIMGs were identified and their risk-correlation coefficients were calculated to determine the prognosis of BLCA patients. The risk score was calculated as follows: riskScore = NRP2*0.0101119 + CD3D*−0.1990949 + GNLY*−0.1241769 + LCK*−0.0519549 + VIM*0.1464182 + RBP1*0.1038418 + PDGFRA*0.1589969 + ZAP70*−0.12644895 + TNC*0.0693184. Data from TCGA was selected as the training group, the risk score of each BLCA case in this group was calculated, and all cases were classified into the high-risk group (203 patients) and the low-risk group (204 patients) based on the median risk score of 0.4886 (Supplementry Data S1; Supplementary Figure S1C). The correlation analysis indicated that the risk score had significant negative correlation with the survival time of BLCA patients which gradually decreased with the increase of the risk score (Supplementary Figures S1D,E). The Kaplan-meier curve showed that the difference in overall survival (OS) between the high-risk group and the low-risk group was significant, and patients in the low-risk group had a longer overall survival time than those in the high-risk group (Log-rank test, p < 0.0001) (Figure 3C). In order to evaluate the predictive power and accuracy of PMB-based risk characteristics, the ROC curves of the training group were drawn, and the AUC values of 1-, 3- and 5-year survival were 0.688, 0.719, and 0.706, respectively (Supplementary Figure S1F). The accuracy of the prognostic model was verified by the calibration chart, which suggested that the predicted value of the prognostic model was in good consistence with the actual value (Figure 3H). Besides, GSE31684 and GSE39281 were used as the external validation group, and we combined their data (GSECD) using R “sva” package to further confirm the accuracy and feasibility of the prognostic model, and the number of deaths in the high-risk group increased significantly (Supplementary Data S2). Then, the Pierce correlation analysis and Kaplan-Meier curves suggested that the constructed PMB-based risk characteristics still had good predictive power in the external validation group (Figure 3D).

In order to investigate the correlation between immunotherapy and bladder cancer, 14 immune checkpoint inhibitors inlcuding BTLA, GITR, TNFRSF14, IDO, LAG-3, PD-1, PD-L1, PD-L2, CD28, CD40, CD80, CD137, CD27, and Ctla-4 were selected for analysis. It was found that the risk score had negative correlation with the BTLA, CD27, CD40, CD80, and TNFRSF14 expression, which had significant differences in different risk groups (Supplementary Figure S2), indicating that tumor immunosuppression might lead to an increased risk score of patients. TYK2 and ACE2 were also differentially expressed in different risk groups, and with the increase of the risk score, their expression decreased (Figure 4A). In the TCGA-BLCA cohort, the TMB of patients in the high-risk group was significantly lower than that in the low-risk group (p = 0.009) (Figure 4B). In order to find the potential correlation between TMB and the prognosis of patients, according to the TMB cutoff value of 4.632, we divided the patients into the high TMB group, and the low TMB group (Supplementary Data S3). We found that the survival time of patients in the high TMB group was significantly lower than that of the low TMB group (p < 0.001) (Figure 4C). In order to evaluate the outcomes of patients more comprehensively, we investigated whether the combination of the risk score and TMB could be a more accurate prognostic marker. We integrated PMB-based risk characteristics with TMB, stratified all samples into the H-TMB/high risk, H-TMB/low risk, L-TMB/high risk, and L-TMB/low risk groups. Figure 4D suggested that differences between groups were significant (log-rank test, p < 0.0001), and in the H-TMB/low risk group, patients had the longest overall survival. The above results together suggested that the risk score had positive correlation with the degree of malignant tumor.

Due to the significant correlation between the risk score and the degree of malignant tumor, univariate and multivariate Cox regression analysis for age, sex, and stage as covariates was conducted to test the potential possibility of the risk score as an independent prognostic factor for BLCA patients, of which the results showed that the PMB based risk characteristics had a p value less than 0.001, confirming that the PMB based risk characteristics could be used to predict the prognosis of BLCA patients (Table 1). Combined with the above factors, we constructed a nomogram (Figure 5A) to expand the clinical application and usability of PMB. The total score of each patient was obtained by calculating and summing the score for each prognostic parameter. The higher the total score was, the worse the patient’s clinical outcome was. The ROC curve showed that the nomogram had a good predictive ability for the survival rate, with a high accuracy, and the AUC values of 1-, 3-, and 5-year survival were 0.744, 0.770, and 0.782, respectively (Figure 5B). In addition, the calibration chart indicated that the nomogram performed similarly with the ideal model (Figure 5C).

GSEA revealed that immune-associated pathways in the Immunity_H group were highly active, including the signaling pathway of T cell receptor, the pathway of antigen processing and presentation, cytokine involved immune response, and hematopoietic cell lineage. Additionally, various pathways of immune-associated disease were identified in the Immunity_H roup, including asthma, primary immune deficiency, graft-versus-host disease, allograft rejection, thyroid disease related to autoimmune, and immunity to leishmania infection (Figure 6A). In order to clarify the role of the multi-dimensional regulatory network of immune molecules in the occurrence and development of bladder cancer, we firstly explored the upstream mechanism of PIMG. By combining differential expression analysis with data from the CISTROME database, we identified transcription factors significantly associated with the BLCA prognosis. For the Immunity_H subtype, a total of 7 up-regulated transcription factors were identified. Figure 6B showed the regulatory network of BLCA TF-PIMGs. PPI analysis was further conducted and we confirmed the significant correlation between BLCA TF and PIMG (Figure 6C).