Identification of a novel immune microenvironment signature predicting survival and therapeutic options for bladder cancer

Results

Identification of differentially expressed IRGs in BLCA tissues

Data on 403 UBC cases were downloaded from TCGA (TCGA-BLCA dataset) and randomly split into a training group (n=202) and the testing group (n=201). Another independent UBC dataset (GSE13507) downloaded from GEO constituted an independent testing set (n=165). Detailed demographics and clinicopathologic characteristics of patients are listed in Table 1. After applying cutoff thresholds of |log2FC|>1 and FDR<0.05, 1,924 differentially expressed genes between normal urothelium tissues and bladder tumors from TCGA and GTEx were identified using the Wilcoxon signed-rank test. Of these, 870 were upregulated and 1,054 downregulated (Figure 1A, 1B). 158 differentially expressed IRGs (DEIRGS) were extracted from these gene lists (Figure 1C).

Figure 1. Analysis of differentially expressed immune-related genes. (A) Heatmap and (B) volcano plot of differentially expressed genes (DEGs) in UBCs compared with normal healthy bladder samples in datasets from TCGA and GTEx. (C) Venn diagram of the 158 intersecting genes between DEGs and IRGs in UBCs.

Table 1. Clinical data in the training and testing sets.

Variables	Group	TCGA-BLCA			GSE13507 Testing cohort (n = 152)
		Total cohort (n = 403)	Training cohort (n = 202)	Testing cohort (n =201)
OS (months)
	Mean	25.5	24.0	27.1	48.38
	Range	0.4-168.3	0.4-168.3	0.6-168.0	1.03-136.97
Vital status
	Living	248	127	121	96
	Dead	155	75	80	69
T stage
	Ta/1	4	2	2	104
	T2	117	59	58	31
	T3	192	95	97	19
	T4	57	28	29	11
	Unknown	33	18	15
N stage
	Negative	235	123	112	149
	Positive	126	61	65	15
	Unknown	42	18	24	1
M stage
	M0	195	99	96	158
	M1	11	5	6	7
	Unknown	197	98	99	0
Grade
	Low	20	10	10	105
	High	380	190	190	60
	Unknown	3	2	1	0
Age (years)
	Mean	67.9	68.1	67.7	65.2
	Range	34-89	37-88	34-89	24-88
Gender
	Male	298	146	152	135
	Female	105	56	49	30

Functional enrichment analyses of DEIRGS

To explore the potential function of the DEIRGS, a functional enrichment analysis was performed. Gene Ontology (GO) terms enrichment analysis was done and the 10 most enriched terms for biological process (BP), cellular component (CC), and molecular function (MF) identified (Figure 2A). Among them, “positive regulation of cell migration”, “regulation of chemotaxis” and “cell proliferation” were considered related to cancer progression. Kyoto Encyclopedia of Gene and Genomes (KEGG) analysis identified the 7 pathways enriched for in the DEIRGS, such as “IL-17 signaling pathways”, “cytokine-cytokine receptor interaction”, and “MAPK signaling pathway” (Figure 2B). To further illustrate the relationships between the enriched terms, terms of 158 DEIRGs were extracted and presented as a network plot (Figure 2C). Protein-protein interaction enrichment analysis identified 6 closely connected network components (Figure 2D).

Figure 2. Enrichment and protein-protein interaction analysis of differentially expressed IRGs. (A) GO analysis of the immune-related genes. (B) Circular plot of the top seven most significant KEGG pathways. (C) Network plot of enriched terms: colored by cluster ID, where terms with a similarity > 0.3 are connected by edges. (D) Protein-protein interaction network and MCODE components identified in the gene lists.

Construction of a 10-gene risk signature using the UBC training group

To explore the prognostic value of the above 158 genes, we performed univariate Cox regression analysis in the UBC training group and identified 14 DEGs as significantly associated with UBC survival (p<0.05). Analysis of the 14 DEIRGs using LASSO to minimize overfitting (Supplementary Figure 2A, 2B) identified 10 IRGs as risk genes (Table 2). The 10 genes were then used to construct an immune-related model using the following formula: risk score = 0.2427 * AHNAK + (0.3584 * CALR) + (-0.0863 * PLA2G2A) + (-0.1880 * OAS1) + (-0.0401 * CDH1) + (0.2079 * PDGFRA) +

Table 2. Coefficients and LASSO Cox model results of each gene in 10-IRG risk signature.

Gene	HR	95% CI	p value	LASSO regression coefficient
AHNAK	1.372	1.083-1.739	0.008	0.243
CALR	1.82	1.169-2.832	0.008	0.358
PLA2G2A	0.902	0.816-0.996	0.043	-0.086
OAS1	0.726	0.583-0.904	0.004	-0.188
CDH1	0.851	0.730-0.991	0.038	-0.040
PDGFRA	1.348	1.087-1.672	0.006	0.208
SEMA3F	0.811	0.671-0.981	0.031	-0.001
RAC3	1.336	1.077-1.657	0.008	0.228
IRF5	0.702	0.545-0.905	0.006	-0.070
CARD11	0.745	0.632-0.879	0.0005	-0.119
CI, confidence interval; HR, hazard ratio.

