Tumor immune microenvironment-based classifications of bladder cancer for enhancing the response rate of immunotherapy.

Immunotherapy is a potential way to save the lives of patients with bladder cancer, but it only benefits approximately 20% of them. A total of 4,028 bladder cancer patients were collected for this study. Unsupervised non-negative matrix factorization and the nearest template prediction algorithms were employed for the classification. We identified the immune and non-immune classes from The Cancer Genome Atlas Bladder Urothelial Carcinoma (TCGA-BLCA) training cohort. The 150 most differentially expressed genes between these two classes were extracted, and the classification reappeared in 20 validation cohorts. For the activated and exhausted subgroups, a stromal activation signature was assessed by the NTP algorithm. Patients in the immune class showed highly enriched signatures of immunocytes, while the exhausted subgroup also exhibited activated transforming growth factor (TGF)-β1, and cancer-associated extracellular matrix signatures. Patients in the immune-activated subgroup showed a lower genetic alteration and better overall survival. Anti-PD-1/PD-L1 immunotherapy was more beneficial for the immune-activated subgroup, while immune checkpoint blockade therapy plus a TGF-β inhibitor or an EP300 inhibitor might achieve greater efficacy for patients in the immune-exhausted subgroup. Novel immune molecular classifier was identified for the innovative immunotherapy of patients with bladder cancer.

Introduction

Bladder cancer is the 10th most frequent tumor globally and exhibits a high rate of recurrence. The major challenge of clinical care of bladder cancer is the short-term recurrence of non-muscle-invasive bladder cancer (NMIBC), as well as the shortened overall survival of muscle-invasive bladder cancer (MIBC) patients, especially those with distant metastases, the 5-year survival rate of whom less than 10%., In the tumor mass, the normal cells, blood vessels, and cytokines that surround and support the vitality of tumor cells compose the tumor microenvironment (TME). Crosstalk exists between the tumor and the TME. Tumor cells can alter the TME, and the TME can promote the growth and spread of tumors.

From the molecular side, bladder cancer is composed of a multitude of heterogenetic characteristics, including gene mutations, gene copy number alterations, and neoantigens, as well as the infiltration of immunocytes. Several teams have established molecular classifications among bladder cancer. Mo et al. generated an 18-gene tumor signature in MIBC patients that can reflect the urothelial differentiation and predict clinical outcomes; basal and differentiated groups are the two groups with the highest and lowest risk scores, respectively. Damrauer et al. developed BASE47, a 47 gene-based classifier, for the separate of luminal-like or basal-like subtypes of MIBC tumors. Robertson et al. further generated a consensus hierarchical clustering of luminal-papillary, luminal-infiltrated, luminal, basal/squamous, and neuronal subtypes. However, most of the molecular classifiers only focused on the clinical outcomes, not the tumor immune microenvironment. Therefore, our goals were to provide comprehensive insight into the immune response of bladder cancer patients with diverse inner molecular features and generate the classifier to screen the patients suited for immunotherapy.

The non-negative matrix factorization (NMF) algorithm is a multiplicative updates algorithm; it can decompose a non-negative matrix V into two non-negative matrices, W and H. Similar to principal component analysis (PCA) or independent component analysis (ICA), the NMF algorithm can also use a limited number of components to reflect the original observed data, which might contain huge volumes. NMF has been applied to reveal biomarkers, classify tumor subtypes, and predict the prognosis of tumors in recent works.9, 10, 11, 12

We enrolled 4,028 patients with bladder cancer from independent cohorts. The NMF and nearest template prediction (NTP) algorithms were applied to distinguish patients with different immunophenotypes in The Cancer Genome Atlas Bladder Urothelial Carcinoma (TCGA-BLCA) training cohort and reappeared in the validation cohorts. The novel definition of these immunophenotypes could provide illumination for the immunotherapy of patients with bladder cancer.

Results

Identification of the immune module and derivation of the immune class of bladder cancer

4,028 bladder cancer patients were involved, along with the matched overall survival data, clinicopathological information, and gene expression profiles (Figure 1). We performed virtual microdissection using the NMF algorithm in the TCGA-BLCA training cohort. To obtain the robust immune module, we preset the respective module numbers as five to 10. When the total module number was nine, the first module strongly enriched patients with high immune enrichment scores, which were defined as the immune module (Figure 2A). The top 150 weighted genes in the immune module were defined as exemplar genes that reflected the characteristics of the immune module (Table S1). According to the ontological analysis of biological processes, these genes are associated with the activation of immunocytes, T helper 1 (Th1)/Th2 cell differentiation, T cell receptor signaling, and B cell receptor signaling (all p < 0.05; Table S2). Subsequently, we redefined the 408 bladder patients into immune-enriched or non-immune-enriched groups via the consensus clustering analysis of the 150 exemplar genes (Figure 2B). Furthermore, multidimensional scaling (MDS) random forest (RF) was further applied to define a more precise classification for the immune and non-immune classes (Figure 2C). In Figure 2D, the distributions of the 408 bladder cancer patients among the NMF modules, immune module weight, exemplar gene clustering, final immune classes, and immune enrichment score are shown.

