The tumor microenvironment and prognostic role of autophagy- and immune-related genes in bladder cancer

Abstract

BACKGROUND:

Autophagy-related genes and immune-related genes contribute significantly to the initiation and prognosis of bladder cancer (BLCA).

OBJECTIVE:

We aimed to explore differentially expressed autophagy-related genes (DEARGs) and immune-related genes (DEIRGs) in BLCA to create a prognostic risk assessment model and gain some insights into BLCA’s molecular underpinnings.

METHODS:

The prognostic DEARGs and DEIRGs were evaluated for BLCA through The Cancer Genome Atlas (TCGA) database ( $n=$ 399) and GSE13507 dataset ( $n=$ 165). The BLCA risk model was constructed and verified. The immune score, stromal score, and estimate score in different risk groups were calculated by the ESTIMATE algorithm. Immune infiltration levels were assessed by a single sample gene set enrichment analysis (GSEA) algorithm.

RESULTS:

In the risk model, AURKA, ACTC1, MYLK, PDGFD, PDGFRA and TNC were significantly associated with the overall survival. The pathways in cancer, T cell receptor signaling pathway and B cell receptor signaling pathway were significantly gathered in the high-risk group. Moreover, the risk score was significantly correlated with infiltrating immune cells, expression of critical immune checkpoints and mismatch repair genes including MSH6, MLH1, and MSH2.

CONCLUSIONS:

In this study, three DEARGs (AURKA, ACTC1, MYLK) and three DEIRGs (PDGFD, PDGFRA, TNC) were demonstrated to be potential prognostic biomarkers for BLCA patients through bioinformatics methods, which might be novel therapeutic targets and prognostic markers for BLCA, in follow up studies, we will combine experiments to verify this.

Keywords

Bladder cancer autophagy immune risk model prognosis

1. Introduction

Bladder cancer (BLCA) is the major malignant tumor of the urinary system and one of the tenth most common new malignant tumors worldwide [1]. Approximately 5% of BLCA patients have metastatic disease at the time of diagnosis [2]. Metastatic bladder cancer is mainly treated by chemotherapy based on platinum and gemcitabine, partly by radiotherapy and immunotherapy [3, 4]. Despite these aggressive therapies, the survival of patients with BLCA still faced severe challenges [5] and the prognosis of BLCA patients with recurrence and metastasis is still poor [6]. It is needed to further clarify the mechanism of disease development and develop new diagnosis and treatment measures to improve clinical diagnosis and treatment of BLCA. Therefore, there is an urgent need to explore and determine promising prognostic indicators in BLCA, to further improve the survival rate and prognosis of BLCA.

Autophagy plays an important role in inhibiting tumorigenesis and tumor development by maintaining genomic integrity, reducing proliferation and inhibiting inflammation [7]. Autophagy is closely related to the occurrence and development of cancers [8]. Some autophagy genes such as Jun, myc and itga3 can be used as potential prognostic markers of BLCA [9]. Numerous studies have reported the influence of the immune microenvironment on BLCA development and immunotherapy including programmed cell death protein 1 (PD-1)/programmed cell death-ligand 1 (PD-L1) blockade, which was extensively applied to the treatment of BLCA patients [10, 11]. The medicines of immunotherapy and regulating autophagy bring new opportunities for patients with BLCA [12, 13].

However, at present, effective prognostic markers of immunotherapy are lacking. So, it is a new idea for clinical treatment of BLCA to study the relationship between autophagy, immunity and BLCA, which effectively predicts the prognosis of BLCA. Although there have been reports on autophagy or immune genes related to BLCA prognosis. There is no study on autophagy combined with immune genes as prognostic markers of BLCA and their relationship with immunotherapy response. It is significant for the treatment and prognosis of BLCA to explore the genes and construct a gene model related to the prognosis of BLCA.

In this study, 32 differentially expressed autophagy-related genes (DEARGs) or immune-related genes (DEIRGs) were identified from the cancer genome atlas (TCGA) dataset and GSE13507 dataset of BLCA. We identified three DEARGs (AURKA, ACTC1, MYLK) and three DEIRGs (PDGFD, PDGFRA, TNC) that were associated with the overall survival (OS) of BLCA patients and suggested as potential prognostic markers of BLCA. This may provide a foundation for developing highly specific prognostic tools and more robust therapeutic strategies for BLCA in the future.