(-0.0014 * SEMA3F) + (0.2281 * RAC3) + (-0.0702 * IRF5) + (-0.1192 * CARD11).

Next, the risk score in each UBC case was calculated and ranked based on its expression of the 10 risk genes. The UBC cases in the training cohort were divided into a high-risk group (n =101) and a low-risk group (n=101) based on the median risk score (Figure 3A). The survival status of the UBC cases was visualized on dot plots (Figure 3B). Kaplan-Meier analysis revealed significantly different survival rates between the two groups (p<0.05; Figure 3C). The AUC value for the risk model was 0.780 at 3 years and 0.754 at 5years for overall survival (OS; Figure 3D). Expression levels of the risk genes in the 2 groups were visualized by heatmaps (Figure 3E).

Figure 3. Identification of the 10-IRG risk model in the training set. (A, B) Risk score distribution and survival status of low-risk (green) (n=101) and high risk (red) UBC cases (n=101) in the training group. (C) Survival analysis of the high-risk UBC (n=101) and low-risk UBC (n=101) cases. (D) Time-dependent ROC curve analysis showing the 3-year (red) and 5-year (blue) survival of the 10 IRGs expression patterns. (E) Heat map showing expression level of risk genes in the high and low-risk BLCA patients in the training set.

Internal and external validation of the 10-gene risk model

Next, we tested the robustness of the 10-gene signature on 2 independent datasets. First, we repeated the analysis as described above in the internal TCGA testing group and observed lower survival rates for high-risk UBC cases relative to low-risk ones (p<0.05; Figure 4A). ROC curve analysis of overall survival revealed that the value of the AUC was 0.667 for 3 years and 0.679 for 5 years (Figure 4B). Risk genes expression levels in the internal testing group were visualized on heatmap (Figure 4C). Next, we validated the 10-gene risk model in the external GEO testing dataset. The cases were first grouped into a low-risk group (n=82) and a high-risk group (n=83) based on median risk scores. Interestingly, cases with high-risk had shorter median survival relative to low-risk ones (p<0.05; Figure 4D). In the GEO dataset, the AUC values were 0.652 at 3 years and 0.637 at 5 years (Figure 4E). The expression of risk genes in the GEO testing cohort were visualized in heatmap (Figure 4F). These data indicate that the immune-related risk model accurately predicts UBC prognosis.

Figure 4. Validation of the immune-based risk model in the internal TCGA and external GEO testing cohorts. (A, D) Kaplan-Meier curve analysis of overall survival of high-risk UBC (n=101 for TCGA; n=83 for GSE13507) and low-risk UBC (n=100 for TCGA; n=82 for GSE13507) cases. (B, E) Time-dependent ROC curve analysis showing the 3-year and 5-year survival of the10 IRGs signature. (C, F) Heatmaps showing expression of the selected genes in the immune-based risk model.

Relationship between the 10-gene risk model and clinicopathological features

To assess the clinical value of the immune gene-related signature, we evaluated the association between the 10-gene risk model and UBC clinicopathological features. This analysis showed that signature risk score positively correlates with grade, clinical stage, pathological T stage and metastatic lymphatic status (Figure 5A–5D). In the independent GEO cohort, the 10-IRG risk score was significantly elevated in high grade and high muscle-invasive UBC cases relative to low grade and superficial UBCs (Figure 5E–5G). No difference was observed between risk score and N stage (Figure 5H). These data suggested that our IRG risk model closely correlates with clinical UBC features.

Figure 5. Association between the 10-gene risk model and clinicopathological factors. (A–D) Correlation of the risk score and clinicopathological factors including (A) grade, (B) clinical stage, (C) T stage, (D) N stage in the TCGA-BLCA patient cohort. (E–H) Association between the risk score and clinicopathological features including (E) histological grade, (F) invasive status, (G) T stage and (H) N stage in the GEO UBC patient cohort.

Considering different clinical characteristics of bladder cancer patients, we performed a stratified analysis and subgroup analysis was stratified by age, gender, N stage, pathological stage, T stage and metastatic status (Figure 6A–6L). As shown in the Kaplan-Meier curves, we found that the survival rates of the high-risk groups had significantly decreased than the low-risk groups in any group. Thus, the predict power of overall survival by the 10-gene risk model was unrestricted to specific subgroups.

Figure 6. Risk-stratified analysis of the 10-gene signature for overall survival of bladder cancer patients. Kaplan-Meier curves for overall survival of patients in Age (≤65-year-old) subgroup (A), Age (>65-year-old) group (B), Female subgroup (C), Male subgroup (D), N0 subgroup (E), N1-N3 subgroup (F), clinical stage I-II subgroup (G), clinical stage III-IV subgroup (H), T1-2 subgroup (I), T3-4 subgroup (J), M0 subgroup (K), M1 subgroup (L).

Prognosis analysis and predictive accuracy of the 10-gene risk signature