Figure 1

Flow chart of the analyses performed in this study

NMF, non-negative matrix factorization; TCGA-BLCA, The Cancer Genome Atlas Bladder Urothelial Carcinoma; TILs, tumor-infiltrating lymphocytes; CNA, copy number alteration; TMB, tumor mutation burden.

Figure 2

Recognition of the immune classes by the non-negative matrix factorization (NMF) algorithm

(A) Nine modules were generated from the NMF algorithm, and the module gathered patients with high immune enrichment score were recognized as the immune module. (B) Heatmap showing the top 150 exemplar genes expression among immune-enriched and non-immune-enriched clusters, divided by consensus clustering. (C) The multidimensional scaling random forest further modified the clusters to immune and non-immune classes. (D) The distributions of patients in different NMF modules, immune module weight, exemplar gene clustering, final immune classes, and immune enrichment score.

Several immune-associated signatures (Table S3) were collected to help confirm the classification of the immune or non-immune classes, and the enrichment score of each signature for each patient was determined by single-sample gene set enrichment analysis (ssGSEA). We observed the increased enrichment of immunocytes in the immune class as compare with the non-immune class, including T cells (as reflected by the signatures of 13 T cell signature, T cells, CD8+ T cells, and T. NK. Metagene), B cells (as reflected by the signatures of B cell clusters and B.P. metagene), macrophages, tertiary lymphoid structure (TLS), cytolytic activity score (CYT), and interferon (IFN) signatures (all p < 0.05; Figure 3A). We also analyzed the activated KEGG signaling pathways by GSEA, revealing that immune cell pathways (including T cell-, B cell-, natural killer cell-, and leukocyte-associated pathways), immune response pathways (including chemokine signaling pathways, antigen processing presentation, cell adhesion molecules, and complement coagulation cascades), and proinflammatory pathways (including FC-Epsilon-RI-, NOD-like receptor-, and FC gamma R-mediated phagocytosis pathways) were all activated in the immune class (Figure S1). From the results of Figures 2, 3A (top panel), and S1 and Tables S1, S2, and S3, we microdissected the immune and non-immune classes in the TCGA-BLCA cohort, and activated immune-associated signatures and signaling pathways were observed in the immune class.

Figure 3

The diverse immune characteristics and heterogeneity of genetic phenotypes of non-immune class, immune-activated subgroup, and immune-exhausted subgroup

(A) Division and characterization of three immunophenotypes. CYT, cytolytic activity score; TITR, tumor-infiltrating Tregs; MDSC, myeloid-derived suppressor cell; TLS, tertiary lymphoid structure; C-ECM, cancer-associated extracellular matrix. (B) Difference of tumor-infiltrating lymphocyte abundance. (C) Difference in the PD-L1 mRNA expression level. (D) Difference in gene copy number alterations, including amplification and deletion, among arm levels and focal levels. (E) Difference the tumor mutation burden. (F) Difference in tumor neoantigens. (G) Specific mutant genes in the immune-activated subgroup. (H) Specific mutant genes in the immune-exhausted subgroup. WT, wild-type; IM-Act, immune-activated subgroup; IM-Exh, immune-exhausted subgroup.

Tumor immune microenvironment immunophenotypes are distinguished by the activation of stromal cells

Fibroblasts, mesenchymal stromal cells (MSCs), and the extracellular matrix (ECM) are the key components of the tumor stroma and support and connect tumor cells. Especially during the late stages of tumors, genetic and epigenetic alterations of tumor cells are driven by activated stromal components. MSCs act as inherent regulators of tumors and can secrete inhibiting soluble factors and alter cell surface markers to suppress the immune microenvironment. MSCs can regulate the expression of PD-L1 and impact the proliferation and induction of T regulatory cells (Tregs)., MSCs suppress immune processes by decreasing the expression of proinflammatory factors, including IFN-γ, tumor necrosis factor (TNF)-α, and interleukin (IL)-1β, or by promoting the expression of type 2 factors IL-10 and IL13.17, 18, 19 For this reason, the previously defined stromal activated signature was used for the further separation of patients with high immunocyte infiltration to immune-activated and immune-exhausted immunophenotypes, which could reflect the immune response status. A total of 11.0% (45/408) of bladder cancer patients in the TCGA-BLCA cohort were recognized as the immune-activated subgroup, while the remaining 110 patients (27.0%, 110/408) belonged to the immune-exhausted subgroup, with the activated stromal phenotype (Table 2; Figure 3A). Cancer-associated ECM (C-EC) regulated by the fibroblasts can recruit immunosuppressive cells, transforming growth factor (TGF)-β is an accepted immunosuppressor in the immune microenvironment, and Tregs and MDSCs can reflect the immune-exhausted status in the TME.20, 21, 22, 23 These signatures were evaluated by ssGSEA, and we revealed that the tumor-infiltrating Tregs (TITRs), WNT/TGF-β, TGF-β1-activated, and C-ECM signatures were higher in the immune-exhausted subgroup than in the immune-activated subgroup (all p < 0.05; Figure 3A; Figure S2). TIM-3 and LAG3 are reported to be associated with immune exhaustion status,, and we also found similar results in the immune-exhausted subgroups; increased TIM-3 (p = 0.008) and LAG3 (p = 0.218) were observed in the immune-exhausted subgroup (Figure S2). Based on the results from Figure 3A (bottom panel) and Figure S2, we separated the immune class into the immune-activated and immune-exhausted subgroups. The stromal enrichment score, TITR, myeloid-derived suppressor cell (MDSC), and WNT/TGF-β signatures increased in the immune-exhausted subgroup, as validated by the immune-exhausted markers TIM-3 and LAG3.

