Comprehensive immune transcriptomic analysis in bladder cancer reveals subtype specific immune gene expression patterns of prognostic relevance.

Recent efforts on genome wide profiling of muscle invasive bladder cancer (MIBC) have led to its classification into distinct genomic and transcriptomic molecular subtypes that exhibit variability in prognosis. Evolving evidence from recent immunotherapy trials has demonstrated the significance of pre-existing tumour immune profiles that could guide treatment decisions. To identify immune gene expression patterns associated with the molecular subtypes, we performed a comprehensive in silico immune transcriptomic profiling, utilizing transcriptomic data from 347 MIBC cases from The Cancer Genome Atlas (TCGA). To investigate subtype-associated immune gene expression patterns, we assembled 924 immune response genes and specifically those involved in T-cell cytotoxicity and the Type I/II interferon pathways. A set of 157 ranked genes was able to distinguish the four subtypes in an unsupervised analysis in an original training cohort (n=122) and an expanded, validation cohort (n=225). The most common overrepresented pathways distinguishing the four molecular subtypes, included JAK/STAT signaling, Toll-like receptor signaling, interleukin signaling, and T-cell activation. Some of the most enriched biological processes were responses to IFN-γ, antigen processing and presentation, cytokine mediated signaling, hemopoeisis, cell proliferation and cellular defense response in the TCGA cluster IV. Our novel findings provide further insights into the association between genomic subtypes and immune activation in MIBC and may open novel opportunities for their exploitation towards precise treatment with immunotherapy.

INTRODUCTION

Urothelial bladder cancer (UBC) is the fifth most common cancer worldwide [1] and is one of the most management intensive cancers in North America [2]. Although the majority of incident cases of UBC are non-invasive at presentation, muscle invasive bladder cancer (MIBC) represents very aggressive disease with rapid progression to metastases [3] and poor overall survival despite intensive local and systemic therapy. Current standards for localized MIBC include radical cystectomy with or without perioperative cisplatin-based chemotherapy [3]. Unfortunately, many suffer early disease recurrence and, despite palliative chemotherapy, median survival rates are generally less than one year [4]. The optimal management of patients with higher risk UBC is ambiguous with a significant need for better prediction tools and enhanced therapeutics [5].

MIBCs are highly heterogeneous tumours. Recent investigations based on molecular profiling of specimens from large UBC cohorts have led to their classification into molecular subtypes that display distinct genomic and transcriptomic features, resembling those seen in breast cancer [3], [6-8]. Interestingly, these subtypes may exhibit distinct associations with treatment response and survival [8, 9]. Although different groups have classified UBC into two [8], three [3], four [6], or five [7] subtypes, there is a consensus that the top separation occurs as the basal and luminal subtypes [10]. Basal tumours, enriched with EGFR and hypoxia-inducible factor 1 expression, are often metastatic at presentation, possess squamous and sarcomatoid histological features, and have epithelial-to-mesenchymal transition cell biomarkers [11]. In comparison, luminal cancers have papillary features and commonly FGFR3, ERBB2, and ERBB3 activating mutations [11]. The Cancer Genome Atlas Network (TCGA) bladder analysis working group classified bladder tumours into four clusters named I, II, III, and IV [6]. Clusters I and II correspond to the luminal subtype, while III and IV represent the basal subtype [12]. Tumours in Cluster I are enriched in FGFR3 overexpression due to mutations and amplification and show better overall survival, whereas those in cluster II, designated “p53-like” tumours, express active p53 gene signatures and are resistant to neoadjuvant cisplatin-based combination chemotherapy [3]. Cluster IV shares features with the claudin-low subtype of breast cancer, express immune checkpoint molecules, and were actively immunosuppressed, despite having an enriched immune gene signature [13]. In particular, Kardos et al. [13] demonstrated that immune infiltration was not correlated with predicted neoantigen burden, but from unopposed NF-kB activity from downregulated PPARγ signaling.