We then analyzed the relationship between OS, clinicopathological parameters, including gender, age, clinical stage, pathological features (T and N status), and risk score of the immune-related risk model in the TCGA cohort. Univariate Cox regression analysis revealed that the above clinical variables correlate with UBC survival, except gender and risk score (p<0.05; Figure 7A). Multivariable Cox regression analysis showed that this risk model acts as an independent prognostic indicator in the TCGA (Figure 7B). Next, ROC analysis was done to evaluate the risk signature specificity and sensitivity in predicting OS relative to clinicopathological parameters. Interestingly, risk score of the 10 genes exhibited better AUC relative the clinicopathological parameters (Figure 7C), indicating that the risk model independently predicts UBC survival. In the c-index analysis (Table 3), the risk model showed better predictive ability than that of the TNM stage in the training, validation, and entire cohorts.

Figure 7. Prognostic significance and predictive accuracy of the immune-related risk model. (A, B) Univariate (A) and multivariate (B) Cox regression analyses of the risk scores of the clinical features for overall survival in TCGA-BLCA dataset. (C) ROC curve showing the specificity and sensitivity of the 10-gene signature risk score, age, gender clinical stage, T stage and lymph nodes status in predicting the OS of all UBC patients in the TCGA-BLCA dataset.

Table 3. Harrell's concordance indexes of the risk score and stage in different cohorts.

Cohort	Risk score	TNM Stage
TCGA-Training	0.75 (0.72-0.78)	0.67 (0.61-0.74)
TCGA-Validation	0.67 (0.63-0.71)	0.62 (0.58-0.66)
TCGA-Entire	0.73 (0.69-0.75)	0.65 (0.63-0.67)
GSE13507-Validation	0.66 (0.62-0.70)	0.59 (0.55-0.63)

Nomogram analysis of predictive accuracy

To predict UBC survival in the TCGA cohort, we constructed a nomogram by integrating risk scores of the immune-related signature and clinical information, including age, tumor stage, T stage and N stage. On the nomogram, the score on the point scale is easily obtained for each variable. Thus, survival can conveniently be estimated at 1, 3, 5 and 10 years by calculating the overall points (Figure 8A). The nomogram’s accuracy was evaluated using calibration curve analysis, which revealed consistency between the predictive curves at 1, 3, 5 and 10 years, and the ideal curve (Figure 8B–8E), indicating excellent performance by our nomogram. These data suggest that the nomogram is an effective tool for predicting UBC survival.

Figure 8. Construction and validation of a nomogram. (A) Nomogram for predicting the 1, 3, 5, and 10-year survival of UBC patients in the TCGA-BLCA cohort. (B–E) Calibration plot evaluating the consistency of the predicted value of the model and the probability of survival at 1, 3, 5, and 10 years obtained using the aforementioned nomogram.

Correlation of tumor-infiltrating immune cells (TIICs) with the 10-IRG risk model

TIICs have been suggested as independent predictors of survival in various cancer and are regulated by immune-related genes [27, 28]. We therefore evaluated the correlation between immune cells fractions in UBC tissues and risk models using CIBERSORT [29]. This analysis revealed that among the 22 immune cells, activated CD4+ memory T cell, macrophages and neutrophils positively correlated with the 10-IRG risk status, while CD8+ T cells, follicular helper T cells, regulatory T cells, and monocytes negatively correlate with the immune-related risk score (Figure 9A). To cross examine the data, further analysis showed that macrophages and neutrophils were positively correlated with the 10-IRG risk status, while CD8+ T cells were negatively correlate with the immune-related risk score with both methods for deconvolution (Figure 9A and Supplementary Figure 3). Moreover, a correlation heatmap demonstrated that different immune cells fractions have a weak/moderate correlation (Supplementary Figure 4). The proportion of CD8+ T cells, dendritic cells resting, and macrophages (Figure 9B–9D) significantly correlate with survival. These data may reveal the poor clinical outcome in the high-risk cohort in part.

Figure 9. Alteration of immune cell infiltration in BLCA samples with different risk status. (A) Violin plot demonstrating the TILCs associated with the risk model. High- and low-risk groups are represented by red and green violin, respectively. (B–D) Kaplan-Meier curve analysis of overall survival for various immune cells infiltration. (B) CD8 T cells; (C) dendritic cells resting; (D) macrophages M0.

Relationship of immune-related risk model with mutation profile, molecular subtype and immunotherapy markers

To assess if the immune-related risk signature was associated with specific tumor mutations, the correlation between the signature and mutations was analyzed in TCGA BLCA dataset containing somatic mutation profiles. The overall landscape of mutation profile was illustrated in Supplementary Figure 5. Alterations in the mutation landscape in high or low risk group were as follows: 6 genes were mutated in >19% of tissues with high risk score: TP53 (56%), TTN (39%), KMT2D (30%), ARID1A (28%), RB1 (25%) and KDM6A (20%). While eight genes were mutated in >19% of tissues with low risk score: TTN (42%), TP53 (38%), KDM6A (31%), MUC16 (30%), PIK3CA (22%), FGFR3 (22%), KMT2D (28%) and ARID1A (21%) (Figure 10A). Notably, higher rates of FGFR3 mutation occurred in bladder cancer tissues with low risk score compared with tumors with high risk score. However, high TP53 and RB1 mutation rates occurred in the high-risk subgroup. Next, we evaluated the relationship between the risk signature and UBC molecular subtypes. Interestingly, for luminal papillary UBC, strikingly more patients had low risk scores relative to those with high risk scores. For basal-squamous and neuronal UBC most patients had high risk scores (Figure 10B). Moreover, patients with high risk exhibited higher PD1 and PD-L1 expression (p<0.05; Figure 10C, 10D). These results suggest that high-risk patients in the 10-IRG signature are promising candidates for immune checkpoint inhibitors or other treatment management.