Table 2

Summary of the clinicopathological parameters of the TCGA-BLCA, GSE32894, and E-MTAB-1803 cohorts

	TCGA-BLCA (n = 408)	GSE32894 (n = 308)	E-MTAB-1803 (n = 70)
Age
≤70	230	143	42
>70	178	165	28
Sex
Male	301	228	59
Female	107	80	11
Stagea
Ta	–	116	–
T1	11	97	–
T2	191	85	24
T3	157	7	28
T4	43	1	18
Gradeb
G1/low	21	48	–
G2	–	103	4
G3/high	384	154	66
Smokingc
No	109	–	–
Yes	286	–	–
Status
Alive	229	199	28
Dead	177	25	42

Six samples lacked T stage data in the TCGA database, and two samples lacked data in GEO: GSE32894.

Three samples lacked grade data in the TCGA database, and three samples lacked date in GEO: GSE32894.

Two samples lacked alive status data in the TCGA database, and 84 samples lacked data in GEO: GSE32894.

Heterogeneity of genetic phenotypes among the immune classes

To confirm the infiltration of immunocytes among the immune and non-immune classes, as distinguished by the mRNA expression profiles of the exemplar genes, we compared the tumor-infiltrating lymphocyte (TIL) abundance of 408 bladder cancer patients, which was anteriorly estimated by hematoxylin and eosin (H&E) staining, and observed that the TIL abundance was higher in the immune class than in the non-immune class (p < 0.001; Figure 3B), consistent with the definition of these two groups. Furthermore, we also observed a higher PD-L1 expression in the immune class than in the non-immune class (p < 0.001; Figure 3C). Gene copy number alteration (CNA), tumor mutant burden (TMB), and neoantigens were reported to exhibit crosstalk with tumor immune activation. Patients in the non-immune class showed increased levels of deletion at both the arm and focal levels (pArm-del < 0.001, pFocal-del = 0.007) but not CNA amplification (pArm-Amp = 0.733, and pFocal-Amp = 0.065) (Figure 3D), which reflected the positive association of immune infiltration and gene CNA deletion. With the online tool TIMER, we twice confirmed that the association between immune infiltration and gene CNA deletion; deep deletion and arm-level deletion of PD-1, PD-L1 and CTLA-4; and the three major immune checkpoints were linked with decreased immunocyte infiltration, especially for CD4+ T cells, neutrophils, and dendritic cells (Figure S3).

The TMB in the immune class was higher than that in the non-immune class (p = 0.01; Figure 3E), while the neoantigen level exhibited no difference (p = 0.109; Figure 3F). We further compared the specific gene mutations in the immune subgroups (Figure S4A). Mutations of TP53 (53.5% versus 43.1%, p = 0.051), TTN (52.9% versus 39.5%, p = 0.011), PIK3CA (28.0% versus 17.0%, p = 0.007), and RB1 (26.0% versus 13.0%, p < 0.001) appeared more frequently in the immune class than in the non-immune class (Figure S4B). ERBB2 (p = 0.035), KMT2A (p = 0.013), PKHD1 (p = 0.007), and MDN1 (p = 0.015) were the specific mutations noted in the immune-activated subgroup (Figure 3G), and patients with EP300 (p = 0.020), HMCN1 (p = 0.014), AKAP9 (p = 0.003), and MACF1 (p = 0.016) mutations were more highly enriched in the immune-exhausted subgroups (Figure 3H). The mutation of EP300 lead to the increased expression of EP300 (Figure S4C). Based on the results from Figures 3B–3H, S3, and S4, we can conclude that the immune class exhibits lower copy number deletion, higher TIL abundance, higher TMB, and higher PD-L1 level. The specific mutant genes in the immunophenotypes are diverse.

Reappearance of the three immunophenotypes in external cohorts

External cohorts with mRNA expression profiles were collected to recapitulate the three immunophenotypes defined by the NMF algorithm with respect to the microdissected and activated stroma signature (Figure 1; Tables 1 and 2). The top 150 increased differentially expressed genes (DEGs) between the immune and non-immune classes (Table S4) were chosen as the seed genes to regenerate the immune subclasses in the external cohorts using the GenePattern module NMFConsensus, and then immune-activated and immune-exhausted subgroups were further separated via the NTP method.

