The Pattern of ICI in the TME of BLCA
The workflow is displayed in Fig. 1. Firstly, a total 569 common gene expression data were extracted from TCGA and GEO cohorts. Then, the CIBERSORT and ESTIMATE algorithms were applied to assess the level of immune cells (Filter conditions: p<0.05) in BLCA patients. (Supplementary Table 1 & 2). The correlation coefficient heatmap shows that the significant positive correlation of CD8 T cells with activated memory CD4 T cells and immune score. On the contrary, the significant negative correlation of CD8 T cells with Macrophages M0 (Fig. 3C). Based on results of immune cell infiltration, BLCA patients were divided into two different ICI subtypes by R package "conesusclusterplus" unsupervised clustering. The consensus matrix was the crispest when K = 2 (Fig. 2A-2D), namely ICI cluster A and B (Supplementary Table 3). The heatmap enables visualization of the amounts of expression of immune cells of distinct ICI clusters (Fig. 3B). Moreover, two independent ICI subtypes showed significant difference in the overall survival rate (p = 0.002; Fig. 3A).
To better explain and comprehend 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 5 years. Meanwhile, it was marked by increased infiltration CD8 T cells, activated memory CD4 T cells and resting Mast cells, ect. In addition, the ICI cluster A has higher immune score than the ICI cluster B. On the contrary, the ICI cluster B confirmed a poor prognosis (median survival duration roughly 3 years) and performed a large rise in the amount of Macrophages M0 (Fig. 3D).
Identified the subtypes of immune related gene
To understand the underlying biological properties of distinct immunophenotypes, we used R package “limma” to carry out differential analysis to identify the transcriptome differences between two subtypes. To determine the DEGs, unsupervised clustering was implemented by the "limma" package (Supplementary Table 4). We classified BLCA patients into gene clusters A–B by DEGs (Fig. 4A-4D; Supplementary Table 5). The positive correlation of DEGs values with the clusters signature, were coded as ICI gene signature A, while the remaining DEGs were coded as ICI gene signature B. At the same time, to remove interference or duplicated genes, we employed the "Boruta" method to minimize the dimensionality of gene signatures A and B. The transcriptome properties of DEGs are shown in a heat map created with the R package "pheatmap." The R package "clusterProfiler" was applied to execute GO enrichment analysis on the signature genes. Fig. 5B-5E illustrates the biological processes that have been significantly enriched, and Supplementary Table 6 contains a thorough explanation.
Next, we evaluated the prognosis of gene clusters A–B combined with survival information. It showed that two independent gene clusters had significant difference in overall survival (p < 0.001; Fig. 5F). The gene cluster A was characterized with a better outcome (median survival duration roughly 5.5 years), whereas the gene cluster B was linked to a poor outcome. (median survival duration roughly 2 years). As displayed in Fig. 5G, the gene cluster A showed obvious increase in the infiltration of some immune cells, like 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.
Generation of ICI Score
In develop a quantification of the ICI environment in BLCA patients, the PCA algorithm was applied to generate two overall scores: 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 "survival" package (Supplementary Table 7). The alluvial diagram described the correlation among the gene clusters, the ICI score and survival outcomes (Fig. 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. (Fig. 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 one (Fig. 6C). A detailed enrich information description was provided in Supplementary Table 8.
Next, we investigated at how the ICI score affected a patient's prognosis. It showed that two independent ICI score groups had remarkably difference in the overall survival rate (p < 0.001; Fig. 6D). The high ICI score group showed a good prognosis (median survival duration roughly 5.3 years), whereas the low one had the unfavorable outcome (median survival duration roughly 1.2 years).
Lastly, "ggplot2" package was applied to evaluate the relation between ICI score and survival status. We found that two independent ICI score groups had significant difference in survival status. The majority of BLCA in high ICI score group were alive, on the contrary, the majority of BLCA in low one were dead (p = 0.0059; Fig. 6E & F).
TMB & ICI score were applied to assess the prognosis of TCGA-BCLA cohort patients
Since BLCA was reported to have high degree of somatic changes, subsequently, we determined the distribution of somatic mutations and combined with ICI score to evaluate the prognosis of patients. Firstly, the total mutation burden and mutation distribution of TCGA-BCLA were obtained by analyzing mutation annotation files. Meanwhile, we sorted the patients into high and low TMB groups. As demonstrated in Fig. 7A, we discovered that high TMB group was related to 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 high ICI score had the best outcomes in BLCA (p <0.001; Fig. 7B). In conclusion, The ICI score could be utilized as a possible predictor irrespective of TMB, which could effectively predict the response of immunotherapy.
In addition, twenty driver genes with the greatest mutation frequency were chosen for research development. We evaluated at the allocation of driver genes in two ICI score groups. The result showed that the alteration frequency of TP53, KMT2D, PIK3CA, KMT2C and FLG was considerably different between two ICI score groups (Fig. 7C-7D). Moreover, we discovered that the TP53 mutation frequency was higher in 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 tumor.
Analysis of Clinical traits in two ICI score groups
To be able to clarify the function 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 score could be applied to different clinicopathological features. Next, we analyzed patients' age and gender. Results showed that ICI score could effectively forecast OS in all groups from the age and gender clinical characteristics (Fig. 8A-8D).