TCGA and GEO databases were used to gather transcriptome and clinical data. In general, we collected two groups of cohort samples of BLCA: GSE13507 and TCGA-BLCA. The exclusion criteria were as follows: 1) Not tumor tissue sample. 2) The transcriptome sequencing data or clinical information of the samples were incomplete. 3) Not common gene expression sample data. Finally, 569 samples were included. We converted the fragments per kilobase million (FPKM) values to the transcripts per million (TPM) values by using the “limma” R package for the TCGA-BLCA database. We combined TCGA and GEO expression data to get new common gene expression samples’ data for later analysis.

The “CIBERSORT” R package, the LM22 signature, and 1,000 permutations were used to evaluate infiltration levels for various immune cells in BLCA. ESTIMATE calculated the immune score and stromal score in BLCA patients. In addition, we acquired the correlation between different immune cells by using the “corrplot” R package. The hierarchical agglomerative clustering of BLCA was implemented by different ICI patterns of each sample. The number of clusters was determined by the consensus clustering algorithm. We performed the “ConsensusClusterPlus” R package and repeated it 1,000 times to ensure the stability of classification.

In order to find genes linked with ICI patterns, we classified patients into distinct groups based on ICI. DEGs among different groups were screened by means of the R package “limma.” The significant criteria of |log FC| > 1 and p (adjust) < 0.05 were used to determine DEGs.

In order to further analysis, an unsupervised clustering method for DEG analysis was applied to divide the patients into different groups. Positive and negative DEG correlations with cluster signatures were classified as ICI gene signatures A and B, and the “Boruta” algorithm was applied to reduce their dimensionality. Using the PCA, the gene signature score of patients was derived. Finally, we used a procedure analogous to the grading index gene expression to determine the ICI score. ICI score = ∑PC1A-∑PC1B.