Table 1

Summary of the detailed information of the enrolled bladder cancer cohorts

Dataset	Data array	Patients	Reference
TCGA-BLCA	RNA sequencing	408	https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Prostate%20Cancer%20(PRAD)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443
E-MTAB-4321	RNA sequencing	476	https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-4321/
IMvigor210	Illumina HiSeq 2500	348	http://research-pub.gene.com/IMvigor210CoreBiologies/
GSE32894	Illumina HumanHT-12 V3.0 expression beadchip	308	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32894
GSE83586	Affymetrix Human Gene 1.0 ST Array	307	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE83586
GSE87304	Affymetrix Human Exon 1.0 ST Array	305	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE87304
GSE128702	Affymetrix Human Exon 1.0 ST Array	256	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE128702
GSE13507	Illumina human-6 v2.0 expression beadchip	164	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE13507
GSE129871	Illumina HiSeq 2000 (Homo sapiens)	158	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE129871
GSE120736	Illumina HumanHT-12 V4.0 expression beadchip	145	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE120736
GSE39016	Affymetrix Human Exon 1.0 ST Array	141	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39016
GSE128701	Affymetrix Human Exon 1.0 ST Array	136	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE128701
GSE124035	Affymetrix Human Exon 1.0 ST Array	133	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE124305
GSE86411	Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip	132	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE86411
GSE48276	Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip	116	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48276
GSE128192	Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip	112	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE128192
GSE31684	Affymetrix Human Genome U133 Plus 2.0 Array	93	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE31684
GSE134292	Illumina HiSeq 4000 (Homo sapiens)	80	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE134292
GSE93527	Affymetrix Human Transcriptome Array 2.0	79	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE93527
E-MTAB-1803	Affymetrix GeneChip Human Genome U133 Plus 2.0	70	https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-1803/
GSE69795	Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip	61	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE69795

In the GSE32894 cohort, 60.7% (187/308) of patients belonged to the non-immune class, with the lower enrichment of immune-associated signatures; as for the remaining 121 patients, compared with the signatures of stromal enrichment, 42 patients divided into the immune-activated subgroup and 79 belonged to the immune-exhausted subgroup. High enrichment scores of TITR, MDSC, WNT/TGFβ, TGFβ-1-activated, and C-ECM signatures were observed in the immune-exhausted subgroup (all p < 0.01; Table 3; Figure S5).

Table 3

The distribution of three newly defined immunophenotypes in all enrolled cohorts

Dataset	No. of patients	Immunophenotype distribution, n (%)
		Immune activated	Immune exhausted	Non-immune
TCGA-BLCA	408	45 (11.03)	110 (26.96)	253 (62.01)
E-MTAB-4321	476	74 (15.55)	111 (23.32)	291 (61.13)
IMvigor210	348	85 (24.43)	142 (40.8)	121 (34.77)
GSE32894	308	42 (13.64)	79 (25.65)	187 (60.71)
GSE83586	307	63 (20.52)	93 (30.29)	151 (49.19)
GSE87304	305	59 (19.34)	85 (27.87)	161 (52.79)
GSE128702	256	72 (28.13)	88 (34.38)	96 (37.50)
GSE13507	164	23 (14.02)	36 (21.95)	105 (64.02)
GSE129871	158	26 (16.46)	27 (17.09)	105 (66.46)
GSE120736	145	21 (14.48)	35 (24.14)	89 (61.38)
GSE39016	141	16 (11.35)	31 (21.99)	94 (66.67)
GSE128701	136	42 (30.88)	34 (25.00)	60 (44.12)
GSE124035	133	32 (24.06)	54 (40.6)	47 (35.34)
GSE86411	132	22 (16.67)	36 (27.27)	74 (56.06)
GSE48276	116	24 (20.69)	29 (25.00)	63 (54.31)
GSE128192	112	26 (23.21)	36 (32.14)	50 (44.64)
GSE31684	93	14 (15.05)	34 (36.56)	45 (48.39)
GSE134292	80	13 (16.25)	16 (20.00)	51 (63.75)
GSE93527	79	13 (16.46)	15 (18.99)	51 (64.56)
E-MTAB-1803	70	13 (18.57)	19 (27.14)	38 (54.29)
GSE69795	61	9 (14.75)	19 (31.15)	33 (54.10)

In the other cohorts, we also replicated the three immunophenotypes, and the results are displayed in Table 3 and Figures S5 and S6. In these cohorts, the distribution of immune-activated subgroups ranged from 11.3% to 30.9%, while the proportion of immune-exhausted subgroups ranged from 17.1% to 40.8%. We also observed increased scores for the immune enrichment signature and immune signaling signature in the 18 validation cohorts, as well as the other immunocyte signatures. As expected, the immune-exhausted subgroup showed increased enrichment scores for Tregs, TITR, MDSC, WNT/TGFβ, and C-ECM signatures. With the combined results from Tables 1, 2, 3, and S4 and Figures S5 and S6, our results suggest that the NMF and NTP algorithms could stably and precisely divide bladder patients into immune-activated, immune-exhausted, and non-immune phenotypes. The specific immune characteristics can reappear in all the enrolled bladder cancer cohorts.