Given the urgent need of alternative approaches in MIBC treatment, there has been a growing interest in immunotherapies, such as those targeting the immune checkpoints: CTLA-4, PD-L1, and PD-1 [14]. Atezolizumab, a PD-L1 inhibitor, has been recently approved by the FDA for bladder cancer that progressed during or following chemotherapy [15]. Evolving evidence based on the success of immune checkpoint blockade therapies in melanoma and non-small cell lung cancer has confirmed the significance of the pre-treatment tumour immune state as a strong prognostic and response predictive indicator [14, 16]. An important feature, key to the success of immunotherapy, is the spatial organization of cytotoxic CD8+ tumour infiltrating lymphocytes (TILs) in the epithelial and stromal compartments and their activation status [17]. Higher density of CD3+ and CD8+ TILs have been associated with increased disease-free and overall survival in melanoma, head and neck, breast, bladder, urothelial, ovarian, colorectal, prostatic, and lung cancer; however, their activation status determines their prognostic significance in most cancers [18, 19]. In particular the expression of interferons (IFNs), which play a central role in anti-tumour immune responses, are emerging as prognostic and predictive biomarkers of both chemotherapy and immunotherapy [20]. Higher infiltration with CD4+ and regulatory subsets of TILs and higher CD68 to CD3 ratios are associated with poor prognosis in bladder cancer [21-23]. In particular PD-L1, IDO, FOXP3, TIM3, and LAG3 are expressed in T-cell-inflamed, and β-catenin, PPAR-γ, and FGFR3 in non-T-cell-inflamed urothelial tumours [17]. Although the pre-treatment expression of PD-1/PD-L1 initially showed some predictive value, it has recently failed to perform as a good biomarker in the recent clinical trials due to their transient nature of expression [23-25].

As reported in other cancers sites, it is likely that the pre-existing tumour immune landscape in UBC could be an additive determinant of response to chemotherapy as well as immune-based therapies leading to more precise prognostication, patient stratification, and informed treatment decisions [26]. To our knowledge there are no previous studies in MIBC that have evaluated the association between immune transcriptomic alterations, specifically those mediated by IFNs and cytotoxic pathway genes, and their potential associations with their distinct molecular sub-populations. In the current study, we performed a comprehensive in silico immune transcriptomic profiling of MIBC using the publicly available TCGA global transcriptomic datasets in order to determine whether the known molecular subtypes of MIBC are associated with specific immune gene signatures. The findings from our study may not only provide insights into the association between genomic subtypes and immune activation, but may also open novel opportunities for improving the management of MIBC.

RESULTS

We aimed to determine whether the previously defined four TCGA MIBC clusters exhibit differences in their immune gene expression patterns that could be of potential significance in informing treatment decisions for immunotherapies and other combinatorial treatment approaches.

Clinicopathological features of TCGA MIBC cohort

The TCGA cohort as reported earlier, consisted of chemotherapy-naïve, muscle-invasive, high-grade urothelial tumors (T2-T4a, Nx, Mx) [6]. Inclusion criteria reviewed by five expert genitourinary pathologists involved: tumour nuclei ≥ 60% of total, ≤20% tumour necrosis in the specimen, and variant histology ≤50% [6].

Immune gene expression patterns across MIBC clusters

First, the 122 samples previously used to identify the four clusters by TCGA [6] were treated as a discovery group to determine immune gene expression profiles across clusters. A total of 377 genes derived from the NanoString™ panel discriminating among the clusters were identified using a feature selection technique (Supplementary Table 1). The performance of these genes to accurately distinguish the four TCGA clusters was evaluated by clustering of cohort 1 (Figure 1). This set of genes was then used to assign samples in the validation set to the four clusters. Similar to cohort 1, the genes were able to distinguish the four clusters in cohort 2 by supervised and unsupervised analysis (Figure 2A and 2B). Similar unsupervised analysis was done using the top 5% of genes derived from all immune panels (n=157) (Supplementary Table 2) on both cohorts (Figure 3A and 3B). A recent updated analysis of the current TCGA bladder tumour cohort shows that clusters I-IV remained stable [28], supporting our classification approach in cohorts 1 and 2.

Figure 1

Distinct immune gene expression levels in cohort 1 (n=122) between the four TCGA bladder cancer subtypes based on the top 20% (377 NanosString panel genes) using the feature selection algorithm. Red indicates high expression, and blue indicates low expression.

Figure 2

Cohort 2 (n=225) assigned to clusters based on Euclidian distance to the cluster centroids generated from the cohort 1 (n=122)

Supervised (A) and unsupervised (B) analysis based on the samples and 377 NanoString™ panel genes.

Figure 3

Unsupervised analysis of both the cohort 1 (A) and cohort 2 (B) using the top 5% (n=157) genes. Unsupervised grouping shows gradient of under expression in cluster I to overexpression in cluster IV.

