Research Paper Volume 13, Issue 2 pp 2780—2802

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

Yilin Yan1, , Zhengnan Huang1, , Jinming Cai1, , Pengfei Tang1, , Fang Zhang1, , Mingyue Tan1,2, , Bing Shen1, ,

  • 1 Department of Urology, Shanghai General Hospital, Shanghai Jiaotong University School of Medicine, Shanghai 200080, China
  • 2 Department of Urology, Shuguang Hospital, Shanghai University of Traditional Chinese Medicine, Shanghai 200021, China

Received: May 19, 2020       Accepted: November 15, 2020       Published: December 19, 2020      

https://doi.org/10.18632/aging.202327
How to Cite

Copyright: © 2020 Yan et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Few studies have investigated the potential of tumor immune microenvironment genes as indicators of urinary bladder cancer. Here, we sought to establish an immune-related gene signature for determining prognosis and treatment options. We developed a ten-gene tumor immune microenvironment signature and evaluated its prognostic capacity on internal and external cohorts. Multivariate Cox regression and nomogram analyses revealed the prognostic risk model as an independent and effective indicator of prognosis. We observed lower proportions of CD8+ T cells, dendritic cells, regulatory T cells, higher proportions of macrophages and neutrophils in high UBC risk group. UBC tissues with high-risk score tend to exhibit high TP53 and RB1 mutation rates, high PD1/PD-L1 expression and poor-survival basal squamous subtypes, while those with low-risk score tend to have high FGFR3 mutation rates and luminal papillary subtypes. Unexpectedly, we found a highly significant positive correlation between glycolytic genes and risk score, highlighting metabolic competition in tumor ecosystem and potential therapeutic avenues. Our study thus revealed a tumor immune microenvironment signature for predicting prognostic and response to immune checkpoint inhibitors against bladder cancer. Prospective studies are required to further test the predictive capacity of this model.

Introduction

Globally, urinary bladder cancer (UBC) is one of the most common malignancies of the urinary tract, with an estimated 550,000 new cases and 200,000 deaths [1]. UBC is generally subdivided into 2 classes, non-muscle invasive bladder cancer (NMIBC), which accounts for about 75% of UBC cases, or muscle-invasive bladder cancer (MIBC) that makes up about 25% of the cases. NMIBC is often treatable but is associated with frequent relapse. Up to 15% of NMIBC cases progress into MIBC, which is more aggressive [2]. Despite recent advances in cancer diagnosis and treatment, UBC is still associated with a high rate of metastasis [3]. Thus, there is an urgent need to better understand the molecular basis of UBC development and progression so as to improve outcomes.

Currently the standard treatment of MIBC is cisplatin-based chemotherapy [4]. Recent studies have demonstrated that UBC is characterized by high heterogeneity and genomic instability [5, 6]. UBC resulting from various driver events responds differently to treatment. For example, tumors characterized by a mesenchymal related signature appear more sensitive to immune checkpoint inhibitors and resistant to cisplatin-based chemotherapy [7].

Cancer immunotherapy has shown promise against various cancers [810]. Immune checkpoint inhibitor-based immunotherapy targeting programmed cell with cytotoxic T lymphocyte antigen 4 (CTLA4), programmed cell death 1 (PD1) and death-ligand 1 (PD-L1) have shown significant efficacy against advanced cancer [11, 12]. While UBC has been reported to respond to immune checkpoint inhibitors, the response to the outcomes have been poor [13]. Despite low efficacy, pembrolizumab (anti-PD1) is approved for treating metastatic UBC [14]. The tumor immune microenvironment (TME) is thought to influence UBC progression and understanding its characteristics may unveil biomarkers for immunotherapy effectiveness and precision treatment.