Favorable response to anti-PD-L1 therapy was observed in the immune-activated subgroup

To evaluate the response to immunotherapy of the bladder cancer patients with newly defined immunophenotypes, we collected the gene expression profile and clinical outcomes of 348 patients from the IMvigor210 cohort, a large phase II trial investigating the clinical response of PD-L1 blockade with atezolizumab. The observed endpoint is the response to anti-PD-L1 treatment, including complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). The objective response rate (ORR) included patients with CR and PR, while the disease control rate (DCR) included patients with CR, PR, and SD. Patients were first separated into immune (237/348, 65.2%) and non-immune (121/348, 34.8%) classes by the 150 DEGs generated from the TCGA-BLCA training cohort. Then, 85 patients (85/348, 24.4%) with a negative result of stromal activation prediction in the immune class were recognized as the immune-activated subgroup, while the remaining 142 patients (142/348, 40.8%) belonged to the immune-exhausted subgroup (Figure 4A). It is obvious that patients in the immune class showed an abundance of immune signatures, including predefined immune enrichment score, immune cell subset score, immune signaling score, as well as the signatures of several immunocytes (all p < 0.05). For the separation of the immune-activated and immune-exhausted subgroups, MDSC, TITR, Tregs, WNT/TGFβ, TGFβ-1, and C-ECM signatures showed a higher enrichment score in the immune-exhausted subgroup, consistent with the characteristics in the TCGA-BLCA training cohort. We observed that more than half of the patients in the immune-activated subgroup benefited from anti-PD-L1 therapy; the ORR was 24.66%, and the DCR was 52.06%. However, most patients in the immune-exhausted subgroup showed PD results after anti-PD-L1 treatment (p = 0.028; Figure 4B). Moreover, if we only focused on the treatment results of CR and PD, we found that patients in the immune-activated subgroup showed a higher response rate than other subgroups (20.5% versus 15.1% versus 4.8%; Figure 4C). In general, the results in Figure 4 show that the newly defined classifier can distinguish the appropriate patients to receive anti-PD-L1 therapy, of whom in the immune-activated subgroup.

Figure 4

The immunophenotypes in the IMvigor210 cohort reflected the different responses to anti-PD-L1 immunotherapy

(A) Reappearing three immunophenotypes in IMvigor210 cohort. (B) The distribution of the best confirmed overall response to anti-PD-L1 treatment in the three immunophenotypes. (C) The distribution of complete response and progressive disease in the three immunophenotypes. Non-IM, non-immune class; IM-Act, immune-activated subgroup; IM-Exh, immune-exhausted subgroup.

Immune-activated subgroup shows a favorable prognosis and benefits more from anti-PD-1 therapy

To verify the excellent effectiveness of immune therapy in the immune-activated subgroup, we used subclass mapping (SubMap) analysis to compare the gene expression distribution of the responders in the IMvigor210 cohort and the MD Anderson melanoma cohort, with the immune-activated subgroups of TCGA-BLCA, GSE32894, and E-MTAB-1803 cohorts. We found that patients belonging to immune-activated subgroups contained a more similar gene distribution to the responders than non-responders of anti-PD-L1 treatment in the IMvigor210 cohort but not those patients in the non-immune-activated subgroup (all Bonferroni corrected p < 0.05) (Figure 5A). Moreover, we also observed consistent results in the comparison with the MD Anderson melanoma cohort, and patients in the immune-activated subgroup can profit more from anti-PD-1 treatment (all Bonferroni-corrected p < 0.05) (Figure 5B). In addition, we revealed that patients belonging to the immune-activated subgroup illustrated a better prognosis, while patients separated into the immune-exhausted subgroup suffered from a poor prognosis (Figure 5C). Taken together with the results from Figure 5, we concluded that patients in the immune-activated subgroup can profit more from anti-PD-1 or anti-PD-L1 therapy and had the longest average overall survival time.

Figure 5

Predicting the response to immune checkpoint blockade therapy and revealed diverse overall survival outcomes in the three immunophenotypes

(A) Prediction results of the response to anti-PD-L1 therapy. (B) Prediction results of the response to anti-CTAL-4 and anti-PD-1 therapy. (C) Different overall survival outcomes in the three immunophenotypes.

Comparison of the three immunophenotypes with the proposed molecular subtypes and reappearance in pan cancer