The correlation mutation information of patients in the TCGA-BLCA cohort was obtained from the TCGA data portal (https://portal.gdc.cancer.gov/repository). In order to determine tumor mutation burden, we calculated the total number of non-synonymous mutations in BLCA. We got 20 driver genes through the R package “maftool,” which had the highest mutation frequency in BLCA patients. Finally, we evaluated whether differences in the mutation frequency of genes between two ICI score groups.

The “clusterprofiler” package was utilized for gene annotation and enrichment analysis of ICI distinctive genes. Gene ontology (GO) terms were screened by a stringent cut-off (p < 0.05). Furthermore, we used a gene set enrichment analysis (GSEA) to find pathways that were up- and down-regulated between two ICI score groups. The parameter settings were Gene sets database = “Kyoto Encyclopedia of Genes and Genomes (KEGG),” n Perm = 1,000, and p < 0.05.

The corresponding clinical data from the TCGA and GEO databases were retrieved and manually organized. We verified the predictive value of the ICI score with distinct clinical features (such as age and gender) by the R package” survival.”

All data were analyzed by R software (version 4.0.4). The proportion of 22 types of immune cells in BLCA was calculated by the “CIBERSORT” algorithm. The “ConsensusClusterPlus” R package was used to divide BLCA patients into two types. DEGs between two ICI phenotypes were filtered through the R package “limma.” The ICI score was calculated by the PCA algorithm. Further analysis, TMB was obtained in TCGA-BLCA by “TMB.pl.” The prognosis of BLCA patients was evaluated by the R package “survival.” The “clusterprofiler” package was used for gene annotation and enrichment analysis of ICI distinctive genes. The predictive value of ICI scores with different clinical characteristics (such as age and gender) was verified by the R package “survival.” p < 0.05 was considered statistically significant.

The workflow is displayed in Figure 1. First, a total of 569 common gene expression data were extracted from TCGA and GEO cohorts. Then, the CIBERSORT and ESTIMATE algorithms were applied to evaluate the levels of immune cells (filter conditions: p < 0.05) in BLCA patients. (Supplementary Tables 1, 2). The correlation coefficient heatmap shows the significant positive correlation of CD8 T cells with activated memory CD4 T cells and immune score. On the contrary, there was a significant negative correlation between CD8 T cells with macrophages M0 (Figure 3C). Based on the results of immune cell infiltration, unsupervised clustering was performed by using the “ConsesusClusterPlus” package of R software to divide BLCA patients into two distinct ICI subtypes. The consensus matrix was the crispest when K = 2 (Figures 2A–D), namely, ICI clusters A and B (Supplementary Table 3). The heatmap enables visualization of the expression levels of immune cells of distinct ICI clusters (Figure 3B). Moreover, two independent ICI subtypes showed a significant difference in the overall survival rate (p = 0.002; Figure 3A).

To further describe and understand the biological and clinical distinctions among these inherent features, we analyzed the immune cell composition of two ICI subtypes. Between two ICI subtypes, ICI cluster A had a better outcome with a median duration of roughly five years. Meanwhile, it was marked by increased infiltration of CD8 T cells, activated memory CD4 T cells, resting mast cells, etc. In addition, the ICI cluster A has a higher immune score than the ICI cluster B. On the contrary, the ICI cluster B confirmed a poor prognosis (median survival duration of roughly three years) and showed a large rise in the number of macrophages M0 (Figure 3D).

In order to reveal the underlying biological properties of distinct immunophenotypes, we used the R package “limma” to carry out differential analysis to identify the transcriptome differences between two subtypes. Unsupervised clustering was implemented using the “limma” package of R software to obtain the differentially expressed genes (DEGs) (Supplementary Table 4). We classified BLCA patients into gene clusters A–B by DEGs (Figures 4A–D; Supplementary Table 5). The positive correlation of DEG values with the clusters signature was coded as ICI gene signature A, and the rest of the DEGs were coded as ICI gene signature B. At the same time, to reduce noise or redundant genes, we employed the “Boruta” method to reduce the dimension of gene signatures A and B. The transcriptome properties of DEGs are shown in a heatmap created with the R package “pheatmap (Figure 5A).” The R package “clusterProfiler” was applied to execute GO enrichment analysis on the signature genes. The significantly enriched biological processes are summarized in Figures 5B–E, and a detailed description is provided in Supplementary Table 6.

Next, we evaluated the prognosis of gene clusters A–B combined with survival information. It showed that two independent gene clusters had a significant difference in overall survival (p < 0.001; Figure 5F). The gene cluster A was characterized by a better outcome (median survival duration of roughly 5.5 years), whereas the gene cluster B had unfavorable outcomes (median survival duration of roughly two years). As displayed in Figure 5G, the gene cluster A showed an obvious increase in the infiltration of some immune cells, such as CD8 T cells and naive B cells. Meanwhile, the gene cluster B performed a higher macrophages M0 infiltration. Finally, some differentially expressed target genes were analyzed in two gene clusters by the “limma” package.

In order to obtain the quantitative index of the ICI landscape of BLCA patients, we computed two aggregate scores by the PCA algorithm: the ICI score A from ICI signature gene A and the ICI score B from ICI signature gene B. We obtained the sum of individual scores using ICI scores A and B of each sample in the study. Finally, we obtained the ICI score, which is a predictive signature score. The TCGA-BLCA and GSE13507 patients were separated into high and low ICI score groups using the “survival” package (Supplementary Table 7). The alluvial diagram described the correlation among the gene clusters, the ICI score, and survival outcomes (Figure 6A). CD274, CTLA4, HAVCR2, LAG3, and PDCD1 were chosen as immune-checkpoint–relevant signatures, and CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF were selected as immune-activity–related signatures, to investigate the immunological activation and tolerant state of the TCGA-BLAC and GSE13507 cohorts. With the exception of TBX2, the ICI score was shown to have a substantial negative correlation with the expression quantity of immune-checkpoint–relevant and immune-activity–relevant genes (Figure 6B). Moreover, GSEA analysis results revealed that fatty acid metabolism and PPAR signaling pathways were considerably enriched in the high ICI score group, whereas proteasome and NOD-like receptor signaling pathways were substantially enriched in the low score one (Figure 6C). Detailed enriched information and description is provided in Supplementary Table 8. Next, we evaluated the impact of the ICI score on the prognosis of patients. It showed that two independent ICI score groups had remarkable difference in the overall survival rate (p < 0.001; Figure 6D). The high ICI score group showed a good prognosis (median survival duration of roughly 5.3 years), whereas the low one had an unfavorable outcome (median survival duration of roughly 1.2 years). Lastly, the “ggplot2” package was applied to evaluate the relation between ICI score and survival status. We found that two independent ICI score groups had a significant difference in survival status. The majority of BLCA in the high ICI score group were alive; on the contrary, the majority of BLCA in the low one was dead (p = 0.0059; Figures 6E,F).

Since BLCA was reported to have a high degree of somatic changes, subsequently, we determined the distribution of somatic mutations and combined it with the ICI score to evaluate the prognosis of patients. First, the total mutation burden and mutation distribution of TCGA-BCLA were obtained by analyzing mutation annotation files. Meanwhile, we divided the patients into high- and low-TMB groups. As demonstrated in Figure 7A, we discovered that the high-TMB group was related to a better outcome than the low one (p < 0.001). Considering the contraindication value of TMB and ICI score for prognosis, we subsequently studied the synergistic effect of ICI score in the prognostic classification of BLCA. The results reveal that there was a substantial difference in survival between the high- and low-TMB groups depending on ICI score subtypes. Among them, the high TMB combined with a high ICI score had the best prognosis in BLCA (p < 0.001; Figure 7B). In conclusion, the ICI score could be utilized as a possible predictor different from TMB, which could effectively predict the response to immunotherapy. In addition, we screened out 20 driver genes with the highest mutation frequency for further analysis. We analyzed the distribution of driver genes in high and low ICI score groups. The result showed that the alteration frequency of TP53, KMT2D, PIK3CA, KMT2C, and FLG was considerably different between the two ICI score groups (Figures 7C,D). Moreover, we discovered that the TP53 mutation frequency was higher in the low ICI score group. The result proved once again that the group with the low ICI had a poor prognosis. These results may provide new ideas for the study of the mechanism of ICI in tumors (Şenbabaoğlu et al., 2016; Chen and Mellman, 2017).

To be able to clarify the role of ICI score in BLCA, the relationship between ICI score and clinical characteristics was researched. The stratified survival analysis was used to observe whether ICI scores could be applied to different clinicopathological features. Next, we analyzed patients’ age and gender. Results showed that the ICI score could effectively forecast OS in all groups from the age and gender clinical features (Figures 8A–D).