However, from the immunotherapy perspective, recent studies have discussed the role of tumor glycolysis in therapeutic resistance [15]. Enhancing glycolysis has been reported to be utilized and manipulated by bladder tumor cells to promote tumor progression [16]. The glycolytic switch of cancer cell, namely “aerobic glycolysis” or “Warburg effect”, not only promotes tumor growth and metastasis but also contributes to the acidic microenvironment by releasing lactate into the extracellular environment, limiting the development of an effective antitumor immune response. Sporadic reports have demonstrated that inhibition of tumor glycolysis could promote antitumor immune response [17]. Thus, glycolysis blockade in cancer cells could provide a novel approach to overcoming immunotherapy resistance.

Data-mining the deposited information of multi-omics data using bioinformatic methods have become essential tools for a deeper understanding of tumor biology [18, 19]. Transcriptome profiling has uncovered various gene signatures that may predict UBC risk [20, 21]. Multiple immune-related genes (IRGs) risk models have been proposed for diagnosis and prognosis of various cancers, including colorectal cancer [22], melanoma [23] and clear-cell renal cell cancer [24]. However, most studies have focused on diagnostic gene signatures for UBC [25, 26]. Therefore, there is an urgent need to develop a modified prognostic risk model of UBC to predict outcome and treatment options.

Herein, we evaluated publicly available UBC datasets, aiming to identify an IRG signature associated UBC prognosis (Supplementary Figure 1). We identified 158 differentially expressed IRGs in UBC samples and constructed an IRG risk model using Cox and Lasso regression analyses. We then validated the IRG signature by evaluating its prognostic accuracy and the relationship of the signature with clinicopathological features. Changes in tumor immune cell infiltration, mutation profile, molecular subtypes and glycolytic biological pathway associated with the risk model of UBC were also explored. Finally, potential upstream regulatory transcription factors for the ten genes in the signature were analyzed. Our study may help provide a deeper understanding and more specific personalized therapies for UBCs.

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).

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.

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.

VariablesGroupTCGA-BLCAGSE13507 Testing cohort (n = 152)
Total cohort (n = 403)Training cohort (n = 202)Testing cohort (n =201)
OS (months)
Mean25.524.027.148.38
Range0.4-168.30.4-168.30.6-168.01.03-136.97
Vital status
Living24812712196
Dead155758069
T stage
Ta/1422104
T2117595831
T3192959719
T457282911
Unknown331815
N stage
Negative235123112149
Positive126616515
Unknown4218241
M stage
M01959996158
M111567
Unknown19798990
Grade
Low201010105
High38019019060
Unknown3210
Age (years)
Mean67.968.167.765.2
Range34-8937-8834-8924-88
Gender
Male298146152135
Female105564930

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).

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.

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.

GeneHR95% CIp valueLASSO regression coefficient
AHNAK1.3721.083-1.7390.0080.243
CALR1.821.169-2.8320.0080.358
PLA2G2A0.9020.816-0.9960.043-0.086
OAS10.7260.583-0.9040.004-0.188
CDH10.8510.730-0.9910.038-0.040
PDGFRA1.3481.087-1.6720.0060.208
SEMA3F0.8110.671-0.9810.031-0.001
RAC31.3361.077-1.6570.0080.228
IRF50.7020.545-0.9050.006-0.070
CARD110.7450.632-0.8790.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).

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.

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.

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.

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 5A5D). 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 5E5G). 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.

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.

Figure 5. Association between the 10-gene risk model and clinicopathological factors. (AD) 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. (EH) 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 6A6L). 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.

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).

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.

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.

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.

CohortRisk scoreTNM Stage
TCGA-Training0.75 (0.72-0.78)0.67 (0.61-0.74)
TCGA-Validation0.67 (0.63-0.71)0.62 (0.58-0.66)
TCGA-Entire0.73 (0.69-0.75)0.65 (0.63-0.67)
GSE13507-Validation0.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 8B8E), indicating excellent performance by our nomogram. These data suggest that the nomogram is an effective tool for predicting UBC survival.

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.

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. (BE) 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 9B9D) significantly correlate with survival. These data may reveal the poor clinical outcome in the high-risk cohort in part.

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.

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. (BD) 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.

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.

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).

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.

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.

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.

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.

Discussion