The newly defined immunophenotypes are the supplement to the proposed molecular subtypes, which can help to precisely depict the molecular characters of bladder cancer patients. Six immune molecular features were previously defined by Thorsson et al., reflecting wound healing (C1), IFN-γ-dominant (C2), inflammatory (C3), lymphocyte-depleted (C4), immunologically quiet (C5), and TGF-β-dominant (C6) specific features. We found that most patients in the immune-activated subgroup (75.00%) belonged to the IFN-γ group, which is associated with a strong CD8 signal and a high proliferation rate, and that approximately 68.22% of patients in the immune-exhausted subgroup also belonged to the IFN-γ group (Figure S7A). Another classical molecular classifier was defined by Kamoun et al.; luminal papillary, luminal nonspecified, luminal unstable, stroma-rich, basal/squamous, and neuroendocrine-like subtypes cover almost every aspect of bladder cancer. In the current study, we revealed that most immune-exhausted patients and half of the immune-activated patients belonged to the basal/squamous classes (Figure S7B), with TP53 being the most frequently mutated gene, consistent with previous results. Unquestionably, the Ba/sq subclass was associated with poor prognosis, and patients in the immune-exhausted subgroup also exhibited poor prognosis. Along with the results shown in Figure S7, we can conclude that the newly defined bladder cancer immunophenotypes have a high consistency with the proposed molecular subtypes.

We also tried to reproduce the three immunophenotypes in other cancer types of the TCGA dataset. With our recently developed R package MOVICS, we first selected the specific DEGs for each immunophenotype and then distinguished the non-immune, immune-activated, and immune-exhausted subgroups in pan cancer cohorts. Finally, we successfully distinguished the three immunophenotypes in genitourinary system tumors, including papillary renal cell carcinoma (KIRP), prostate cancer (PRAD), and testicular germ cell tumor (TGCT) (Figure S8). Patients in the immune-exhausted subgroup showed a poor prognosis (Figure S9), and the prognosis results of PRAD are consistent with our recent work. In addition, we also found that the classification system is suitable for adrenocortical carcinoma (ACC), brain lower-grade glioma (LGG), mesothelioma (MESO), stomach adenocarcinoma (STAD), and uveal melanoma (UVM) (Figures S8 and S9).

Dimensionality reduction of the 150 DEGs for the distinguishment of immune and non-immune classes

The definition of the three immunophenotypes was a two-step process. First, we separated the patients into immune and non-immune classes by NMF consensus via the 150 DEGs from the TCGA training cohort. Second, we partitioned patients in the immune class into immune-activated and immune-exhausted subtypes by the NTP method based on the expression of 47 stromal-associated genes defined by Moffitt et al. The smaller the panel is, the easier it is to implement clinically. Therefore, we used the random forest algorithm to generate a predictive model of immune and non-immune classes via a smaller number of genes. A total of 21 genes were selected (Table S5), the functions of which were enriched in several immune-associated biological pathways (Figure S10A). The 21 genes can 100% predict the distribution of patients in two immune classes in the TCGA-BLCA training cohort, with the predicted accuracy of 0.85 in the E-MTAB-1803 validation cohort, 0.97 in the GSE32894 validation cohort, and 0.93 in the IMvigor210 validation cohort (Figures S10B and S10C).

Discussion

Bladder cancer is a heavy health burden worldwide, especially in Europe and North America. There are approximately 550,000 new cases and 200,000 specific deaths caused by bladder cancer each year. The incidence rate of bladder cancer varies around the world, with the highest rate in southern Europe and the lowest rate in Middle Africa. In the United States, bladder cancer is the 6th most common type of cancer, with approximately 81,000 new cases and 18,000 new deaths reported in 2020. Most patients are diagnosed with early-stage bladder cancer, also known as NMIBC, in which processes do not involve the muscle layer and can be removed easily through transurethral resection (TUR) or with intravesical therapy with bacillus Calmette-Guérin (BCG) or by another chemotherapeutic treatment.32, 33, 34 Recurrence is extremely common in NMIBC; approximately 70% of patients will suffer from the health burden of bladder cancer again within 10 years, and one-third of patients will progress to an advanced stage, also called MIBC. The standard care for MIBC is radical cystectomy with or without neoadjuvant chemotherapy or chemoradiation. Even after treatment, almost 50% of MIBC patients will experience recurrence and death from the metastatic stage within 3 years.