Figure 10. Variations in the mutation landscape and molecular features of the high and low-risk UBCs in TCGA-BLCA cohort. (A) Mutation landscapes of UBC samples in the low- and high-risk groups. (B) Association between the 10-genes risk model and molecular subtypes of UBC patients in the TCGA-BLCA dataset. (C, D) Normalized PD1 and PD-L1 gene expression profile in the low- and high-risk groups.

Identification of ten-gene risk score correlated biological pathways

In order to investigate biological function of high- and low risk score in bladder cancer, GSEA was performed and unsurprisingly, the results revealed that high-risk score was closely correlated with cancer related pathways including ADHERENS JUNCTION, FOCOL ADHESION, CELL CYCLE, DNA REPLICATION (Figure 11A). Moreover, patients with recurrence/ progression tended to have a higher risk score than those with disease free (Figure 11B). Unexpectedly, genes involved in glycolysis were significantly enriched in patients with a high-risk score (Figure 11C). Significant positive correlation between core glycolytic genes, including SLC2A3(GLUT3), PFKFB3, PKM2, LDHA and SLC16A1/4(MCT1/4), and risk score was observed (Figure 11D, 11E). These results indicate that this ten-gene risk score is closely related to tumor progression and glycolytic status of UBC. Interestingly, we observed a trend for higher expression of Ki-67 in high-risk groups versus the low-risk group (Figure 11F).

Figure 11. Enrichment analysis of cancer and glycolysis related pathways associated with risk score in TCGA-BLCA database. (A) Gene set enrichment analysis (GSEA) on the association of risk score with adherens junction, focal adhesion, cell cycle and DNA replication. (B) Risk score of UBC patients with disease free and progression/recurrence in TCGA-BLCA dataset. (C) GSEA plots depicting enrichment of glycolysis pathways and risk score. (D) The correlation between the expression of glycolytic genes and the risk score. (E) The identified core candidate genes (red) involved in enzymatic glycolysis are shown. (F) The representative IHC images of ki-67 expression in UBC tissues with low- or high-risk.

Regulatory network among differentially expressed TFs and IRGs

To elucidate the mechanisms involved in the regulation of these immune-related genes, we examined differentially expressed transcription factors and correlation analysis of DEGs between the 2 gene lists (Figure 12A). Construction of the regulatory network revealed KLF5, SOX4, GRHL2 and GATA3 as the nodes with highest connectivity (Figure 12B). Interestingly, further investigation via Tumor Immune Estimation Resource (TIMER) showed these transcriptional factors, particularly GATA3, were closely related to the immune cell infiltration of tumors (Figure 12C). These results indicate that the 4 genes may involve in the upstream regulation of the ten IRGs and influences the immune status of the tumor microenvironment.

Figure 12. Construction of the TF-IRG network in UBC group. (A) Diagram showing the construction of regulatory network of TFs and genes in the immune-related risk signature. (B) Regulatory network of the differentially expressed TFs and IRGs. (C) Correlation of GATA3, GRHL2, SOX4 and KLF5 expression with immune cell infiltration level by TIMER in UBCs.

Author Contributions

BS designed the experiments, supervised the progress throughout this study and revised the manuscript; YLY constructed the figures and also performed experiments; ZNH, JMC, PFT, FZ and MYT provided technical support. All authors reviewed and approved the final manuscript.

Acknowledgments

We are grateful to the contributors of data to public databases including TCGA, GTEx, GEO and UCSC Xena.

Conflicts of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflicts of interest.

Funding

This work was supported by the National Natural Science Foundation (81972370 to MT), Shanghai Pujiang Program of China (No.18PJD042 to MT), Wu Jieping Medical Foundation (320.6750.16051 to BS), Shanghai Songjiang Municipal Science and Technology Commission Natural Science Foundation (17SJKJGG10 to BS) and Shanghai Specialized Research Fund for Integrated Chinese and Western Medicine in General Hospitals (ZHYY-ZXYJHZX-1-201705 to BS). The founding sponsor had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the paper; and in the decision to publish the results.

References