Mounting evidence has shown that not all patients benefit from chemotherapy [30]. However, despite growing interest in immunotherapy, only a small proportion UBC patients respond to immunotherapy [31]. Therefore, personalized therapies are proposed as means to optimize outcomes. The complexity of the response to immunotherapy makes it impossible that a single gene could be sufficient to predict prognosis and response to immunotherapy. Recent studies have reported multiple gene expression signatures, including immune-related signature, that can predict cancer outcomes patients with UBC [21, 25]. However, very few of these are prognostic or capable of guiding treatment. Here, we have identified an immune-related risk model for UBC by analyzing TCGA and GTEx datasets. Using the 2 datasets minimizes bias from having too few normal samples in the TCGA dataset.

The risk model comprises 10 differentially expressed IRGs (DEIRGs) that correlate with prognosis. We demonstrate that the risk model based on the immune-related genes is an accurate predictor of UBC survival by validating it in 2 independent UBC datasets. Multivariate Cox regression analysis revealed that the risk score is an independent predictor for UBC prognosis. Our study also reveals correlation between the risk model and several clinicopathological factors including invasive status, clinical stage, and T and N stage. Nomogram analysis using risk score, clinical characteristics and patient information revealed a high predictive performance by the risk model. These results indicate that the signature can serve as an accurate prognostic tool.

Among these risk genes, 5 DEIRGs (AHNAK, CALR, PDGFRA, RAC3, PLA2G2A) correlated with high risk and 5 (OAS1, CDH1, SEMA3F, IRF5, CARD11) were protective factors. AHNAK1 is a scaffold protein, highly expressed by CD4+ T cells, and is a critical component for calcium signaling. AHNAK1 deficiency resulted in a reduced calcium influx upon TCR crosslinking [32]. Calreticulin (CALR) is a Ca2+-binding endoplasmic reticulum (ER) protein that contributes to the initiation of adaptive anticancer immunity in the context of immunogenic cell death (ICD), an immune response that eradicates cells experiencing damage beyond recovery in support of organismal fitness [33]. Phospholipase A2 group IIA (PLA2G2A) is an antimicrobial molecular which shows bactericidal activity against bacteria. Upregulation of PLA2G2A in cancer cell significantly suppressed L. monocytogenes infection. It has been suggested that the anti-bacterial activity of Pla2g2a is involved in tumor inhibition [34]. OAS1 is a member of OAS family which acts as antiviral enzymes induced by interferon. It was reported that neutrophils, which contribute to tumor metastasis through multiple pathways, were the top tumor immune infiltrating cell type associated with OAS in in breast cancer [35]. CDH1 deletions have been shown to cause cancers with immune cell infiltration, activation of targetable immune checkpoint pathways and gene expression related to T-regulatory (Treg) cell signaling [36]. Platelet-derived growth factor receptor alpha (PDGFRA) belongs to the type III receptor tyrosine kinase family. It has been reported that imatinib downregulates PD-L1 and IRF1 expression through the inhibition of KIT and PDGFRA, thus contributing to counteract the suppressed adaptive immune response against GIST [37]. SEMA3F is a member of Semaphorins and is found to regulate immune cell trafficking during the development of thymocytes, and it is also found to regulate NK-cell migration and NK–DC interactions [38]. RAC3 was reported to show antitumor activity through maintaining the maturation and effector function of NK cells via modulation of several critical T-bet-dependent genes such as IRF5 [39]. Transcription factor interferon regulatory factor 5 (IRF5) has been shown to regulate the expression of genes involved in the stimulation of the immune system. The GM-CSF-IRF5 signaling axis has been reported to promote antitumor immunity through activation of CD8 T cell infiltrates [40]. CARD11-BCL10-MALT1 (CBM) signaling mediates TCR-induced NF-κB activation in Tregs and controls the conversion of resting Tregs to effector Tregs and disruption of the CBM complex is sufficient to prime the tumor environment for successful immune checkpoint therapy [41]. AHNAK was identified as a novel candidate biomarker [42] for UBC diagnosis by liquid-based cytology and its prognostic value was demonstrated in other studies [43]. Somatic CDH1 mutations are frequent in UBC and promote UBC progression [44]. CDH1 deletions have been shown to cause cancers with immune cell infiltration, activation of targetable immune checkpoint pathways and gene expression related to T-regulatory (Treg) cell signaling, suggesting immune checkpoint potential therapeutic options [36]. Although the remaining 8 genes have not been associated with UBC, their roles and their prognostic value in UBC merit investigation.