Several studies investigated the diagnostic markers, prognostic signatures, or therapeutic targets for malignancy of tumors based on the TME, as well as in bladder cancer. BCG is the earliest immune therapy approved for bladder cancer treatment, which can stimulate an immunologic reaction and induce proinflammatory cytokines and direct cell-to-cell cytotoxicity. BCG is still the standard therapy for NMIBC, which reflects that bladder cancer patients could benefit from immunotherapy. The blockade of immune checkpoints is also applied in the treatment of bladder cancer. Two PD-1 inhibitors (pembrolizumab and nivolumab) and three PD-L1 inhibitors (atezolizumab, durvalumab, and avelumab) have been approved by the US Food and Drug Administration (FDA) for the treatment of bladder cancer (http://www.accessdata.fda.gov/scripts/cder/daf/). In the IMvigor210 clinical study, atezolizumab was used to block PD-L1: the ORR for the IMvigor210 cohort 1 was only 23%, and it was only 15% for cohort 2., The ORR values for nivolumab and durvalumab are similar, ranging from 17% to 24.4%.39, 40, 41, 42 Therefore, a full understanding of the immunophenotypes of bladder cancer is essential and could serve as guidance in choosing patients to receive the appropriate immunotherapy.

For the enrolled 4,028 bladder cancer patients, we used NMF and NTP algorithms to investigate a robust classification of three immunophenotypes. First, we identified the immune-activated subgroup, immune-exhausted subgroup, and non-immune class in the TCGA-BLCA training cohort. The top weighted 150 exemplar genes selected from the immune module represent the immune features of bladder cancer patients and further divide the entire cohort into immune and non-immune classes. The other 150 DEGs among the immune and non-immune classes were extracted as the input profile for the reappearance of the classification in validation cohorts. Stromal activation status was evaluated via the NTP algorithm to further separate the immune class into immune-activated and immune-exhausted subgroups. The features of these three immunophenotypes were illuminated by several verified signatures of immunocytes or immune signaling pathways. Patients in the immune classes showed highly enriched signatures of T cells, B cells, IFN, and CYT (L.Q.M. Chow et al., 2016, J. Clin. Oncol., abstract),, while the exhausted subgroup also showed increased signatures of TITR, WNT/TGF-β, TGF-β1-activated, and C-ECM signatures;45, 46, 47 this was not the case for the immune-activated subgroup. B Based on our results, we revealed that only approximately 11% to 30.9% of bladder cancer patients belong to the immune-activated subgroup, and can benifit from immunotherapy.

Clinical outcome is an important factor on which we focused with respect to the newly defined immunophenotypes. With the clinical information of the TCGA-BLCA, GSE32894, and E-MTAB-1803 cohorts, we found that patients belonging to the immune-activated subgroup exhibited the most favorable overall survival, while patients belonging to the immune-exhausted subgroup showed the worst clinical outcome of a shortened overall survival time. Immune exhaustion, which is mostly focused on the exhaustion of T cells, is reflected by altered inflammatory and tissue microenvironments, lymphocytes, and inhibitory signals from cytokines. These alterations in the TME could lead to the escape from immune recognition by blocking immune checkpoints, which is related to unfavorable overall survival for patients. We predicted the potential response to immunotherapy of the bladder cancer patients by comparing their mRNA expression profiles with those of melanoma samples receiving anti-CTLA-4 or anti-PD-1 checkpoint therapy. As expected, patients in the immune-activated subgroup could benefit from anti-PD-1 therapy, but this was not the case for patients in the non-immune-activated subgroup, which included the combined immune-exhausted and non-immune classes.

To further understand the molecular diversity among the three immunophenotypes, we compared the CNA, TMB, and gene mutations. Recent studies have reported the association of CNA with increased immune infiltration and the outcome of immune checkpoint blockade therapy., Patients in the immune class showed a lower CNA burden with gene deletion at the arm and focal levels. The association was verified twice by the deletion copy numbers of PD-1 and PD-L1, and CTLA-4 was positively associated with the decreased infiltration level of immunocytes. Tripathi et al. found that antigen presentation through the major histocompatibility complex (MHC) class I pathway was suppressed in tumors with high chromosomal instability, also known as high CNA, which plays a pivotal role in immune evasion. In addition, Lu et al. revealed that patients treated with immune checkpoint blockade therapy could receive a durable clinical benefit and achieve better survival if they exhibited a lower CNA burden. Gene mutation is another key component upon which we focused with respect to the three immunophenotypes. We extracted the specific mutant genes for each subgroup. The proportions of mutant TP53, TTN, PIC3CA, and RB1 were higher in the immune class than that in the non-immune class. Nusrat et al. (2019, J. Clin. Oncol., abstract) reported that colorectal cancer patients with mutant PIK3CA had higher median densities of CD3+ and CD8+ cells, as well as a high rate of clinical benefit from immunotherapy (50% versus 8.6%). Furthermore, we observed a high rate of ERBB2 mutations in the immune-activated subgroup. ERBB2 amplification or overexpression is a biomarker of anti-ERBB2-targeted therapy in breast cancer, while the activity of trastuzumab and infiltration of immune cells can both be regulated by the activated function of ERBB2 via the inducing of CCL2 and PD-1 ligands. The V659E mutation of the ERBB2 gene has also been reported to be associated with altered sensitivity to afatinib and lapatinib treatment in vitro., The highest mutation proportion of EP300 was observed in the immune-exhausted subgroup, and the mutation leads to the increased expression of EP300. The function of effector T cell responses to tumor antigens can be dampened by the Tregs and cause the immunosuppressive status of the TME; this is associated with poor prognosis. Tregs mediated immune suppressive function was reported to be impaired by the inhibition of EP300 as well; therefore, it is considerable to combine the immune checkpoint blockade therapy plus an EP300 inhibitor for the treatment of immune-exhausted patients. The limitation of the current study is the lack of the prospective validation of the newly defined immunophenotypes. In the future, we will use the immunophenotype classification system to choose appropriate bladder cancer patients to receive the immunotherapy.

We defined and validated a novel classifier among 4,028 bladder cancer patients to separate them into immune-activated, immune-exhausted, and non-immune subgroups. Patients in the immune-activated subgroup could benefit more from a single treatment of anti-PD-1 immunotherapy; for patients in the immune-exhausted subgroup, immune checkpoint blockade therapy plus a TGF-β inhibitor or an EP300 inhibitor might have greater efficacy. In summary, our novel classifier provides a better understanding of enhanced immunotherapy to treat bladder cancer patients.

Materials and methods

Bladder cancer patient cohorts

For the TCGA-BLCA cohort, we obtained the level three gene expression profiles of 408 patients from the TCGA Data Portal (https://www.cancer.gov/tcga/), and only genes expressed in at least 50% of the samples were retained for the subsequent analyses. For further external validation, we enrolled the Gene Expression Omnibus (GEO): GSE32894, GSE83586, GSE87304, GSE128702, GSE13507, GSE129871, GSE120736, GSE39016, GSE128701, GSE124035, GSE86411, GSE48276, GSE31684, GSE134292, GSE93527, and GSE69795 cohorts; the gene expression profiles and clinical information were collected from the GEO (https://www.ncbi.nlm.nih.gov/geo/). For the E-MTAB-4321 and E-MTAB-1803 cohorts, the gene expression profiles and paired clinical features were downloaded from ArrayExpress (https://www.ebi.ac.uk/arrayexpress/). Detailed information of these datasets is displayed in Tables 2 and 3.

Bioinformatic analyses

The mRNA expression profiles of the TCGA-BLCA training cohort were microdissected using the NMF algorithm. The identification of an immune class was achieved as reported by Sia et al. The immune module was selected, with the gathering of patients with the highest immune enrichment score, which was estimated by ssGSEA. Then, the top 150 exemplar genes with the highest weight in the immune module were extracted as the key genes to dichotomize the immune and non-immune classes, which were further modified by the MDS random forest method. In addition, immune-activated and immune-exhausted subgroups were recognized by the stromal activation signature by the NTP method. To depict the characteristics of these three immunophenotypes, several published immune gene sets were manually collected to reflect the immune status, and the scores of each gene set for each patient were also calculated via ssGSEA (Table S3). The different genetic types among the immune and non-immune classes were evaluated, including TIL abundance, PD-L1 expression, CNAs, TMB, neoantigens, and mutant genes. To recreate the immunophenotypes in validation cohorts, the expression profiles of the top 150 DEGs were used to dichotomize the immune classes via the NMF consensus method, and the immune class was subsequently divided into the activated and exhausted subgroups. To predict the immunotherapy response of bladder cancer patients, we collected the expression profile of the IMvigor210 cohort, along with the response results of anti-PD-L1 treatment from IMvigor210CoreBiologies. The MD Anderson melanoma cohort that received anti-CTLA-4 or anti-PD-1 therapy was also considered for the prediction of immunotherapy response. SubMap analysis was used to compare the similarity of gene expression between the newly defined bladder cancer immunophenotypes and the responders of anti-CTLA-4, anti-PD-1, and anti-PD-L1 therapy in the IMvigor210 cohort or MD Anderson melanoma cohort. We developed a predictive model to estimate immune classes based on the random forest algorithm using R package varSelRF. To this end, we obtained the DEGs using R package limma between immune and non-immune classes in the TCGA training cohort to select informative genes, which served as input to random forest model (foldchange > 2 or < 0.5 with false discovery rate [FDR] < 0.05). A backward elimination procedure was applied to search for optimal markers where an out-of-bag (OOB) error was used as the minimization criterion. In this procedure, variable elimination was conducted by setting the dropping fraction of each iteration at 0.2, and the gene combination was identified to develop the final predictive model when random forest reached the smallest OOB error rate. The predictive performance was assessed in TCGA, EMATB-1803, GSE32894, and IMvigor210 cohorts using receiver operating characteristic (ROC) curves. The details of the enrolled cohorts, as well as the specific method used for each step, are provided in the Supplemental methods.

Statistical analysis

For normally distributed continuous data, t test was used for the comparisons, while non-normally distributed data were compared via Wilcoxon rank-sum test. The comparisons between factors with more than two groups were performed by Kruskal-Wallis test. The difference in survival between different groups was analyzed by Kaplan-Meier plots and log-rank test. Chi-square test was used to illustrate the correlations between newly defined immunophenotypes and proposed molecular subtypes. Two-side p value less than 0.05 was considered as the statistically significant difference. All analyses were performed by GenePattern and R version 4.0.2 (http://www.r-project.org).

Availability of data and materials

All data used in this work can be acquired from the GDC portal (https://portal.gdc.cancer.gov/), GEO (https://www.ncbi.nlm.nih.gov/geo/), and ArrayExpress (https://www.ebi.ac.uk/arrayexpress/).