Differential pre-existing expression patterns of interferon associated genes

The four cluster patterns were also noted for the top 20% ranking IFN-γ associated genes upon hierarchical clustering analysis in both training and validation cohorts (Figure 4A and 4B). A gradient of under-expression of IFN-γ associated genes in cluster I to overexpression in cluster IV is observed in both. A similar pattern was also noted in the top 20% ranking IFN-α (Figure 5A and 5B) and cytotoxic genes (Figure 6A and 6B). Most importantly, key IFN response genes and downstream T-cell recruiting target chemokine genes, CXCL9, CXCL10, and CXCL11, and their common receptor CXCR3, showed increased expression in clusters III and IV. Similarly, others in the list included the key players in IFN response such as IFITM2, CCL5, IRF4 and others.

Figure 4

Supervised heat map of top 20% of IFN-γ associated pathway genes in both cohort 1 (A) and cohort 2 (B).

Figure 5

Supervised heat map of top 20% of IFN-α associated pathway genes in both discovery (A) and validation (B) groups.

Figure 6

Supervised heat map of top 20% of cytotoxic associated pathway genes in both cohort 1 (A) and cohort 2 (B).

Antigen processing pathways are overrepresented in T-cell inflamed MIBC clusters

We determined the Gene Ontology functional annotations of the differentially expressed genes that distinguished the four clusters, using both the 377 and 157 genes as input gene lists. Using the overrepresentation statistic in PANTHER, we calculated the probability of highly populated protein classes and gene ontology classes among the two gene lists (Table 2a and 2b). The most enriched GO biological processes in the 377-gene list were response to IFN-γ, antigen processing and presentation, cytokine mediated signaling, hemopoeisis, cell proliferation and cellular defense response (Supplementary Table 3). The top five overrepresented pathways included JAK/STAT signaling pathway, Toll-like receptor signaling pathway, interleukin signalling pathway, and T-cell activation (Figure 7). Interestingly, similar analysis using the top ranking 157 genes as input list, revealed only the B-cell activation, T-cell activation, and inflammation mediated by chemokine and cytokine signaling pathways as the only three overrepresented pathways. Response to IFN-γ, hemopoiesis, macrophage activation and cell proliferation were the most overrepresented GO biological processes in the top ranked 157 genes (Supplementary Table 4).

Table 2a

Analysis Type:	PANTHER Overrepresentation Test (release 20160715)
Annotation Version and Release Date:	PANTHER version 11.1 Released 2016-10-24
Analyzed List:	377genes (Homo sapiens)
Reference List:	Homo sapiens (all genes in database)
Bonferroni correction:	TRUE
Bonferroni count:	158
PANTHER Pathways	Homo sapiens - REFLIST (20972)	377genes (385)	377genes (expected)	377genes (over/under)	377 genes (fold Enrichment)	377genes (P-value)
JAK/STAT signaling pathway (P00038)	17	6	0.31	+	19.23	1.50E-04
Toll receptor signaling pathway (P00054)	60	21	1.1	+	19.07	5.00E-18
Interleukin signaling pathway (P00036)	98	23	1.8	+	12.78	4.53E-16
T cell activation (P00053)	96	22	1.76	+	12.48	4.03E-15
B cell activation (P00010)	72	15	1.32	+	11.35	1.83E-09
p38 MAPK pathway (P05918)	42	7	0.77	+	9.08	2.49E-03
Apoptosis signaling pathway (P00006)	122	20	2.24	+	8.93	5.24E-11
Inflammation mediated by chemokine and cytokine signaling pathway (P00031)	261	36	4.79	+	7.51	3.53E-18
Blood coagulation (P00011)	47	6	0.86	+	6.95	4.23E-02
VEGF signaling pathway (P00056)	72	8	1.32	+	6.05	1.08E-02
Ras Pathway (P04393)	76	8	1.4	+	5.73	1.56E-02
CCKR signaling map (P06959)	173	17	3.18	+	5.35	6.02E-06
p53 pathway (P00059)	88	8	1.62	+	4.95	4.19E-02
Angiogenesis (P00005)	176	16	3.23	+	4.95	4.23E-05
Integrin signalling pathway (P00034)	192	17	3.52	+	4.82	2.59E-05
EGF receptor signaling pathway (P00018)	139	11	2.55	+	4.31	1.07E-02
Gonadotropin-releasing hormone receptor pathway (P06664)	235	15	4.31	+	3.48	6.36E-03
Unclassified (UNCLASSIFIED)	18333	232	336.55	-	0.69	0.00E+00