Given that multiple genes mentioned above modulate the tumor microenvironment in cancers and that TIICs may contribute to UBC progression [45, 46], we explored the proportions of TIICs in UBC and the correlation of risk status with various immune cell subtype. The proportion of CD8+ T cells, Tregs, and dendritic cells were lower in the high-risk group according to the 10-gene risk model. High CD8+ T cells and dendritic cells proportions correlated with better survival. In contrast, the proportion of macrophages and neutrophils were high in the high-risk set, in which the high proportion of macrophages indicated poorer OS. Moreover, high CD8+ T cells and macrophage infiltration correlated with longer and shorter UBC survival, respectively, which is consistent with previous studies [47, 48]. Although neutrophils do not correlate with OS, they have been reported to correlate with negative prognostic in some cancer types [49]. Nonetheless, the function of Tregs and monocytes in UBC prognosis remains controversial. Changes in tumor immune cells may partly explain the prognostic value of this model.

To investigate the mechanisms underlying this risk model, we analyzed the UBC mutation profiles. Interestingly, we uncovered a subgroup with low-risk scores that featured higher FGFR3 mutation rates. Moreover, for the luminal papillary subtype most patients have low risk scores [5]. The luminal papillary subtype cohort was characterized by papillary morphology, lower stage, lower risk for progression and enrichment with FGFR3 mutations. These results suggest tyrosine kinase inhibitors against FGFR3 as a therapeutic option in patients with low risk scores. Patients with high risk score accounted for the bulk of basal-squamous and neuronal subtypes, which had high rates of TP53 and RB1 mutations, suggesting the worse survival, making it significant to recognize clinically. In addition, high levels of immune markers and mesenchymal expression signature in basal-squamous subtype suggested resistance to cisplatin-based chemotherapy, and sensitivity to immune checkpoint inhibition. However, despite high PD1/L1 expression in high-risk tumors, glycolytic pathways also simultaneously enriched. Cancer cells are known to reprogram their metabolism by upregulating glucose uptake. Enhanced glycolysis within the TME results in an accumulation of lactic acid [50] and glucose-deprived niche can directly inhibit glycolysis in immune cells, thus hindering antitumor immune responses. Taken together, these results suggest that the immune-related risk signature may help personalize therapies. Glycolysis inhibition coupled with immunotherapy design may be a rational approach to improve the efficacy of ICIs.

TFs bind to promoter regions of specific genes to modulate their expression. To identify the modulators of these immune-related genes, we constructed a regulatory network of DEGs between the TF and IRG genes. This analysis revealed KLF5, SOX4, GRHL2 and GATA3 as nodes with highest connectivity, suggesting them as core genes regulators. KLF5 is a known oncogene in bladder cancer that promotes angiogenesis and invasion [43]. Low GATA3 levels have been reported to regulate self-renewal and tumorigenicity in bladder cancer stem cells by modulating STAT3 expression [51].

A limitation of this study is its reliance on public datasets, whose accompanying information is limited. Thus, these findings need to be validated in future clinical studies. Moreover, the proportion of Asian patients in our study was small. Future studies should incorporate more Asian UBC samples.

In conclusion, we have developed and validated a reliable 10-IRG signature model for predicting UBC patient outcomes and potential therapeutic options. Our data suggest that changes in mutation landscape, TIME and glycolytic status might be probable causes for this signature’s prognostic capacity. These findings may unveil novel, potential immune biomarkers and therapeutic strategies. Additional studies are needed to confirm the prognostic value of this signature.