2. Materials and methods

2.1 Data source and preprocessing

We collected the transcriptome data of 424 samples (include 19 normal and 405 BLCA samples) from the TCGA database (https://portal.gdc.cancer.gov/). 399 BLCA patients with complete clinical information in the 405 BLCA samples were selected serve as a training set. GSE13507 dataset containing 256 samples (include 68 normal and 188 BLCA samples) was downloaded from the gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). 165 BLCA patients with complete clinical information in the 188 BLCA samples were selected as a validation set. The autophagy-related genes (ARGs) were downloaded from the Human Autophagy Database (HADb) and Molecular Signatures Database (MSigDb) in terms of “GO_AUTOPHAGY”, then integrated into a set containing 743 ARGs after removing duplication. 456 immune-related genes (IRGs) were downloaded from the ImmPort database.

2.2 Identification of DEARGs and DEIRGs

The differentially expressed genes (DEGs) between normal and BLCA samples in the TCGA dataset and GSE13507 dataset were determined respectively by using the ‘limma’ R package [14, 15]. The $|$ log2 Fold Change (FC) $|$ $>$ 1 and $P$ -value $<$ 0.05 were regarded as the cut-off criteria for DEGs. Heatmaps were drawn to visualize the differences in gene expression levels between the two groups. DEARGs/DEIRGs were identified by taking the intersection of DEGs from the TCGA dataset, GSE13507 dataset, and ARGs/IRGs.

Figure 1.

Differentially expressed ARGs/IRGs. (A) Heatmap demonstrating DEGs identified in the TCGA dataset between BLCA and normal samples, with red representing high expression and green representing low expression; (B) Heatmap demonstrating DEGs identified in the GSE13507 dataset between BLCA and normal samples, with red representing high expression and green representing low expression; (C) Venn graph of 32 DEARGs/DEIRGs, of which 14 ARGs and 18 IRGs, were identified by taking the intersection of DEGs from the TCGA dataset, GSE13507 dataset, and ARGs/IRGs.

2.3 Construction and verification of the prognostic signature

The univariate Cox regression model was used to identify DEARGs/DEIRGs that were significantly associated ( $P<$ 0.2) with the overall survival (OS) of the BLCA patients in the training set. Subsequently, the candidate DEARGs/DEIRGs were subjected to multivariate Cox regression and stepwise analysis to evaluate their contribution as prognostic factors in patient survival. The risk score calculating formula was: Riskscore $=$ h0 (t) * exp ( $\beta$ 1X1 $+$ $\beta$ 2X2 $+\ldots+$ $\beta$ nXn). In this formula, $\beta$ refers to the regression coefficient, and the value of HR can be obtained by taking its inverse logarithm. X1 represented the gene expression level. h0 (t) is the benchmark risk function. The training set was stratified into high- and low-risk groups by the median risk score using the ‘survminer’ package in R software. Next, the K-M survival curve and log-rank test were performed to evaluate the OS between high- and low-risk groups. The area under the curve (AUC) of receiver operating characteristic (ROC) was calculated using the ‘survival ROC’ package in R software. In addition, the risk plot was illustrated using the ‘pheatmap’ package in R software. We used the validation set to perform the same procedure to validate the risk model.

2.4 Correlation analysis of prognostic gene signature and clinicopathological characteristics

The risk scores between subgroups were compared based on the following clinicopathological features: age ( $<$ 60 or $\geqslant$ 60 years), gender (male or female), stage (I, II, III, and IV), T stage (T0, T1, T2, T3, T4), N stage (N0, N1, N2, N3), and M stage (M0, M1).

2.5 Gene set enrichment analysis

The gene set enrichment analysis (GSEA) was carried out in high- and low-risk groups with c2.cp.kegg.v7.0. symbols derived from MSigDb as the gene set via GSEA software (3.18.0) [16]. The significantly enriched KEGG pathway was screened with a NOM $p$ -value $<$ 0.05 as the threshold.

2.6 Evaluation of immune microenvironment

The immune-, stromal- and estimate-scores of each TCGA-BLCA sample in the low- and high-risk groups were calculated by applying the ‘Estimation of Stromal and Immune cells in Malignant Tumours using Expression data’ (ESTIMATE) algorithm. The enrich score of 22 immune cells was calculated with single-sample gene set enrichment analysis (ssGSEA) [17]. A box plot was drawn to visualize the differences in 22 immune cells infiltration between the low- and high-risk groups.

2.7 Statistical analysis

All analyses were conducted using the R programming language, and the data from different groups were compared by the Wilcoxon test. If not specified above, a $P$ -value less than 0.05 was considered statistically significant.

3. Results

3.1 Identification of DEARGs and DEIRGs in BLCA

The expression level of genes was compared between normal and BLCA samples in the TCGA dataset and GSE13507 dataset respectively. We found that 3134 DEGs were identified in the TCGA dataset, of which 1310 were up-regulated and 1824 were down-regulated(Supplemental Table 1). A total of 492 DEGs were identified in the GSE13507 dataset, including 74 up-regulated and 418 down-regulated genes(Supplemental Table 2). The heatmaps showed the expression levels of these DEGs (Fig. 1A and B). A total of 32 DEARGs/DEIRGs, of which 14 ARGs and 18 IRGs, were identified by taking the intersection of DEGs from the TCGA dataset, GSE13507 dataset, and ARGs/IRGs (Fig. 1C, Supplemental Table 3).

Figure 2.

Univariate (A) and multivariate (B) Cox regression analysis in terms of OS for patients with BLCA.

Figure 3.

The evaluation of the BLCA risk model. (A) The Kaplan-Meier curves of OS for patients with high-risk scores (red line) and low-risk scores (blue line) in the training set; (B) The Kaplan-Meier curves of OS for patients with high-risk scores (red line) and low-risk scores (blue line) in the validation set; (C) The AUC values of 1-, 3-, and 5-year survival in the training set; (D) The AUC values of 1-, 3-, and 5-year survival in the validation set; (E) Distribution of patients with high-risk scores (red color) and low-risk scores (green color) in the training set; (F) Survival status of patients with BLCA (red dots stand for the deceased patients and the green dots stand for the survivors) in the training set; (G) Distribution of patients with high-risk scores (red color) and low-risk scores (green color) in the validation set; (H) Survival status of patients with BLCA (red dots stand for the deceased patients and the green dots stand for the survivors) in the validation set; (I) The heatmap of the six survival-associated DEARGs and DEIRG expression profiles in the training set; (J) The heatmap of the six survival-associated DEARGs and DEIRG expression profiles in the validation set.

Figure 4.

The relationships between the riskscore and clinicopathological factors. (A) age; (B) gender; (C) pathological stage; (D) T stage; (E) M stage and (F) N stage in the BLCA patients from the TCGA database. ${}^{****}P$ -value $\leqslant$ 0.0001, ${}^{***}P$ -value $\leqslant$ 0.001, ${}^{**}P$ -value $\leqslant$ 0.01, ${}^{*}P$ -value $\leqslant$ 0.05, ns $P$ -value $>$ 0.05.

Figure 5.

The prognostic values of the risk score in forest plot of the hazard ratios with high-risk genes (HR $>$ 1, red) and low-risk genes (HR $<$ 1, green). (A) The univariate analysis indicated that the pathological stage, T stage, N stage, and risk score were associated with the prognosis of BLCA patients; (B) The multivariate Cox regression analyses demonstrated that the risk score can be used independently to predict the prognosis of BLCA patients.

Figure 6.

The KEGG pathways significantly enriched in low- and high-risk groups. (A) Pathways in cancer; (B) T cell receptor signaling pathway; (C) B cell receptor signaling pathway; (D) Oxidative phosphorylation; (E) Ribosome; (F) Glycerophospholipid metabolism.

Figure 7.

Comparison of the immune microenvironment and immunotherapy efficiency between low-risk group (blue) and high-risk group (yellow) in TCGA dataset. (A) immune score; (B) stromal score; (C) estimate score; (D) enrichment scores of 22 immune-associated gene sets using the ssGSEA algorithm; (E) expression of immune checkpoint molecules. ${}^{****}P$ -value $\leqslant$ 0.0001, ${}^{***}P$ -value $\leqslant$ 0.001, ${}^{**}P$ -value $\leqslant$ 0.01, ${}^{*}P$ -value $\leqslant$ 0.05, ns $P$ -value $>$ 0.05.

Figure 8.

Correlations of risk score with MMR genes (MSH6, MLH1, PMS2, MSH2). (A) MMR genes expression in low- (yellow) and high-risk groups (blue) in TCGA; (B) MMR genes expression in low- (yellow) and high-risk groups (blue) in GSE13507 dataset; (C) Pearson correlation of risk score with MMR genes in TCGA; (D) Pearson correlation of risk score with MMR genes in GSE13507 dataset. ${}^{****}P$ -value $\leqslant$ 0.0001, ${}^{***}P$ -value $\leqslant$ 0.001, ${}^{**}P$ -value $\leqslant$ 0.01, ${}^{*}P$ -value $\leqslant$ 0.05, ns $P$ -value $>$ 0.05.

3.2 Construction and verification of the BLCA risk model

Through univariate Cox regression analysis, 15 out of the 32 DEARGs/DEIRGs were identified as closely related to BLCA patients’ OS in the training set ( $P<$ 0.2) (Fig. 2A). Then, 15 DEARGs/DEIRGs were further submitted to stepwise multivariate Cox regression analysis, and three DEARGs (AURKA, ACTC1, MYLK) and three DEIRGs (PDGFD, PDGFRA, TNC) were identified as prognostic genes (Fig. 2B).

After erecting the risk model, the patients were divided into high- and low-risk groups based on the median value of the risk score. Then K-M survival analysis was performed in the training set and validation set. The results suggested that patients with high-risk scores had remarkably worse survival rates than those with low-risk scores (Fig. 3A and B). As shown in Fig. 3C and D, the AUC values of 1-, 3-, and 5-year survival in the training set and validation set were greater than 0.6 respectively, exhibiting a decent accuracy. Furthermore, we ranked all BLCA patients according to their risk scores to evaluate the survival rate distribution. The scatterplot showed the survival status of patients with a different risk scores (Fig. 3E–H). We described the expression level of prognostic genes with the different risk scores of samples using heatmaps (Fig. 3I and J). As shown in Fig. 3I and J, the expression of AURKA, PDGFD, PDGFRA, and TNC was higher in the high-risk group. In summary, these results indicated the risk score showed satisfactory performance in predicting the OS for BLCA.

3.3 Stratification analysis of the risk score

For determining the relationship between the risk score with the clinicopathological characteristics, we analyzed the corresponding clinical information of BLCA patients from the TCGA database. The results revealed that risk scores had significant differences in different age groups, gender groups, pathological stage groups, and T stage groups (Fig. 4A–D), while the risk score was not associated with M and N stages (Fig. 4E and F). The result of the stratification analysis showed that the prognostic signature could be a potential prognostic marker for some subgroups of BLCA patients.

3.4 Independent prognostic value of the risk model

We performed univariate and multivariate Cox regression analyses to assess whether the risk score was independent of other clinical parameters (age, gender, pathological stage, TNM stage) as a prognostic factor for BLCA patients. The univariate analysis indicated that the pathological stage, T stage, N stage, and risk score were associated with the prognosis of BLCA patients (Fig. 5A). The multivariate Cox regression analyses demonstrated that the risk score can be used independently to predict the prognosis of BLCA patients (Fig. 5B).

3.5 Gene set enrichment analyses between low- and high-risk groups

To unravel the molecular mechanism underlying the six-gene signature, GSEA was conducted in the TCGA database. The KEGG pathways enriched in low- and high-risk groups were shown in the Supplementary Table 4. We noted that some cancer and immune-related pathways, including pathways in cancer, T cell receptor signaling pathway and B cell receptor signaling pathway were significantly enriched in the high-risk group as defined (Fig. 6A–C). In the low-risk group, the significantly enriched pathways were mainly related to oxidative phosphorylation, ribosome, and glycerophospholipid metabolism (Fig. 6D–F).

3.6 Comparison of the immune microenvironment and immunotherapy efficiency between low- and high-risk groups

To further explore the correlation between the risk score and immune status, the ESTIMATE algorithm was used to calculate the immune score, stromal score, and estimate score for the tumor tissues. The results show that the immune score, stromal score, and estimate score were higher in the high-risk group compared with the low-risk group (Fig. 7A–C). Moreover, we estimated the enrichment scores of 22 immune-associated gene sets between low- and high-risk groups using the ssGSEA algorithm. As shown in Fig. 7D, the low-risk BLCA patients had significantly higher enrich scores of Th17 cells than high-risk BLCA patients and the high-risk BLCA patients had significantly higher enrich scores of aDC, B cells, CD8 T cells, cytotoxic cells, DC, Eosinophils, iDC, Macrophages, Mast cells, Neutrophils, NK CD56dim cells Tcm, NK cells, T cells, T helper cells, Tem, TFH, Tgd, Th1 cells, and Th2 cells than low-risk BLCA patients, These findings suggested that there were significant differences in immune microenvironment between the two risk groups.

Immune checkpoint molecules are essential for immune function and have diverse clinical implications in immunotherapy. Therefore, we contrasted the expression of immune checkpoint molecules between the two risk groups. As shown in Fig. 7E, the expression of immune checkpoint molecules PDCD1, CD274, CTLA4, IDO2, IDO1, LAG3, TIGIT and BLTA in the high-risk group was greater than that in the low-risk group in the TCGA dataset.

3.7 Correlations of risk score with MMR genes

Mismatch repair (MMR) proteins are crucial for DNA mismatch repair and play important role in maintaining the integrity and stability of genetic material. As shown in Fig. 8A and B, three MMR genes, MSH6, MLH1, and MSH2 present higher expression in high-risk groups in TCGA and GSE13507 datasets. Subsequently, we analyzed the Pearson correlation of risk score with MMR genes (MSH6, MLH1, PMS2, MSH2) and showed that risk score was significantly correlated with MSH6, MLH1, and MSH2 (Fig. 8C and D).

4. Discussion

At present, the molecular mechanism of BLCA has not yet been described, while growing evidence has revealed that ARGs, IRGs and immune cell infiltration, play critical roles in the pathogenesis and progression of BLCA [18, 19]. it is a new idea for clinical treatment of BLCA to study the relationship between autophagy, immunity and BLCA, which effectively predicts the prognosis of BLCA. In this study, our goal was to identify common DEARGs and DEIRGs in BLCA. We had identified three DEARGs (AURKA, ACTC1, MYLK) and three DEIRGs (PDGFD, PDGFRA, TNC) were associated with the development and OS of BLCA patients.

In our study, the BLCA risk score model was constructed based on the above three DEARGs and three DEIRGs, and verified to be effective to predict the prognosis of BLCA patients. Compared with other risk score models of BLAC [20, 21], our model included not only immune genes but also autophagy genes. Independent prognostic analysis shown that our risk score could be used as a clinically independent prognostic factor, which provided a reference for doctors to decide clinical treatment strategies. In GO and KEGG enrichment analysis, cancer related pathways were enriched in the high-risk group, which further illustrated the effectiveness of our risk scoring model. Moreover, the enrichment of immune related pathways in high-risk group indicated that immunity played an important role in affecting the prognosis of BLCA patients.

The prognosis of BLCA patients with different infiltration of immune cells was variant, and these differences affected development of BLCA and sensitivity of immune checkpoint inhibitors (ICIs) treatments [10]. The number of tumor infiltrating lymphocytes (TILs) was closely related to the prognosis of BLCA patients [22]. In this study, we found that the TILs and expression of critical immune checkpoints were significantly different between high-risk group and low-risk group. The immune score, stromal score, and estimate score were higher and the expression of many immune checkpoint molecules was greater in the high-risk group than those in the low-risk group. These results suggested that immunotherapy would be more suitable for BLCA patients in the high-risk group. Tumor-associated macrophages affected progress of solid tumors as part of tumor microenvironmental cells [23]. The macrophages directly affected immune response to BLCA and were the main inhibitors of the tumor microenvironment that promoted growth and metastasis of BLCA [24]. The paracrine interaction between tumor-associated macrophages and BLCA cells led to enhanced BLCA cell growth, migration, and invasiveness, and moreover, increased CCL5 cytokine production and IL-8 in tumor microenvironment [25]. In our study, the enrich score of macrophages was significantly higher in high-risk score group than that in low-risk group. A higher abundance of macrophages in the tumor immune microenvironment suggested higher amounts of abnormal expression of immune genes; this, in turn, participated in the occurrence and progression of BLCA and the processes of invasion and metastasis [24, 26].

MMR played important role in maintaining the integrity and stability of genetic material. In this study, there were significant differences in MMR genes between the high- and low-risk groups. The risk score was significantly correlated with multiple MMR genes. At present, MMR defect or high instability of microsatellite was one of the markers predicting the efficacy of PD1 antibody [27]. BLCA lacking MMR genes was usually immunogenic and showed extensive infiltration of T cells, making it highly responsive to ICIs [28, 29].

Viewed in total, this study found biomarkers of autophagy- and immune-related genes that could be used as prognostic factors for BLCA. The risk model established was highly correlated with autophagy and immunity and provided valuable targets and prognostic indicators for BLCA patients. These DEARGs and DEIRGs in the model played an important role in the occurrence and development of BLCA, and may be new targets for precise treatment of BLCA. This study provides a theoretical basis for judging the prognosis and treatment response of BLCA patients. However, the study was conducted based on public databases and did not include experimental verification. Moreover, this study had another limitation of small sample size. Notwithstanding these limitations, this study can still provide a theoretical basis for further study of prognosis-related biomarkers of BLCA, and we will continue to pay attention to the role of these genes.

5. Conclusion

AURKA, ACTC1, MYLK, PDGFD, PDGFRA and TNC can be used as novel prognostic markers for BLCA, and they are related to immune cell infiltration and immune checkpoint. The BLCA prognosis model based on DEARGs and DEIRGs is reliable, and is expected to help clinical physicians to predict the prognosis of BLCA. At present, this study is limited to small sample size and bioinformatics clinical data from public databases, in subsequent studies, we will verify it in more samples and combination with clinical trials. Experiments of in vitro and in vivo were lacked to validate our data mining results, which are also needed in the future.

Author contributions

Conception: Cheng Chen and Zhenhua Gao.

Interpretation or analysis of data: Cheng Chen, Zhenhua Gao and Peng Gu.

Preparation of the manuscript: Cheng Chen, Zhenhua Gao and Jianheng Chen.

Revision for important intellectual content: Jianheng Chen, Xiaodong Liu and Jihong Shen.

Supervision: Xiaodong Liu and Jihong Shen.

Supplementary data

The supplementary files are available to download from https://dx-doi-org.web.bisu.edu.cn/10.3233/CBM-220058.

Footnotes

Acknowledgments

This study was partly supported by grants from Joint Special Funds for Applied Basic Research of Yunnan Provincial Science and Technology Department and Kunming Medical University (Grant No. 202001AY070001-206 and 202101AY070001-014), and the National Natural Science Foundation of China (Grant Nos. 81960133, 81860451, and 81802548).

References

Siegel

R.L.

Miller

K.D.

and Jemal

, Cancer statistics, 2020, CA-Cancer Journal for Clinicians 70 (2020), 7–30. doi: 10.3322/caac.21590.

Bethesda, Cancer Stat Facts: Bladder Cancer, NIH NCI: Surveillance, epidemiology, and end results program; 2021. URL: https://seer.cancer.gov/statfacts/html/urinb.html.

Alfred Witjes

Lebret

Compérat

E.M.

Cowan

N.C.

Santis

M.D.

Bruins

H.M.

Hernández

Espinós

E.L.

Dunn

Rouanne

Neuzillet

Veskimäe

van der Heijden

A.G.

Gakis

and Ribal

M.J.

, Updated 2016 EAU guidelines on muscle-invasive and metastatic bladder cancer, European Urology 71 (2017), 462–475. doi: 10.1016/j.eururo.2016.06.020.

Martinez Rodriguez

R.H.

Buisan Rueda

and Ibarz

, Bladder cancer: Present and future, Medicina Clinica (Barc) 149 (2017), 449–455. doi: 10.1016/j.medcli.2017.06.009.

Guo

Yao

Wang

and Li

, Non-coding RNA in bladder cancer, Cancer Letters 485 (2020), 38–44. doi: 10.1016/j.canlet.2020.04.023.

Aggen

D.H.

and Drake

C.G.

, Biomarkers for immunotherapy in bladder cancer: A moving target. Journal for ImmunoTherapy of Cancer 5 (2017), 94. doi: 10.1186/s40425-017-0299-1.

Cheng

Ren

Hait

W.N.

and Yang

J.M.

, Therapeutic targeting of autophagy in disease: Biology and pharmacology, Pharmacological Reviews 65 (2013), 1162–1197. doi: 10.1124/pr.112.007120.

Yun

C.W.

and Lee

S.H.

, The roles of autophagy in cancer, International Journal of Molecular Sciences 19 (2018), 3466. doi: 10.3390/ijms19113466.

Wang

S.S.

Chen

S.H.

Pang

J.S.

Cai

K.T.

Yan

H.B.

Huang

Z.G.

and He

R.Q.

, Identification and validation of an individualized autophagy-clinical prognostic index in bladder cancer patients, Onco Targets and Therapy 12 (2019), 3695–3712. doi: 10.2147/OTT.S197676.

10.

Zhu

Liu

Chen

Cui

Guo

Zhou

Zhu

and Shi

, Profiles of immune infiltration in bladder cancer and its clinical significance: an integrative genomic analysis, International Journal of Medical Sciences 17 (2020), 762–772. doi: 10.7150/ijms.42151.

11.

Patel

V.G.

W.K.

and Galsky

M.D.

, Treatment of muscle-invasive and advanced bladder cancer in 2020, CA-Cancer Journal for Clinicians 70 (2020), 404–423. doi: 10.3322/caac.21631.

12.

Song

Powles

Shi

Zhang

Ingersoll

M.A.

and Lu

Y.J.

, Bladder cancer, a unique model to understand cancer immunity and develop immunotherapy approaches, Journal of Pathology 249 (2019), 151–165. doi: 10.1002/path.5306.

13.

Guo

Yang

Feng

Liu

Ren

and Zhou

, Autophagy modulation in bladder cancer development and treatment (Review), Oncology Reports 42 (2019), 1647–1655. doi: 10.3892/or.2019.7286.

14.

Ritchie

M.E.

Phipson

Law

C.W.

Shi

and Smyth

G.K.

, limma powers differential expression analyses for RNA-sequencing and microarray studies, Nucleic Acids Research 43 (2015), e47. doi: 10.1093/nar/gkv007.

15.

Phipson

Lee

Majewski

I.J.

Alexander

W.S.

and Smyth

G.K

, Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression, Annals of Applied Statistics 10 (2016), 946–963. doi: 10.1214/16-AOAS920.

16.

Subramanian

Tamayo

Mootha

V.K.

Mukherjee

Ebert

B.L.

Gillette

M.A.

Paulovich

Pomeroy

S.L.

Golub

T.R.

Lander

E.S.

and Mesirov

J.P.

, Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles, Proceedings of the National Academy of Sciences of the United States of America 102 (2005), 15545–15550. doi: 10.1073/pnas.0506580102.

17.

Bindea

Mlecnik

Tosolini

Kirilovsky

Waldner

Obenauf

A.C.

Angell

Fredriksen

Lafontaine

Berger

Bruneval

Fridman

W.H.

Becker

Pagès

Speicher

M.R.

Trajanoski

and Galon

, Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer, Immunity 39 (2013), 782–795. doi: 10.1016/j.immuni.2013.10.003.

18.

Poli

Egidi

M.G.

Cochetti

Brancorsini

and Mearini

, Relationship between cellular and exosomal miRNAs targeting NOD-like receptors in bladder cancer: Preliminary results, Minerva Urologica E Nefrologica 72 (2020), 207–213. doi: 10.23736/S0393-2249.19.03297-1.

19.

Schneider

A.K.

Chevalier

M.F.

and Derré

, The multifaceted immune regulation of bladder cancer, Nature Reviews Urology 16 (2019), 613–630. doi: 10.1038/s41585-019-0226-y.

20.

Liu

Yang

Xiang

and Tang

, Development of a prognostic index and screening of prognosis related genes based on an immunogenomic landscape analysis of bladder cancer, Aging (Albany NY) 13 (2021), 12099–12112. doi: 10.18632/aging.202917.

21.

Guo

Yuan

and Zhao

, Identification of immune-related genes that predict prognosis and risk of bladder cancer: Bioinformatics analysis of TCGA database, Aging (Albany NY) 13 (2021), 19352–19374. doi: 10.18632/aging.203333.

22.

Tumeh

P.C.

Harview

C.L.

Yearley

J.H.

Shintaku

I.P.

Taylor

E.J.

Robert

Chmielowski

Spasic

Henry

Ciobanu

West

A.N.

Carmona

Kivork

Seja

Cherry

Gutierrez

A.J.

Grogan

T.R.

Mateus

Tomasic

Glaspy

J.A.

Emerson

R.O.

Robins

Pierce

R.H.

Elashoff

D.A.

Robert

and Ribas

, PD-1 blockade induces responses by inhibiting adaptive immune resistance, Nature 515 (2014), 568–571. doi: 10.1038/nature13954.

23.

Wang

Yung

M.M.H.

Ngan

H.Y.S.

Chan

K.K.L.

and Chan

D.W.

, The impact of the tumor microenvironment on macrophage polarization in cancer metastatic progression, International Journal of Molecular Sciences 22 (2021), 6560. doi: 10.3390/ijms22126560.

24.

Sharifi

Nowroozi

M.R.

Amini

Arami

M.K.

Ayati

and Mohsenzadegan

, A review on the role of M2 macrophages in bladder cancer; pathophysiology and targeting, International Immunopharmacology 76 (2019), 105880. doi: 10.1016/j.intimp.2019.105880.

25.

Huang

C.P.

Liu

L.X.

and Shyr

C.R.

, Tumor-associated macrophages facilitate bladder cancer progression by increasing cell growth, migration, invasion and cytokine expression, Anticancer Research 40 (2020), 2715–2724. doi: 10.21873/anticanres.14243.

26.

Xue

Tong

LiuAnwei

L.F.

Liu

Zeng

Xiong

Yang

Sun

and Xu

, Tumor-infiltrating M2 macrophages driven by specific genomic alterations are associated with prognosis in bladder cancer, Oncology Reports 42 (2019), 581–594. doi: 10.3892/or.2019.7196.

27.

Arora

Velichinskii

Lesh

R.W.

Ali

Kubiak

Bansal

Borghaei

Edelman

M.J.

and Boumber

, Existing and emerging biomarkers for immune checkpoint immunotherapy in solid tumors, Advances in Therapy 36 (2019), 2638–2678. doi: 10.1007/s12325-019-01051-z.

28.

D.T.

Durham

J.N.

Smith

K.N.

Wang

Bartlett

B.R.

Aulakh

L.K.

Kemberlin

Wilt

Luber

B.S.

Wong

Azad

N.S.

Rucki

A.A.

Laheru

Donehower

Zaheer

Fisher

G.A.

Crocenzi

T.S.

Lee

J.J.

Greten

T.F.

Duffy

A.G.

Ciombor

K.K.

Eyring

A.D.

Lam

B.H.

Joe

Kang

S.P.

Holdhoff

Danilova

Cope

Meyer

Zhou

Goldberg

R.M.

Armstrong

D.K.

Bever

K.M.

Fader

A.N.

Taube

Housseau

Spetzler

Xiao

Pardoll

D.M.

Papadopoulos

Kinzler

K.W.

Eshleman

J.R.

Vogelstein

Anders

R.A.

and Diaz

L.J.

, Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade, Science 357 (2017), 409–413. doi: 10.1126/science.aan6733.

29.

Crispen PL

P.L.

and Kusmartsev

, Mechanisms of immune evasion in bladder cancer, Cancer Immunology, Immunotherapy 69 (2020), 3–14. doi: 10.1007/s00262-019-02443-4.