• 1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin. 2020; 70:7–30. https://doi.org/10.3322/caac.21590 [PubMed]
• 2. Dy GW, Gore JL, Forouzanfar MH, Naghavi M, Fitzmaurice C. Global burden of urologic cancers, 1990-2013. Eur Urol. 2017; 71:437–46. https://doi.org/10.1016/j.eururo.2016.10.008 [PubMed]
• 3. Alfred Witjes J, Lebret T, Compérat EM, Cowan NC, De Santis M, Bruins HM, Hernández V, Espinós EL, Dunn J, Rouanne M, Neuzillet Y, Veskimäe E, van der Heijden AG, et al. Updated 2016 EAU guidelines on muscle-invasive and metastatic bladder cancer. Eur Urol. 2017; 71:462–75. https://doi.org/10.1016/j.eururo.2016.06.020 [PubMed]
• 4. Parikh RB, Adamson BJ, Khozin S, Galsky MD, Baxi SS, Cohen A, Mamtani R. Association between FDA label restriction and immunotherapy and chemotherapy use in bladder cancer. JAMA. 2019; 322:1209–11. https://doi.org/10.1001/jama.2019.10650 [PubMed]
• 5. Robertson AG, Kim J, Al-Ahmadie H, Bellmunt J, Guo G, Cherniack AD, Hinoue T, Laird PW, Hoadley KA, Akbani R, Castro MA, Gibb EA, Kanchi RS, et al, and TCGA Research Network. Comprehensive molecular characterization of muscle-invasive bladder cancer. Cell. 2018; 174:1033. https://doi.org/10.1016/j.cell.2018.07.036 [PubMed]
• 6. Kamoun A, de Reyniès A, Allory Y, Sjödahl G, Robertson AG, Seiler R, Hoadley KA, Groeneveld CS, Al-Ahmadie H, Choi W, Castro MA, Fontugne J, Eriksson P, et al, and Bladder Cancer Molecular Taxonomy Group. A consensus molecular classification of muscle-invasive bladder cancer. Eur Urol. 2020; 77:420–33. https://doi.org/10.1016/j.eururo.2019.09.006 [PubMed]
• 7. Wang L, Saci A, Szabo PM, Chasalow SD, Castillo-Martin M, Domingo-Domenech J, Siefker-Radtke A, Sharma P, Sfakianos JP, Gong Y, Dominguez-Andres A, Oh WK, Mulholland D, et al. EMT- and stroma-related gene expression and resistance to PD-1 blockade in urothelial cancer. Nat Commun. 2018; 9:3503. https://doi.org/10.1038/s41467-018-05992-x [PubMed]
• 8. Cheng S, Xu C, Jin Y, Li Y, Zhong C, Ma J, Yang J, Zhang N, Li Y, Wang C, Yang Z, Wang Y. Artificial mini dendritic cells boost T cell-based immunotherapy for ovarian cancer. Adv Sci (Weinh). 2020; 7:1903301. https://doi.org/10.1002/advs.201903301 [PubMed]
• 9. Devalaraja S, To TK, Folkert IW, Natesan R, Alam MZ, Li M, Tada Y, Budagyan K, Dang MT, Zhai L, Lobel GP, Ciotti GE, Eisinger-Mathason TS, et al. Tumor-derived retinoic acid regulates intratumoral monocyte differentiation to promote immune suppression. Cell. 2020; 180:1098–114.e16. https://doi.org/10.1016/j.cell.2020.02.042 [PubMed]
• 10. Shekarian T, Sivado E, Jallas AC, Depil S, Kielbassa J, Janoueix-Lerosey I, Hutter G, Goutagny N, Bergeron C, Viari A, Valsesia-Wittmann S, Caux C, Marabelle A. Repurposing rotavirus vaccines for intratumoral immunotherapy can overcome resistance to immune checkpoint blockade. Sci Transl Med. 2019; 11:eaat5025. https://doi.org/10.1126/scitranslmed.aat5025 [PubMed]
• 11. Schachter J, Ribas A, Long GV, Arance A, Grob JJ, Mortier L, Daud A, Carlino MS, McNeil C, Lotem M, Larkin J, Lorigan P, Neyns B, et al. Pembrolizumab versus ipilimumab for advanced melanoma: final overall survival results of a multicentre, randomised, open-label phase 3 study (KEYNOTE-006). Lancet. 2017; 390:1853–62. https://doi.org/10.1016/S0140-6736(17)31601-X [PubMed]
• 12. Herbst RS, Baas P, Kim DW, Felip E, Pérez-Gracia JL, Han JY, Molina J, Kim JH, Arvis CD, Ahn MJ, Majem M, Fidler MJ, de Castro G Jr, et al. Pembrolizumab versus docetaxel for previously treated, PD-L1-positive, advanced non-small-cell lung cancer (KEYNOTE-010): a randomised controlled trial. Lancet. 2016; 387:1540–50. https://doi.org/10.1016/S0140-6736(15)01281-7 [PubMed]
• 13. Sundahl N, Vandekerkhove G, Decaestecker K, Meireson A, De Visschere P, Fonteyne V, De Maeseneer D, Reynders D, Goetghebeur E, Van Dorpe J, Verbeke S, Annala M, Brochez L, et al. Randomized phase 1 trial of pembrolizumab with sequential versus concomitant stereotactic body radiotherapy in metastatic urothelial carcinoma. Eur Urol. 2019; 75:707–11. https://doi.org/10.1016/j.eururo.2019.01.009 [PubMed]
• 14. Gourd E. Neoadjuvant pembrolizumab in bladder cancer. Lancet Oncol. 2018; 19:e669. https://doi.org/10.1016/S1470-2045(18)30814-3 [PubMed]
• 15. Cascone T, McKenzie JA, Mbofung RM, Punt S, Wang Z, Xu C, Williams LJ, Wang Z, Bristow CA, Carugo A, Peoples MD, Li L, Karpinets T, et al. Increased tumor glycolysis characterizes immune resistance to adoptive T cell therapy. Cell Metab. 2018; 27:977–87.e4. https://doi.org/10.1016/j.cmet.2018.02.024 [PubMed]
• 16. Wan W, Peng K, Li M, Qin L, Tong Z, Yan J, Shen B, Yu C. Histone demethylase JMJD1A promotes urinary bladder cancer progression by enhancing glycolysis through coactivation of hypoxia inducible factor 1α. Oncogene. 2017; 36:3868–77. https://doi.org/10.1038/onc.2017.13 [PubMed]
• 17. Ohashi T, Akazawa T, Aoki M, Kuze B, Mizuta K, Ito Y, Inoue N. Dichloroacetate improves immune dysfunction caused by tumor-secreted lactic acid and increases antitumor immunoreactivity. Int J Cancer. 2013; 133:1107–18. https://doi.org/10.1002/ijc.28114 [PubMed]
• 18. Yang Y, Sui Y, Xie B, Qu H, Fang X. GliomaDB: a web server for integrating glioma omics data and interactive analysis. Genomics Proteomics Bioinformatics. 2019; 17:465–71. https://doi.org/10.1016/j.gpb.2018.03.008 [PubMed]
• 19. Sun D, Xi Y, Rodriguez B, Park HJ, Tong P, Meong M, Goodell MA, Li W. MOABS: model based analysis of bisulfite sequencing data. Genome Biol. 2014; 15:R38. https://doi.org/10.1186/gb-2014-15-2-r38 [PubMed]
• 20. Wang Y, Du L, Yang X, Li J, Li P, Zhao Y, Duan W, Chen Y, Wang Y, Mao H, Wang C. A nomogram combining long non-coding RNA expression profiles and clinical factors predicts survival in patients with bladder cancer. Aging (Albany NY). 2020; 12:2857–79. https://doi.org/10.18632/aging.102782 [PubMed]
• 21. Elamin AA, Klunkelfuß S, Kämpfer S, Oehlmann W, Stehr M, Smith C, Simpson GR, Morgan R, Pandha H, Singh M. A specific blood signature reveals higher levels of S100A12: a potential bladder cancer diagnostic biomarker along with urinary engrailed-2 protein detection. Front Oncol. 2020; 9:1484. https://doi.org/10.3389/fonc.2019.01484 [PubMed]
• 22. Wu J, Zhao Y, Zhang J, Wu Q, Wang W. Development and validation of an immune-related gene pairs signature in colorectal cancer. Oncoimmunology. 2019; 8:1596715. https://doi.org/10.1080/2162402X.2019.1596715 [PubMed]
• 23. Huang R, Mao M, Lu Y, Yu Q, Liao L. A novel immune-related genes prognosis biomarker for melanoma: associated with tumor microenvironment. Aging (Albany NY). 2020; 12:6966–80. https://doi.org/10.18632/aging.103054 [PubMed]
• 24. Wan B, Liu B, Yu G, Huang Y, Lv C. Differentially expressed autophagy-related genes are potential prognostic and diagnostic biomarkers in clear-cell renal cell carcinoma. Aging (Albany NY). 2019; 11:9025–42. https://doi.org/10.18632/aging.102368 [PubMed]
• 25. Lian P, Wang Q, Zhao Y, Chen C, Sun X, Li H, Deng J, Zhang H, Ji Z, Zhang X, Huang Q. An eight-long non-coding RNA signature as a candidate prognostic biomarker for bladder cancer. Aging (Albany NY). 2019; 11:6930–40. https://doi.org/10.18632/aging.102225 [PubMed]
• 26. Lin P, Wen DY, Chen L, Li X, Li SH, Yan HB, He RQ, Chen G, He Y, Yang H. A radiogenomics signature for predicting the clinical outcome of bladder urothelial carcinoma. Eur Radiol. 2020; 30:547–57. https://doi.org/10.1007/s00330-019-06371-w [PubMed]
• 27. Kim K, Park S, Park SY, Kim G, Park SM, Cho JW, Kim DH, Park YM, Koh YW, Kim HR, Ha SJ, Lee I. Single-cell transcriptome analysis reveals TOX as a promoting factor for T cell exhaustion and a predictor for anti-PD-1 responses in human cancer. Genome Med. 2020; 12:22. https://doi.org/10.1186/s13073-020-00722-9 [PubMed]
• 28. Sato J, Kitano S, Motoi N, Ino Y, Yamamoto N, Watanabe S, Ohe Y, Hiraoka N. CD20⁺ tumor-infiltrating immune cells and CD204⁺ M2 macrophages are associated with prognosis in thymic carcinoma. Cancer Sci. 2020; 111:1921–32. https://doi.org/10.1111/cas.14409 [PubMed]
• 29. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–57. https://doi.org/10.1038/nmeth.3337 [PubMed]
• 30. Gómez de Liaño Lista A, van Dijk N, de Velasco Oria de Rueda G, Necchi A, Lavaud P, Morales-Barrera R, Alonso Gordoa T, Maroto P, Ravaud A, Durán I, Szabados B, Castellano D, Giannatempo P, et al. Clinical outcome after progressing to frontline and second-line anti-PD-1/PD-L1 in advanced urothelial cancer. Eur Urol. 2020; 77:269–76. https://doi.org/10.1016/j.eururo.2019.10.004 [PubMed]
• 31. Petrylak DP, Powles T, Bellmunt J, Braiteh F, Loriot Y, Morales-Barrera R, Burris HA, Kim JW, Ding B, Kaiser C, Fassò M, O’Hear C, Vogelzang NJ. Atezolizumab (MPDL3280A) monotherapy for patients with metastatic urothelial cancer: long-term outcomes from a phase 1 study. JAMA Oncol. 2018; 4:537–44. https://doi.org/10.1001/jamaoncol.2017.5440 [PubMed]
• 32. Matza D, Badou A, Kobayashi KS, Goldsmith-Pestana K, Masuda Y, Komuro A, McMahon-Pratt D, Marchesi VT, Flavell RA. A scaffold protein, AHNAK1, is required for calcium signaling during T cell activation. Immunity. 2008; 28:64–74. https://doi.org/10.1016/j.immuni.2007.11.020 [PubMed]
• 33. Schcolnik-Cabrera A, Oldak B, Juárez M, Cruz-Rivera M, Flisser A, Mendlovic F. Calreticulin in phagocytosis and cancer: opposite roles in immune response outcomes. Apoptosis. 2019; 24:245–55. https://doi.org/10.1007/s10495-019-01532-0 [PubMed]
• 34. Okita Y, Shiono T, Yahagi A, Hamada S, Umemura M, Matsuzaki G. Interleukin-22-induced antimicrobial phospholipase A2 group IIA mediates protective innate immunity of nonhematopoietic cells against listeria monocytogenes. Infect Immun. 2015; 84:573–79. https://doi.org/10.1128/IAI.01000-15 [PubMed]
• 35. Schwartz SL, Park EN, Vachon VK, Danzy S, Lowen AC, Conn GL. Human OAS1 activation is highly dependent on both RNA sequence and context of activating RNA motifs. Nucleic Acids Res. 2020; 48:7520–31. https://doi.org/10.1093/nar/gkaa513 [PubMed]
• 36. An Y, Adams JR, Hollern DP, Zhao A, Chang SG, Gams MS, Chung PE, He X, Jangra R, Shah JS, Yang J, Beck LA, Raghuram N, et al. Cdh1 and Pik3ca mutations cooperate to induce immune-related invasive lobular carcinoma of the breast. Cell Rep. 2018; 25:702–14.e6. https://doi.org/10.1016/j.celrep.2018.09.056 [PubMed]
• 37. Pantaleo MA, Tarantino G, Agostinelli C, Urbini M, Nannini M, Saponara M, Castelli C, Stacchiotti S, Fumagalli E, Gatto L, Santini D, De Leo A, Marafioti T, et al. Immune microenvironment profiling of gastrointestinal stromal tumors (GIST) shows gene expression patterns associated to immune checkpoint inhibitors response. Oncoimmunology. 2019; 8:e1617588. https://doi.org/10.1080/2162402X.2019.1617588 [PubMed]
• 38. Alfaro D, García-Ceca JJ, Cejalvo T, Jiménez E, Jenkinson EJ, Anderson G, Muñoz JJ, Zapata A. EphrinB1-EphB signaling regulates thymocyte-epithelium interactions involved in functional T cell development. Eur J Immunol. 2007; 37:2596–605. https://doi.org/10.1002/eji.200737097 [PubMed]
• 39. Hu M, Lu Y, Qi Y, Zhang Z, Wang S, Xu Y, Chen F, Tang Y, Chen S, Chen M, Du C, Shen M, Wang F, et al. SRC-3 functions as a coactivator of T-bet by regulating the maturation and antitumor activity of natural killer cells. Cancer Immunol Res. 2020; 8:1150–62. https://doi.org/10.1158/2326-6066.CIR-20-0181 [PubMed]
• 40. Arnold IC, Artola-Boran M, Gurtner A, Bertram K, Bauer M, Frangez Z, Becher B, Kopf M, Yousefi S, Simon HU, Tzankov A, Müller A. The GM-CSF-IRF5 signaling axis in eosinophils promotes antitumor immunity through activation of type 1 T cell responses. J Exp Med. 2020; 217:e20190706. https://doi.org/10.1084/jem.20190706 [PubMed]
• 41. Di Pilato M, Kim EY, Cadilha BL, Prüßmann JN, Nasrallah MN, Seruggia D, Usmani SM, Misale S, Zappulli V, Carrizosa E, Mani V, Ligorio M, Warner RD, et al. Targeting the CBM complex causes T_reg cells to prime tumours for immune checkpoint therapy. Nature. 2019; 570:112–16. https://doi.org/10.1038/s41586-019-1215-2 [PubMed]
• 42. Lee H, Kim K, Woo J, Park J, Kim H, Lee KE, Kim H, Kim Y, Moon KC, Kim JY, Park IA, Shim BB, Moon JH, et al. Quantitative proteomic analysis identifies AHNAK (neuroblast differentiation-associated protein AHNAK) as a novel candidate biomarker for bladder urothelial carcinoma diagnosis by liquid-based cytology. Mol Cell Proteomics. 2018; 17:1788–802. https://doi.org/10.1074/mcp.RA118.000562 [PubMed]
• 43. Wu Q, Fu C, Li M, Li J, Li Z, Qi L, Ci X, Ma G, Gao A, Fu X, A J, An N, Liu M, et al. CINP is a novel cofactor of KLF5 required for its role in the promotion of cell proliferation, survival and tumor growth. Int J Cancer. 2019; 144:582–94. https://doi.org/10.1002/ijc.31908 [PubMed]
• 44. Al-Ahmadie HA, Iyer G, Lee BH, Scott SN, Mehra R, Bagrodia A, Jordan EJ, Gao SP, Ramirez R, Cha EK, Desai NB, Zabor EC, Ostrovnaya I, et al. Frequent somatic CDH1 loss-of-function mutations in plasmacytoid variant bladder cancer. Nat Genet. 2016; 48:356–58. https://doi.org/10.1038/ng.3503 [PubMed]
• 45. Korpal M, Puyang X, Jeremy Wu Z, Seiler R, Furman C, Oo HZ, Seiler M, Irwin S, Subramanian V, Julie Joshi J, Wang CK, Rimkunas V, Tortora D, et al. Evasion of immunosurveillance by genomic alterations of PPARγ/RXRα in bladder cancer. Nat Commun. 2017; 8:103. https://doi.org/10.1038/s41467-017-00147-w [PubMed]
• 46. Chevalier MF, Trabanelli S, Racle J, Salomé B, Cesson V, Gharbi D, Bohner P, Domingos-Pereira S, Dartiguenave F, Fritschi AS, Speiser DE, Rentsch CA, Gfeller D, et al. ILC2-modulated T cell-to-MDSC balance is associated with bladder cancer recurrence. J Clin Invest. 2017; 127:2916–29. https://doi.org/10.1172/JCI89717 [PubMed]
• 47. Baras AS, Drake C, Liu JJ, Gandhi N, Kates M, Hoque MO, Meeker A, Hahn N, Taube JM, Schoenberg MP, Netto G, Bivalacqua TJ. The ratio of CD8 to treg tumor-infiltrating lymphocytes is associated with response to cisplatin-based neoadjuvant chemotherapy in patients with muscle invasive urothelial carcinoma of the bladder. Oncoimmunology. 2016; 5:e1134412. https://doi.org/10.1080/2162402X.2015.1134412 [PubMed]
• 48. Kobatake K, Ikeda KI, Nakata Y, Yamasaki N, Ueda T, Kanai A, Sentani K, Sera Y, Hayashi T, Koizumi M, Miyakawa Y, Inaba T, Sotomaru Y, et al. Kdm6a deficiency activates inflammatory pathways, promotes M2 macrophage polarization, and causes bladder cancer in cooperation with p53 dysfunction. Clin Cancer Res. 2020; 26:2065–79. https://doi.org/10.1158/1078-0432.CCR-19-2230 [PubMed]
• 49. Dorschner RA, Lee J, Cohen O, Costantini T, Baird A, Eliceiri BP. ECRG4 regulates neutrophil recruitment and CD44 expression during the inflammatory response to injury. Sci Adv. 2020; 6:eaay0518. https://doi.org/10.1126/sciadv.aay0518 [PubMed]
• 50. Ho PC, Bihuniak JD, Macintyre AN, Staron M, Liu X, Amezquita R, Tsui YC, Cui G, Micevic G, Perales JC, Kleinstein SH, Abel ED, Insogna KL, et al. Phosphoenolpyruvate is a metabolic checkpoint of anti-tumor T cell responses. Cell. 2015; 162:1217–28. https://doi.org/10.1016/j.cell.2015.08.012 [PubMed]
• 51. Yang Z, He L, Lin K, Zhang Y, Deng A, Liang Y, Li C, Wen T. The KMT1A-GATA3-STAT3 circuit is a novel self-renewal signaling of human bladder cancer stem cells. Clin Cancer Res. 2017; 23:6673–85. https://doi.org/10.1158/1078-0432.CCR-17-0882 [PubMed]
• 52. Ito K, Murphy D. Application of ggplot2 to pharmacometric graphics. CPT Pharmacometrics Syst Pharmacol. 2013; 2:e79. https://doi.org/10.1038/psp.2013.56 [PubMed]
• 53. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019; 10:1523. https://doi.org/10.1038/s41467-019-09234-6 [PubMed]
• 54. Bader GD, Hogue CW. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. 2003; 4:2. https://doi.org/10.1186/1471-2105-4-2 [PubMed]
• 55. Wang H, Lengerich BJ, Aragam B, Xing EP. Precision lasso: accounting for correlations and linear dependencies in high-dimensional genomic data. Bioinformatics. 2019; 35:1181–87. https://doi.org/10.1093/bioinformatics/bty750 [PubMed]
• 56. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 2017; 77:e108–10. https://doi.org/10.1158/0008-5472.CAN-17-0307 [PubMed]
• 57. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018; 28:1747–56. https://doi.org/10.1101/gr.239244.118 [PubMed]