Table 2b

Analysis Type:	PANTHER Overrepresentation Test (release 20160715)
Annotation Version and Release Date:	PANTHER version 11.1 Released 2016-10-24
Analyzed List:	157 genes (Homo sapiens)
Reference List:	Homo sapiens (all genes in database)
Bonferroni correction:	TRUE
Bonferroni count:	158
PANTHER Pathways	Homo sapiens - REFLIST (20972)	157 genes Input (158)	157 genes Input (expected)	157 genes Input (over/under)	157 genes Input (fold Enrichment)	157 genes Input (P-value)
JAK/STAT signaling pathway (P00038)	17	3	0.13	+	23.42	4.94E-02
Interleukin signaling pathway (P00036)	98	9	0.74	+	12.19	1.20E-05
B cell activation (P00010)	72	6	0.54	+	11.06	3.25E-03
Toll receptor signaling pathway (P00054)	60	5	0.45	+	11.06	1.62E-02
T cell activation (P00053)	96	7	0.72	+	9.68	1.55E-03
Inflammation mediated by chemokine and cytokine signaling pathway (P00031)	261	12	1.97	+	6.1	1.34E-04
Integrin signalling pathway (P00034)	192	8	1.45	+	5.53	1.86E-02
Unclassified (UNCLASSIFIED)	18333	96	138.12	-	0.7	0.00E+00

Figure 7

Bar graph depicting distribution of fold enrichment levels of biological pathways defined by PANTHER based analysis in the 377 genes that show differential expression patterns in the four TCGA MIBC clusters

The enriched categories were obtained upon analysis using the statistical overrepresentation test defined by PANTHER tool [27].

DISCUSSION

Evolving research from correlative studies as well as clinical trials, including those for UBC, have emphasized the value of the pre-existing tumour immune state as a predictor of response to treatment and survival [29, 30]. Furthermore, UBC is associated with a comparatively high mutational burden [31], which could potentially contribute to increased immunogenicity making them more susceptible to novel immunotherapy-based approaches. In order to gain a better understanding of the pre-existing tumour immune landscape in UBC, we conducted a comprehensive in silico immune transcriptomic profiling of the MIBC tumours from the TCGA database.

Four distinct molecular subtypes in MIBC were defined recently based on the TCGA MIBC genome wide profiling datasets [6]. Indeed, variability in subtype nomenclature has been reported [10], which could be attributed to heterogeneity in tissue samples in addition to several other factors such as inclusion of NMIBC cases in classification schemes. However, since the TCGA bladder cohort is enriched for MIBC and gene expression based clusters have been well defined, we specifically used this cohort to address our question on immune gene expression patterns associated with genomic alterations. Based on the distinct immune signature between clusters in cohort 1 (n=122), we were able to assign cohort 2 (n=225) into the associated TCGA clusters using the top 20% of ranked immune genes from the training cohort. Further analysis on the IFN-γ, IFN-α, and cytotoxic genes were then compared in both cohorts. All of these analyses revealed an increased expression of immune-associated genes in Cluster IV and underactive immune environment in Cluster I. Given that specific genetic alterations associate with these molecular subtypes, it seems that anti-tumour immune responses could be partly driven by oncogenic drivers.