Materials and Methods

Patient data and databases

A comprehensive list of IRGs was obtained from Immunology Database and Analysis Portal (ImmPort) database (https://immport.niaid.nih.gov) [51]. UCSC database (https://xena.ucsc.edu/) was utilized to obtain gene expression data of normal bladder samples from the GTEx database. Data on normal and UBC tissues was downloaded from the Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/). To get a unified form, the RNA-seq data in the 2 databases were converted into the log2(x+1) form and then normalized. Clinical data was downloaded and organized from TCGA. The gene expression data and corresponding clinical information of microarray dataset from GSE13507, which used the Illumina human-6 v2.0 expression beadchip, were downloaded from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/).

Differential gene analysis

Analysis of differential gene expression between normal and UBC tissues was done using the Wilcoxon signed-rank test, applying FDR <0.05 and |log2fold change (FC)|>1 as cutoff threshold. Intersection between the DEG list and IRG list or TF list was established to identify differentially expressed immune-related genes (DEIRGs).

Functional enrichment analysis

In order to explore the underlying molecular mechanisms of DEGs, GO term and KEGG pathway analysis were done using clusterProfiler R package and were visualized using “ggplot2” [52] on R. Protein-protein interaction enrichment analysis was done on Metascape (http://www.home-for-researchers.com) [53] using OmniPath8, BioGrid6 and InWeb_IM7. Closely connected network components were identified using the Molecular Complex Detection (MCODE) [54].

Construction and validation of an immune-related prognostic model

UBC cases in the TCGA dataset that had overall survivals >0 days were randomly classified into a training and a testing group. The training group was used to identify the prognostic DEIRGs using univariate Cox regression analysis (p≤0.05), to establish a prognostic immune-related risk model using Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis to minimize overfitting, and to seek out the optimal gene pattern using glmnet on R [55]. Risk scores for each patient in the training set and testing set were calculated using the formula: Risk score = exprgene1* coefficient gene 1+exprgene 2 * coefficient gene 2+ ⋯ + exprgene10 * coefficient gene 10. UBC cases in the training group were grouped into high and low-risk categories based on the median risk score. The signature’s prognostic value was tested on the TCGA and GSE13507 cohorts. To validate the prognostic power of the IRG risk model, area under the curve (AUC) was calculated with the timeROC R package. The discrimination of the risk models was measured and compared by Harrell’s concordance index (c-index). Survival analysis was performed using the survminer R package and visualized using ggplot2 R package.

TF-IRG network construction

To explore potential relationships between TF and IRGs, correlation analysis of DEGs between the 2 gene lists was carried out. Significant results in correlation analysis were used to construct a regulatory network, which was visualized using Cytoscape. Correlation values between TF and IRGs >0.3 were considered significantly correlated. Based on the co-expression network, we analyzed the regulatory relationship of genes and identified core regulatory TF genes.

Estimate of tumor-infiltrating immune cells

RNA-seq gene expression data on bladder cancer dataset (TCGA-BLCA) from TCGA database were used to estimate the relative proportion of infiltrated immune cells utilizing the CIBERSORT R package and TIMER [29, 56].

Mutation analysis

Mutation annotation format (MAF) containing somatic variants data was downloaded from the TCGA database. MAF files were then visualized and summarized from this study using the maftools Bioconductor package [57].

Statistical analysis

Continuous variables were described as mean ±S.D., while categorical variables were presented by frequency (n) and proportion (%). Statistical analysis was performed using the R software or GraphPad Prism version 6.0. p<0.05 was considered statistically significant. Univariate and multivariate analyses were performed using the Cox proportional hazards regression model to evaluate the prognostic effect of our immune-related signature and other clinicopathological features. Time-dependent ROC analysis was utilized to assess the accuracy of the immune-related prognostic model. The log-rank test in Kaplan–Meier survival curves was performed to analyze differences in overall survival (OS).

Supplementary Materials

Supplementary Figures

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 Treg 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]