Cluster I has been previously reported to show higher expression of FGFR3 via mutations, amplifications, and other mechanisms [3]. Interestingly this cluster showed a distinct underactive tumour immune state with reduced expression of IFN genes and genes associated with T-helper type-1 response. It is indeed intriguing that tumours with FGFR3 mutations or overexpression as per previous classification [6], show an increased overall survival, which contradicts the underactive immune state observed here. In contrast, cases in cluster IV showed the most dominant immune response amongst all four clusters. Tumours in this cluster show decreased expression of PPAR-γ and GATA3, and significantly increased expression of IFN and antigen presentation pathway genes, in addition to MHC class II genes and those involved in T-cell cytolytic activity. Previous reports have shown that based on broader classifications, cases in cluster IV belong to the basal subtype, which shows poor overall survival [8], [9]. One potential contributor to these associations is the increased expression of immune checkpoint molecules such as CD274 (PD-L1), IDO1 and the immunosuppressive IL6 in clusters III and IV that potentially lead to increased resistance to cytotoxic killing and poor response to treatment and ultimately poor survival. As previously shown, this cluster also shows higher levels of EMT related genes, indicating more aggressive tumour phenotype [33]. Our recent report demonstrating that higher PD-L1 expression in cancer cells leads to increased drug resistance upon activation by IFN-γ or PD-1 [32] supports this notion. It could thus be speculated that the IFN-γ secreted by activated T-cells, reflected by the increased expression of GZMA in these clusters, could be inducing PD-L1 on the cancer cells with further interaction between these leading to T-cell dysfunction. However, the mechanistic basis of these significant associations needs to be explored further. In other cancers such as melanoma, colorectal, and ovarian, higher expression of IFN pathway genes and of those representing an active immune response is associated with a favourable treatment outcome and overall survival. Furthermore, it is also possible that factors other than anti-tumour immune responses contribute to increased survival rates in tumours with FGFR3 mutations in MIBC.

Increased expression of MHC class II genes CD74, HLA-DMB, and HLA-DQA1 indicate higher tumour antigen processing by the antigen presenting cells in clusters III and IV. This was also confirmed by gene ontology-based analysis, which reflected a dominance of response to IFN-γ, antigen processing and presentation, cytokine mediated signaling, and cell proliferation, NK cell and macrophage activation, and B cell mediated immunity. These enrichments not only confirm the increased active anti-tumour immune response in clusters IV and some of cluster III but also indicate the immunogenic nature of these clusters that could be potentially be driven by higher mutational burden and recognition of immune cells.

Overall, our findings based on comprehensive immune transcriptomic analysis have significant implications in informing treatment decisions based on immune gene expression patterns. Specifically, since immune checkpoint blockade therapy has shown some promise in bladder cancer [15], near future biomarker driven clinical trials will benefit from these findings that emphasize appropriate patient stratification for treatment. Although recent reports have described the presence of T-cell inflamed and non-inflamed MIBC tumours [17, 33], no previous reports have identified associations between immune response and IFN-associated genes with the four molecular MIBC subtypes. Our study is limited by the fact that the TCGA dataset is enriched for MIBC and thus further validations in other cohorts need to be performed; however, these associations are timely and complement the ongoing and future clinical trials based on immune-based therapies. Clinical translation of our findings will most appropriately be addressed by validation of the most significant differentially expressed genes at both transcriptional and proteomic levels in retrospective and prospective pre-treatment bladder tumour specimens. Future investigations by integration of genomic alterations determined by exome and transcriptome sequencing data are key to identifying the genomic determinants of variability in immune response. Finally our study provides an improved understanding of the bladder cancer molecular subtype associated immune gene expression patterns and will significantly impact the design of novel immune based therapies.

MATERIALS AND METHODS

Patient tumour samples

The publicly available global transcriptomic sequencing (RNA-Seq) data from 412 MIBC cases, with the corresponding clinical information was downloaded from TCGA data portal (https://gdc-portal.nci.nih.gov/), now part of the National Cancer Institute’s Genetic Data Commons. The cohort was further divided into two cohorts for downstream analysis. For our training cohort (cohort 1) we used data from the previously defined 129 cases from TCGA that were divided into four clusters based on their integrated analysis of mRNA, miRNA, and protein data [6]. Since our objective was to define immune gene expression patterns in treatment naïve tumours, we excluded patients with any previous therapy. Thus a cohort of 122 MIBC cases, divided into four molecular clusters, was used for in silico immune profiling. The remaining 283 cases constituted the validation cohort (cohort 2) of which 225 had no prior therapy.

Design of immune pathway gene panel for in silico immune profiling

To investigate the presence of subtype associated immune signatures, we assembled a defined set of 924 immune related genes. This curated list (Table 1) primarily consisted of genes involved in IFN-α (97 genes), IFN-γ (200 genes), and cytotoxic (115 genes) pathways as defined by Gene Set Enrichment Analysis (GSEA) in combination with other immune genes. The nCounter PanCancer Immune profiling panel (722 genes), (http://www.nanostring.com/products/pancancer_immune/) was used to derive the immune response related genes.

Table 1

Custom designed immune gene panel of 924 genes, consisting of IFN-α and IFN-γ pathway genes from GSEA and immune response genes defined by the NanoString nCounter PanCancer immune Pathways Panel

Immune gene panel - NanoString nCounter PanCancer immune panel
A2M	C1R	CCND3	CD44	CLEC5A
ABCB1	C1S	CCR1	CD46	CLEC6A
ABL1	C2	CCR2	CD47	CLEC7A
ADA	C3	CCR3	CD48	CLU
ADORA2A	C3AR1	CCR4	CD5	CMA1
AICDA	C4BPA	CCR5	CD53	CMKLR1
AIRE	C5	CCR6	CD55	COL3A1
AKT3	C6	CCR7	CD58	COLEC12
ALCAM	C7	CCR9	CD59	CR1
AMBP	C8A	CCRL2	CD6	CR2
AMICA1	C8B	CD14	CD63	CREB1
ANP32B	C8G	CD160	CD68	CREB5
ANXA1	C9	CD163	CD7	CREBBP
APOE	CAMP	CD164	CD70	CRP
APP	CARD11	CD180	CD74	CSF1
ARG1	CARD9	CD19	CD79A	CSF1R
ARG2	CASP1	CD1A	CD79B	CSF2
ATF1	CASP10	CD1B	CD80	CSF2RB
ATF2	CASP3	CD1C	CD81	CSF3
ATG10	CASP8	CD1D	CD83	CSF3R
ATG12	CCL1	CD1E	CD84	CT45A1
ATG16L1	CCL11	CD2	CD86	CTAG1B
ATG5	CCL13	CD200	CD8A	CTAGE1
ATG7	CCL14	CD207	CD8B	CTCFL
ATM	CCL15	CD209	CD9	CTLA4
AXL	CCL16	CD22	CD96	CTSG
BAGE	CCL17	CD24	CD97	CTSH
BATF	CCL18	CD244	CD99	CTSL1
BAX	CCL19	CD247	CDH1	CTSS
BCL10	CCL2	CD27	CDH5	CTSW
BCL2	CCL20	CD274	CDK1	CX3CL1
BCL2L1	CCL21	CD276	CDKN1A	CX3CR1
BCL6	CCL22	CD28	CEACAM1	CXCL1
BID	CCL23	CD33	CEACAM6	CXCL10
BIRC5	CCL24	CD34	CEACAM8	CXCL11
BLK	CCL25	CD36	CEBPB	CXCL12
BLNK	CCL26	CD37	CFB	CXCL13
BMI1	CCL27	CD38	CFD	CXCL14
BST1	CCL28	CD3D	CFI	CXCL16
BST2	CCL3	CD3E	CFP	CXCL2
BTK	CCL3L1	CD3EAP	CHIT1	CXCL3
BTLA	CCL4	CD3G	CHUK	CXCL5
C1QA	CCL5	CD4	CKLF	CXCL6
C1QB	CCL7	CD40	CLEC4A	CXCL9
C1QBP	CCL8	CD40LG	CLEC4C	CXCR1
CXCR2	FOS	IFI27	IL19	IRAK4
CXCR3	FOXJ1	IFI35	IL1A	IRF1
CXCR4	FOXP3	IFIH1	IL1B	IRF2
CXCR5	FPR2	IFIT1	IL1R1	IRF3
CXCR6	FUT5	IFIT2	IL1R2	IRF4
CYBB	FUT7	IFITM1	IL1RAP	IRF5
CYFIP2	FYN	IFITM2	IL1RAPL2	IRF7
CYLD	GAGE1	IFNA1	IL1RL1	IRF8
DDX43	GATA3	IFNA17	IL1RL2	IRGM
DDX58	GNLY	IFNA2	IL1RN	ISG15
DEFB1	GPI	IFNA7	IL2	ISG20
DMBT1	GPR44	IFNA8	IL21	ITCH
DOCK9	GTF3C1	IFNAR1	IL21R	ITGA1
DPP4	GZMA	IFNAR2	IL22	ITGA2
DUSP4	GZMB	IFNB1	IL22RA1	ITGA2B
DUSP6	GZMH	IFNG	IL22RA2	ITGA4
EBI3	GZMK	IFNGR1	IL23A	ITGA5
ECSIT	GZMM	IGF1R	IL23R	ITGA6
EGR1	HAMP	IGF2R	IL24	ITGAE
EGR2	HAVCR2	IGLL1	IL25	ITGAL
ELANE	HCK	IKBKB	IL26	ITGAM
ELK1	HLA-A	IKBKE	IL27	ITGAX
ENG	HLA-B	IKBKG	IL28A	ITGB1
ENTPD1	HLA-C	IL10	IL29	ITGB2
EOMES	HLA-DMA	IL10RA	IL2RA	ITGB3
EP300	HLA-DMB	IL11	IL2RB	ITGB4
EPCAM	HLA-DOB	IL11RA	IL2RG	ITK
ETS1	HLA-DPA1	IL12A	IL3	JAK1
EWSR1	HLA-DPB1	IL12B	IL32	JAK2
F12	HLA-DQA1	IL12RB1	IL34	JAK3
F13A1	HLA-DQB1	IL12RB2	IL3RA	JAM3
F2RL1	HLA-DRA	IL13	IL4	KIR2DL1
FADD	HLA-E	IL13RA1	IL4R	KIR2DL3
FAS	HLA-G	IL13RA2	IL5	KIR3DL1
FCER1A	HMGB1	IL15	IL5RA	KIR3DL2
FCER1G	HRAS	IL15RA	IL6	KIR3DL3
FCER2	HSD11B1	IL16	IL6R	KIT
FCGR1A	ICAM1	IL17A	IL6ST	KLRB1
FCGR2A	ICAM2	IL17B	IL7	KLRC1
FCGR2B	ICAM3	IL17F	IL7R	KLRC2
FCGR3A	ICAM4	IL17RA	IL9	KLRD1
FEZ1	ICOS	IL17RB	ILF3	KLRF1
FLT3	ICOSLG	IL18	INPP5D	KLRG1
FLT3LG	IDO1	IL18R1	IRAK1	KLRK1
FN1	IFI16	IL18RAP	IRAK2	LAG3
LAIR2	MAPK3	NT5E	RELB	STAT2
LAMP1	MAPK8	NUP107	REPS1	STAT3
LAMP2	MAPKAPK2	OAS3	RIPK2	STAT4
LAMP3	MARCO	OSM	ROPN1	STAT5B
LBP	MASP1	PASD1	RORA	STAT6
LCK	MASP2	PAX5	RORC	SYCP1
LCN2	MAVS	PBK	RPS6	SYK
LCP1	MBL2	PDCD1	RRAD	SYT17
LGALS3	MCAM	PDCD1LG2	RUNX1	TAB1
LIF	MEF2C	PDGFC	RUNX3	TAL1
LILRA1	MEFV	PDGFRB	S100A12	TANK
LILRA4	MERTK	PECAM1	S100A7	TAP1
LILRA5	MFGE8	PIK3CD	S100A8	TAP2
LILRB1	MICA	PIK3CG	S100B	TAPBP
LILRB2	MICB	PIN1	SAA1	TARP
LILRB3	MIF	PLA2G1B	SBNO2	TBK1
LRP1	MME	PLA2G6	SELE	TBX21
LRRN3	MNX1	PLAU	SELL	TCF7
LTA	MPPED1	PLAUR	SELPLG	TFE3
LTB	MR1	PMCH	SEMG1	TFEB
LTBR	MRC1	PNMA1	SERPINB2	TFRC
LTF	MS4A1	POU2AF1	SERPING1	TGFB1
LTK	MS4A2	POU2F2	SH2B2	TGFB2
LY86	MSR1	PPARG	SH2D1A	THBD
LY9	MST1R	PPBP	SH2D1B	THBS1
LY96	MUC1	PRAME	SIGIRR	THY1
LYN	MX1	PRF1	SIGLEC1	TICAM1
MAF	MYD88	PRG2	SLAMF1	TICAM2
MAGEA1	NCAM1	PRKCD	SLAMF6	TIGIT
MAGEA12	NCF4	PRKCE	SLAMF7	TIRAP
MAGEA3	NCR1	PRM1	SLC11A1	TLR1
MAGEA4	NEFL	PSEN1	SMAD2	TLR10
MAGEB2	NFATC1	PSEN2	SMAD3	TLR2
MAGEC1	NFATC2	PSMB10	SMPD3	TLR3
MAGEC2	NFATC3	PSMB7	SOCS1	TLR4
MAP2K1	NFATC4	PSMB8	SPA17	TLR5
MAP2K2	NFKB1	PSMB9	SPACA3	TLR6
MAP2K4	NFKB2	PSMD7	SPINK5	TLR7
MAP3K1	NFKBIA	PTGS2	SPN	TLR8
MAP3K5	NLRC5	PTPRC	SPO11	TLR9
MAP3K7	NLRP3	PVR	SPP1	TMEFF2
MAP4K2	NOD1	PYCARD	SSX1	TNF
MAPK1	NOD2	RAG1	SSX4	TNFAIP3
MAPK11	NOTCH1	REL	ST6GAL1	TNFRSF10B
MAPK14	NRP1	RELA	STAT1	TNFRSF10C
TNFRSF11A	TNFRSF4	TNFSF18	TREM1	VEGFA
TNFRSF11B	TNFRSF8	TNFSF4	TREM2	VEGFC
TNFRSF12A	TNFRSF9	TNFSF8	TTK	XCL2
TNFRSF13B	TNFSF10	TOLLIP	TXK	XCR1
TNFRSF13C	TNFSF11	TP53	TXNIP	YTHDF2
TNFRSF14	TNFSF12	TPSAB1	TYK2	ZAP70
TNFRSF17	TNFSF13	TPTE	UBC	ZNF205
TNFRSF18	TNFSF13B	TRAF2	ULBP2
TNFRSF1A	TNFSF14	TRAF3	USP9Y
TNFRSF1B	TNFSF15	TRAF6	VCAM1

Bioinformatics analysis of RNA-Seq data

We used the upper quartile-normalized RNA-seq data by expectation maximization (RSEM) available for all selected cases at the TCGA data portal. No additional normalization was performed and the expression data were log2 transformed. All downstream data analysis was performed in MATLAB (Mathworks, Inc., Natick, Massachusetts, USA). Focusing on the first set of 122 samples, where clustering information is known, we performed separate analyses of the genes in each of the four immune panels (NanoString™, IFN-α, IFN-γ, and T cell cytotoxicity associated genes). Using a feature selection algorithm, genes were ranked based on their ability to discriminate samples belonging to one cluster from the remaining. The feature selection algorithm uses an ensemble of five different machine-learning techniques (unpublished). The analysis resulted in 16 ranking tables, four tables for each immune panel, where each table ranked genes on their ability to discriminate samples in one cluster from the rest.

In the hierarchical clustering analysis, the top 20% of genes in each ranking group were assembled to represent each immune panel, resulting in 377 genes (Nanostring), 44 genes (IFN-α), 91 genes (IFN-γ), and 62 genes (T cell cytotoxicity). A final feature selection ranking was performed where the combined unique set of 924 genes was used. The top 5% of genes in each of the four ranking groups were then merged, resulting in 157 unique genes.

Gene ontology analysis using PANTHER

We used the Protein Analysis Through Evolutionary Relationships (PANTHER), version 11.0, classification system (http://www.pantherdb.org/) [27] to determine dominant and enriched pathways in the top ranking 377 genes (NanoString panel) that were ranked based on their ability to discriminate samples across the four clusters. We then applied the statistical binomial overrepresentation test, as previously described in PANTHER [27], to derive the most dominant enriched pathways and gene ontology (GO) biological processes in our lists compared to the reference human genome. We performed these tests using both the 377 top ranking NanoString™ genes and 157 top ranked genes from all immune panels combined. The p-values were corrected for multiple testing using Bonferroni correction.

Validation of immune gene signature

The remaining 298 cases, not included in cohort 1, were treated as a validation group. From this cohort, patients with previous BCG therapy were excluded, leaving 225 cases. In order to assign samples in this set to each of the four clusters, only the top ranked 377 genes from the NanoString panel were used. First, for each cluster, two cluster centroids were computed using the expression data from cohort 1 (n=122; total of 8 cluster centroids). The cluster centroids were computed by taking the mean expression of samples in a given cluster (main cluster) and the mean expression of samples that do not belong to that cluster (alternative cluster). Then, for each sample in cohort 2, the Euclidean distance was computed to each of the 8 cluster centroids. A sample was assigned to a cluster with the smallest distance to the main cluster, but only if the distance to the main cluster was smaller than the distance to the alternative cluster. Alternatively, the sample was assigned to the ‘unclassified’ cluster. Using the newly assigned clustering information and ranked list of genes, unsupervised hierarchical clustering was performed.