Online ISSN: 1945-4589
Learn about our FREE Post-Publication Promotion Services
Research Paper Volume 13, Issue 15 pp 19352—19374

Identification of immune-related genes that predict prognosis and risk of bladder cancer: bioinformatics analysis of TCGA database

Liqiang Guo1, , Qiong Wu2, , Zhen Ma3, , Mingzhen Yuan4, *,&, , Shengtian Zhao5, *, ,

  • 1 Shandong Provincial Hospital Affiliated to Shandong First Medical University, Jinan, China
  • 2 Department of Gastroenterology, The Second Hospital of Shandong University, Jinan, China
  • 3 Cheeloo College of Medicine, Shandong University, Jinan, China
  • 4 Department of Urology, Shandong Provincial Hospital Affiliated to Shandong First Medical University, Jinan, China
  • 5 Binzhou Medical University, Binzhou, China
* Equal contribution

Received: May 19, 2020       Accepted: July 6, 2021       Published: July 30, 2021      

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

Copyright: © 2021 Guo 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

Background: Bladder cancer (BLCA) is the major tumor of the urinary system, and immune-related genes (IRGs) contribute significantly to its initiation and prognosis.

Results: A total of 51 prognostic IRGs significantly associated with overall survival were identified. Functional enrichment analysis revealed that these genes were actively involved in tumor-related functions and pathways. Using multivariate Cox regression analysis, we detected 11 optimal IRGs (ADIPOQ, PPY, NAMPT, TAP1, AHNAK, OLR1, PDGFRA, IL34, MMP9, RAC3, and SH3BP2). We validated the prognostic value of this signature in two validation cohorts: GSE13507 (n = 165) and GSE32894 (n = 224). Furthermore, we performed a western blot and found that the expression of these IRGs matched their mRNA expression in TCGA. Moreover, correlations between risk score and immune-cell infiltration indicated that the prognostic signature reflected infiltration by several types of immune cells.

Conclusion: We identified and validated an 11-IRG-based risk signature that may be a reliable tool to evaluate the prognosis of BLCA patients and help to devise individualized immunotherapies.

Methods: Bioinformatics analysis was performed using TCGA and ImmPort databases. Cox regression was used to identify prognostic signatures. Two external GEO cohorts and western blotting of samples were performed to validate the IRG signature.

Introduction

Bladder cancer (BLCA), with about 80,470 new patients in the United States in 2019, is the major malignant tumor of the urinary system [1]. Approximately 25% of patients are muscle-invasive or metastatic bladder cancer when they are initially diagnosed [2, 3]. Surgical resection is the main treatment for localized BLCA [4], whereas systematic chemotherapy is the preferred treatment for advanced and metastatic BLCA [5]. Despite these aggressive therapies, the five-year overall survival of bladder cancer remains less than 20% [6]. Thus, it is critical to explore alternative treatments and to determine promising prognostic indicators in BLCA.

At present, the molecular mechanism of BLCA has not yet been described, while growing evidence has revealed that immune-related genes (IRGs) and immune cell infiltration, play critical roles in the pathogenesis and progression of BLCA [7, 8]. In recent years, immune checkpoint therapies have provided a promising opportunity for the treatment of advanced BLCA patients [9]. Unfortunately, only approximately 20% of platinum-refractory and previously untreated patients may benefit from immunotherapy [10]. Therefore, it is imperative to screen and detect immunotherapy response and prognostic predictors to predict prognosis and risk of bladder cancer.

In this study, we used transcriptome data from TCGA to construct a prognostic signature of 11 differentially expressed IRGs. The prognostic IRG-based signature was further validated in two independent GEO datasets and in proteomics data from our samples. The underlying regulatory mechanisms of IRGs were explored using bioinformatics methods. This immunogenomic signature may be a reliable tool for individualized prediction of prognosis in BLCA patients.

Results

Identification of differentially expressed IRGs

Compared with normal tissues, we identified 4876 differentially expressed genes (DEGs) in BLCA tissues including 3453 upregulated and 1423 downregulated genes (Figure 1A and Supplementary Figure 1A). We extracted 120 upregulated and 140 downregulated IRGs corresponding to those identified in the ImmPort database (Figure 1B and Supplementary Figure 1B). In addition, we performed GO functional enrichment and KEGG pathway analyses in 260 differentially expressed IRGs (DEIRGs). The top ten functional annotations were shown in Supplementary Table 1. The DEIRGs were mostly enriched in cell migration, leukocyte migration, extracellular matrix, receptor complex, receptor ligand activity and cytokine activity (Supplementary Figure 2A–2C). Furthermore, cytokine–cytokine receptor interactions were enriched in the KEGG pathways (Supplementary Figure 2D).

Volcano plot of differentially expressed genes and immune-related genes. (A) Volcano plot of differentially expressed genes between bladder cancer (BLCA) and non-tumor tissues. (B) Volcano plot of differentially expressed immune-related genes between bladder cancer (BLCA) and non-tumor tissues. The green dots represent downregulated genes, the red dots represent upregulated genes, and the black dots represent genes that were not significantly differentially expressed.

Figure 1. Volcano plot of differentially expressed genes and immune-related genes. (A) Volcano plot of differentially expressed genes between bladder cancer (BLCA) and non-tumor tissues. (B) Volcano plot of differentially expressed immune-related genes between bladder cancer (BLCA) and non-tumor tissues. The green dots represent downregulated genes, the red dots represent upregulated genes, and the black dots represent genes that were not significantly differentially expressed.

Identification of prognosis-associated DEIRGs

To determine possible prognosis-associated DEIRGs, univariate Cox regression analysis of DEIRGs were performed in the present study. After screening, we identified 51 DEIRGs that significantly correlated with the overall survival of BLCA patients (Figure 2). Similar to the results for the DEIRGs, we found that these prognosis-associated DEIRGs were most enriched in cell migration, cell proliferation, extracellular matrix, receptor ligand activity and growth factor activity (Supplementary Table 2 and Figure 3A–3C). Human cytomegalovirus infection and the MAPK signaling pathway were most enriched in the KEGG analysis (Figure 3D).

Prognosis-associated DEIRGs. Forest plot of hazard ratios showing the prognostic values of IRGs.

Figure 2. Prognosis-associated DEIRGs. Forest plot of hazard ratios showing the prognostic values of IRGs.

Gene functional enrichment of prognosis-associated DEIRGs. (A) The top ten most significant biological processes in the gene ontology. (B) The top ten most significant cellular components in the gene ontology. (C) The top ten most significant molecular functions in the gene ontology. (D) The top ten most significant KEGG pathways.

Figure 3. Gene functional enrichment of prognosis-associated DEIRGs. (A) The top ten most significant biological processes in the gene ontology. (B) The top ten most significant cellular components in the gene ontology. (C) The top ten most significant molecular functions in the gene ontology. (D) The top ten most significant KEGG pathways.

Significant modular analysis based on a PPI network

A PPI network was established utilizing the 51 prognosis-associated DEIRGs, shown in Figure 4A. In the PPI network, 27 hub genes were identified by the MCODE plugin of Cytoscape. When the k-core = 2, four significant module subgroups were obtained and were showed in different colors, and the most important modules, including THBS1, PGF, SPP1, TGFB3, ELN, OXTR, PROK1, AGTR1, TACR1, and EDNRA, were marked in green (Figure 4B–4E). As shown in Figure 4F, functional annotation indicated the 27 hub genes were mostly related in sprouting angiogenesis, positive regulation of leukocyte chemotaxis, regulation of smooth muscle cell proliferation, MHC class I peptide loading complex, response to testosterone, and embryo implantation.

Significant modular analysis and functional enrichment analysis based on the PPI network. (A) Construction of a PPI network using a total of 51 prognosis-associated DEIRGs. (B) The most significant module subgroup of the hub genes, identified by MCODE plug-in, contains ten genes. (C) Module 2 contains ten hub genes. (D) Module 3 contains four hub genes. (E) Module 4 contains three hub genes. (F) Functional enrichment analysis of the 27 hub genes in the PPI network.

Figure 4. Significant modular analysis and functional enrichment analysis based on the PPI network. (A) Construction of a PPI network using a total of 51 prognosis-associated DEIRGs. (B) The most significant module subgroup of the hub genes, identified by MCODE plug-in, contains ten genes. (C) Module 2 contains ten hub genes. (D) Module 3 contains four hub genes. (E) Module 4 contains three hub genes. (F) Functional enrichment analysis of the 27 hub genes in the PPI network.

Construction of a transcription factor (TF) regulatory network

Next, we investigated the relationships between the DEIRGs in BLCA and the cancer-associated transcription factors (TFs). In total, we identified 77 differentially expressed TFs between BLCA (n = 414) and para-cancer tissues (n = 19) (Figure 5A, Supplementary Figure 3). Subsequently, using correlation scores > 0.4 and p values <0.05 as reference values, combined with 16 TFs and 51 DEIRGs, we constructed a TF regulatory network to illustrate the correlation between TFs and IRGs (Figure 5B).

TF regulatory network. (A) Volcano plot of differentially expressed TFs. The green dots represent downregulated TFs, the red dots represent upregulated TFs, and the black dots represent TFs that were not significantly differentially expressed. (B) Regulatory network of TFs and IRGs; the yellow nodes represent TFs that correlated with the IRGs, the red nodes represent IRGs with hazard ratios p 1 (p 0.4 and p

Figure 5. TF regulatory network. (A) Volcano plot of differentially expressed TFs. The green dots represent downregulated TFs, the red dots represent upregulated TFs, and the black dots represent TFs that were not significantly differentially expressed. (B) Regulatory network of TFs and IRGs; the yellow nodes represent TFs that correlated with the IRGs, the red nodes represent IRGs with hazard ratios < 1 (p < 0.05), the purple nodes represent IRGs with hazard ratios > 1 (p < 0.05) (correlation coefficient > 0.4 and p < 0.05), the green lines indicate negative regulatory relationships, and the red lines indicate positive regulatory relationships.

Construction of a prognostic risk model

Multivariate Cox regression analysis were performed to calculate a risk score for each patient as follows: Risk score = (−0.0067 ×TAP1(expression)) + (0.0004 ×MMP9(expression)) + (0.0616 ×ADIPOQ(expression)) + (0.0359 ×PDGFRA(expression)) + (0.01261 ×AHNAK(expression)) + (0.0295 ×RAC3(expression)) + (0.0066 ×OLR1(expression)) + (0.0263 ×IL34(expression)) + (0.0149 ×NAMPT(expression)) + (0.0172 ×PPY(expression)) + (−0.0720 × SH3BP2(expression)). Based on the risk scores, patients with bladder cancer were divided into high-risk (n = 186) and low-risk (n = 185) groups. Kaplan–Meier analysis revealed that the survival rate significantly favored the low-risk group (p < 0.001) (Figure 6A). The AUC of the ROC curve was 0.745, suggesting moderate effectiveness for the prognostic risk model for monitoring survival (Figure 6B). Figure 6C–6E represented the risk scores, survival status, and heatmap of the ten IRGs between the two groups.

Prognostic risk model of the cohort. (A) Kaplan–Meier curve analysis showed that patients with a high-risk score were correlated with a worse survival outcome (p B) ROC curve analysis of the prognostic risk model. (C) Risk score distribution of patients in the prognostic risk model in the cohort. (D) Survival status scatter plots for patients in the prognostic risk model. (E) Heatmap showing the distribution of the expression of the 11 immune-related genes in the cohort.

Figure 6. Prognostic risk model of the cohort. (A) Kaplan–Meier curve analysis showed that patients with a high-risk score were correlated with a worse survival outcome (p < 0.05). (B) ROC curve analysis of the prognostic risk model. (C) Risk score distribution of patients in the prognostic risk model in the cohort. (D) Survival status scatter plots for patients in the prognostic risk model. (E) Heatmap showing the distribution of the expression of the 11 immune-related genes in the cohort.

Independent prognostic value of the risk model of the cohort

We then performed univariate and multivariate Cox regression analysis to determine the efficacy of risk score derived from our prognostic risk model as an independent predictor after adjusting for other clinical parameters (Table 1). Univariate analysis revealed that age, clinical stage, tumor stage (T), lymph node (N), and risk score were significantly related to the prognosis of BLCA patients (p < 0.05). Multivariate Cox regression analysis demonstrated that the risk score was independently corelated with the overall survival in the cohort (p < 0.001) (Figure 7).

Table 1. Univariate and multivariate Cox regression analyses of the cohort.

Univariate analysisMultivariate analysis
HR95%CIP-valueHR95%CIp-value
Age1.0291.001-1.0570.045451.0220.993-1.0520.13814
Gender0.6160.352-1.0780.089530.8730.466-1.6350.67126
Stage1.8631.293-2.6830.000841.0640.519-2.1810.86637
T1.7691.193-2.6220.004531.4400.857-2.4180.16832
M2.1670.779-6.0250.138220.9070.274-3.0000.87234
N1.5731.196-2.0700.001211.2270.728-2.0680.44174
Risk score1.3981.265-1.5444.33e-111.3271.174-1.5006.28e-06
HR, hazard ratio; CI, confidence interval; T, tumor stage; M, metastasis; N, Lymph nodes.
Univariate and multivariate independent prognostic analysis of the cohort. (A) Univariate Cox regression analysis showed clinical stage, tumor, lymph nodes, and risk score were associated with the prognosis of BLCA patients. (B) Multivariate Cox regression analysis revealed that the risk score was independently associated with OS in the cohort.

Figure 7. Univariate and multivariate independent prognostic analysis of the cohort. (A) Univariate Cox regression analysis showed clinical stage, tumor, lymph nodes, and risk score were associated with the prognosis of BLCA patients. (B) Multivariate Cox regression analysis revealed that the risk score was independently associated with OS in the cohort.

External validation of the prognostic risk model in the GSE13507 and GSE32894 cohorts

To validate the reliability of our IRG-based prognostic risk model, we utilized two external validation cohorts, GSE13507 and GSE32894. Similarly, patients in the high-risk group showed poorer survival than in the low-risk group (Figure 8A, 8B). In addition, the results were also shown with a KM curve and a ROC curve (Figure 8C–8F). The cumulative results imply that the IRG-based prognostic risk characteristics based on IRGs can be used as a reliable prognostic model.

Validation of the risk signature in GSE13507 and GSE32894 datasets. (A, B) Heatmap of the 11 IRGs expression distribution, risk score distribution, and survival status between the low-risk group and high-risk group in the validation cohort. (C, D) Kaplan–Meier curve showed shorter survival time in the high-risk group patients. (E, F) ROC curve illustrated the prognostic value of the risk signature.

Figure 8. Validation of the risk signature in GSE13507 and GSE32894 datasets. (A, B) Heatmap of the 11 IRGs expression distribution, risk score distribution, and survival status between the low-risk group and high-risk group in the validation cohort. (C, D) Kaplan–Meier curve showed shorter survival time in the high-risk group patients. (E, F) ROC curve illustrated the prognostic value of the risk signature.

Validation of the IRGs in the proposed signature by western blot

To further validate the proposed signature at the protein level, we performed western blotting on bladder cancer tissues, adjacent tissues, SV-HUC-1, and T24 cells lines. Comparing the gene levels of RNA-seq from TCGA dataset (data not shown) with the expression levels of WB protein, our results indicated that the protein levels of ADIPOQ, NAMPT, TAP1, OLR1, AHNAK, PDGFRA, IL34, MMP9 and RAC3 were consistent with their RNA expression trends (Figure 9).

Validation of the IRGs by western blot. Validation of the IRGs by western blot on bladder cancer tissues and adjacent tissues, human normal bladder epithelial cells (SV-HUC-1), and a bladder cancer cell line (T24).

Figure 9. Validation of the IRGs by western blot. Validation of the IRGs by western blot on bladder cancer tissues and adjacent tissues, human normal bladder epithelial cells (SV-HUC-1), and a bladder cancer cell line (T24).

Clinical application of the prognostic model

To assess the efficacy of our model in predicting the progression of BLCA, we evaluated the correlations between the risk signatures and clinical parameters including age, sex, pathological grade and clinical stage in the TCGA data cohort (Table 2). As the levels of PPY, NAMPT, ADIPOQ, IL34, RAC3, TAP1, AHNAK, PDGFRA and the risk score increased, the pathological grade of BLCA patients increased (p < 0.01) (Figure 10A–10I). With the increase of the other factors (risk score and PDGFRA and AHNAK levels), the clinical stage and the tumor stage of patients with BLCA also increased (all p < 0.01) (Figure 10J–10O). These results demonstrate that the dysregulation of the expression of the DEIRGs is significantly associated with the development of BLCA.

Table 2. Relationships between prognostic model associated IRGs and clinical variables of patients with BLCA.

VariablesAge (≤65/≥65) t(p)Gender (Female/Male) t(p)Pathological grade (High grade/Low grade) t(p)Clinical stage (3&4/1&2) t(p)Tumor stage (T3&4/T1&2) t(p)Lymph nodes (N1&2/N0) t(p)
TAP11.132(0.259)-0.104(0.917)9.129(3.073e-13)1.298(0.197)1.006(0.316)1.492(0.137)
MMP9-0.974(0.331)-0.888(0.376)1.663(0.097)-1.16(0.247)-1.206(0.229)-1.003(0.318)
ADIPOQ-2.377(0.018)1.167(0.246)3.645(3.115e-04)-0.339(0.735)-0.662(0.509)0.233(0.816)
PDGFRA0.366(0.714)1.446(0.150)4.534(9.997e-05)-3.645(3.199e-04)-4.011(7.503e-05)-1.492(0.137)
AHNAK0.087(0.931)0.901(0.369)6.793(2.305e-07)-3.435(7.172e-04)-2.989(0.003)-2.321(0.021)
OLR1-0.016(0.987)-0.369(0.712)2.12(0.047)-1.851(0.066)-2.358(0.019)-0.388(0.698)
RAC3-1.8(0.073)0.886(0.377)4.227(1.542e-04)-1.774(0.078)-1.059(0.291)-0.581(0.562)
IL34-0.849(0.397)0.352(0.725)5.066(7.503e-07)-1.974(0.049)-2.444(0.015)-2.053(0.042)
NAMPT-0.698(0.486)2.224(0.028)4.5(1.488e-04)-1.933(0.054)-1.62(0.106)0.682(0.495)
PPY-1.385(0.167)-0.174(0.862)4.995(9.603e-07)-1.524(0.128)-1.795(0.074)0.139(0.890)
SH3BP20.983(0.327)-2.401(0.017)-0.586(0.566)2.384(0.019)1.327(0.186)1.974(0.049)
Risk Score-1.426(0.155)1.464(0.146)5.272(2.433e-07)-4.304(2.323e-05)-3.926(1.101e-04)-1.371(0.172)
t, t value determined using the Student t test; p, p value determined using the Student t test.
Relationships between prognostic-model-associated IRGs and clinical variables in the TCGA cohort. (A) PPY expression and pathological grade. (B) NAMPT expression and pathological grade. (C) ADIPOQ expression and pathological grade. (D) IL34 expression and pathological grade. (E) RAC3 expression and pathological grade. (F) TAP1 expression and pathological grade. (G) AHNAK expression and pathological grade. (H) PDGFRA expression and pathological grade. (I) Risk score and pathological grade. (J) AHNAK expression and clinical stage. (K) PDGFRA expression and clinical stage; (L) Risk score and clinical stage. (M) AHNAK expression and tumor stage. (N) PDGFRA expression and tumor stage. (O) Risk score and tumor stage.

Figure 10. Relationships between prognostic-model-associated IRGs and clinical variables in the TCGA cohort. (A) PPY expression and pathological grade. (B) NAMPT expression and pathological grade. (C) ADIPOQ expression and pathological grade. (D) IL34 expression and pathological grade. (E) RAC3 expression and pathological grade. (F) TAP1 expression and pathological grade. (G) AHNAK expression and pathological grade. (H) PDGFRA expression and pathological grade. (I) Risk score and pathological grade. (J) AHNAK expression and clinical stage. (K) PDGFRA expression and clinical stage; (L) Risk score and clinical stage. (M) AHNAK expression and tumor stage. (N) PDGFRA expression and tumor stage. (O) Risk score and tumor stage.

Inferred immune cell fractions of the high- and low-risk groups

To clarify whether our prognostic risk model can reflect the status of tumor immune microenvironment in BLCA patients, we evaluated the correlation between risk scores and immune cell infiltrations estimated by the CIBERSORT algorithm (Figure 11). The correlations between the risk score of the prognostic signature and immune cell infiltration are shown in Figure 12. As the risk score increased, the numbers of neutrophils, M2 macrophages, and CD8+ T cells in BLCA tissues increased (Figure 12).

Summary of the 22 immune cell subtypes estimated by the CIBERSORT algorithm. The bar charts exhibit the cell proportions of BLCA patients and various colors represent the 22 immune cells with annotations below the legend.

Figure 11. Summary of the 22 immune cell subtypes estimated by the CIBERSORT algorithm. The bar charts exhibit the cell proportions of BLCA patients and various colors represent the 22 immune cells with annotations below the legend.

The correlation between the risk score and immune cell infiltration in the cohort. (A) B cells memory. (B) Dendritic cells activated. (C) Dendritic cells resting. (D) Eosinophils. (E) Macrophages M0. (F) Macrophages M1. (G) Macrophage M2. (H) Mast cell activated. (I) Mast cell resting. (J) Monocytes. (K) Neutrophils. (L) NK cells activated. (M) NK cells resting. (N) Plasma cells. (O) T cells CD4 memory activated. (P) T cells CD4 memory resting. (Q) T cells CD4 naïve. (R) T cells CD8. (S) T cells follicular helper. (T) T cells gamma delta. (U) T cells regulatory (Tregs).

Figure 12. The correlation between the risk score and immune cell infiltration in the cohort. (A) B cells memory. (B) Dendritic cells activated. (C) Dendritic cells resting. (D) Eosinophils. (E) Macrophages M0. (F) Macrophages M1. (G) Macrophage M2. (H) Mast cell activated. (I) Mast cell resting. (J) Monocytes. (K) Neutrophils. (L) NK cells activated. (M) NK cells resting. (N) Plasma cells. (O) T cells CD4 memory activated. (P) T cells CD4 memory resting. (Q) T cells CD4 naïve. (R) T cells CD8. (S) T cells follicular helper. (T) T cells gamma delta. (U) T cells regulatory (Tregs).

Discussion

Although cancer immunotherapy has expanded the treatment possibilities for BLCA, only a subset of patients responds to immunotherapy [11–14]. Thus, it is crucial to identify immune-related biomarkers for the progression of BLCA patients to improve the effect of immunotherapy. Recent studies have reported genome-wide profiling investigating the role of multiple immune-related signatures in predicting tumor outcomes [15–17]. However, very few of these studies gained constructive therapeutic implications. Here, we established a robust immune-related risk signature of BLCA by integrated analysis of transcriptional profiles in TCGA. We also conducted external validation on overall survival rate and cancer-specific survival rate through two GEO datasets (GSE13507, n = 256; GSE32894, n = 308). Moreover, we performed western blotting in bladder tissues and cell lines to validate the IRGs at the protein level. As expected, the western blot demonstrated that the protein levels of the IRGs matched their RNA-seq levels derived from TCGA and GEO datasets.

Here we used univariate Cox regression analysis to assess the relationships between 260 DEIRGs and the prognosis of BLCA patients. We found that the expression of 51 DEIRGs significantly correlated with overall survival. To explore the underlying modulators related to these IRGs, a TF regulatory network comprising 77 differentially expressed TFs and DEIRGs that potentially regulate hub IRGs were constructed. Our results reveled that STAT1, TAP1, TAP2 and CXCL10 played central roles in this network, suggesting they are hub gene regulators. The transporter associated with antigen processing (TAP) is necessary for T-cell recognition [18]. Recent studies have demonstrated that TAP plays a critical role in cell differentiation, proliferation, the development and progression of cancer [19–21]. In addition, TAP detection has been recognized as a new independent indicator for the course of chemotherapy and clinical monitoring of several type of cancers [22, 23]. STAT1 is considered as an oncogene that promotes cell adhesion, migration, and proliferation in bladder cancer [24]. Furthermore, the TF–IRG regulatory network emphasized that these hub TFs are closely associated with the prognostic signature genes, such as ADIPOQ, PDGFRA, MMP9 and RAC3.

We then determined the risk-stratification value of these prognosis-related DEIRGs and identified 11 prognostic IRGs as potential prognostic indicators in clinical practice. Furthermore, our results revealed that the prognostic model was independent of clinical characteristics. Moreover, the signature illustrated the ability to discriminate the overall survival and cancer-specific survival in both the GSE13507 and GSE32894 datasets.

Our IRG-based signature highlighted eleven IRGs, ADIPOQ, PPY, NAMPT, TAP1, AHNAK, OLR1, PDGFRA, IL34, MMP9, RAC3, and SH3BP2. Several studies showed that upregulation of MMP9 could promote the proliferation and invasion of bladder cancer [25, 26]. ADIPOQ has been proposed to be a mediator of obesity-associated metabolism and to have direct effects on the development and progression of various types of malignancies [27, 28]. Recent studies have identified NAMPT as a potential bladder cancer biomarker [29]. The subcellular localization of AHNAK exhibited different between BLCA tissues and normal tissues [30]. LOX-1 was upregulated in 57% of bladder cancer cells, and was associated with tumor invasion and metastasis [31, 32]. Furthermore, expression of PDGFRA has been reported in BLCA specimens [33, 34]. Moreover, upregulation of RAC3 in bladder cancer predicted an adverse clinical outcome and increased tumor immune response [35, 36].

Gene functional annotation indicated that our IRGs were involved in cell migration, cell proliferation, cytokine interactions, and chemokine pathways. Cytokines and chemokines, the gene products of transcription factors, played important roles in BLCA progression, and metastasis [37–39]. In addition, a PPI network was conducted to elucidate the regulatory mechanisms influence on the IRGs at the protein level. MMP9, NAMPT, CXCL12, ADIPOQ, CXCL10, TAP1, TAP2, and STAT1 figured prominently in the PPI network. Functional annotation of the core genes derived from the PPI network also showed these hub genes were mainly enriched in sprouting angiogenesis, positive regulation of leukocyte chemotaxis, and regulation of smooth muscle cell proliferation.

Our results showed that after adjusting other clinical characteristics including age, clinical stage, T and M, the risk score could be used as an independent predictor. Multivariate analysis revealed that the risk score was independently associated with the overall survival in the cohort with a considerable hazard ratio. To determine the significance of the predictive values of the DEIRGs, we evaluated the correlations between the signatures and the clinicopathological factors age, sex, pathological grade, clinical stage, and tumor stage. We found that the levels of PPY, NAMPT, ADIPOQ, IL34, TAP1, RAC3, PDGFRA and AHNAK, combined with the risk scores, positively correlated with the progression of BLCA. Thus, combining this signature with other clinical factors may serve as a tool for predicting the prognosis of BLCA patients.

Accumulating evidence indicates that immune infiltration plays vital roles in the prognosis of BLCA patients [40, 41]. Immune infiltration is an important determinant of treatment response and prognosis in BLCA patients, which is further supported by the findings of our present statistical analyses showing that the risk score positively correlated with the infiltration of tumors with CD8+ T cells, M2 macrophages, and neutrophils. The results establish the reliability of our signature to predict the prognoses of patients with BLCA. Limitations of our study include the intrinsically limited information acquired from bioinformatics analysis of transcriptome data, which may not reflect the entirety of the pathologically significant aspects of the antitumor immune response. In addition, as a retrospective study, our results still have a bias because of their heterogeneity, and so further preclinical and clinical investigations are required to identify the specific mechanisms of the effects of IRGs on BLCA.

In conclusion, we identified and validated an 11-IRG-based risk signature. The IRGs were mainly involved in tumor-related functions and pathways. Our immunogenomic signature may be a reliable tool to evaluate the prognosis of BLCA patients and help guide individualized immunotherapies. Nonetheless, further experiments are required to verify our present findings.

Materials and Methods

Clinical samples and data collection

Clinical and transcriptome RNA-seq data of 433 BLCA samples, including 414 BLCA patients and 19 matched normal samples, downloaded from the TCGA cohort. IRGs downloaded from the ImmPort portal database included 2498 IRGs [42]. ImmPort is a curated dataset to promote the reuse of immunological research data generated by intramural and extramural programs of the United States National Institutes of Health, and privately funded investigators [43].

Analysis of DEGs

DEGs were identified via the effective and convenient limma R package. We analyzed differential gene expression of all transcriptional data and IRGs using the cut-off values as follows: FDR < 0.05 and log2 |FC| > 1. The pheatmap R software package was used to display the DEGs and DEIRGs. Subsequently, prognosis-associated DEIRGs were analyzed using univariate and multivariate Cox analysis.

Functional annotations and PPI network

GO enrichment analysis and KEGG pathway analysis were conducted with the clusterProfiler package [44]. In addition, we used STRING database to predict the PPI network of the prognosis-associated DEIRGs and to evaluate the degree of interactions between proteins [39]. Cytoscape (version 3.5) was utilized to visualize the interactive network data [45]. Cytoscape plug-ins including MCODE, ClueGO, and CluePedia were used as previously described [46–48]. Specifically, MCODE was used to screen the most significant modules of hub genes from the PPI networks with selection criteria as follows: node score cut-off = 0.2, degree cut-off ≥ 2, and k-score = 2. The GO and KEGG pathway analyses of the selected hub genes were visualized utilizing ClueGO and CluePedia plug-ins.

Construction of a regulatory network linking TFs and IRGs

Cistrome Cancer is a comprehensive resource for predicting targets of TFs (http://cistrome.org/CistromeCancer/). Target prediction is based on integration of correlations of expression levels among samples in each cancer included in TCGA, as well as the genomic TF binding patterns included in the ChIP-seq data. The Cistrome Cancer database contains a list of 318 TFs [49]. We analyzed differentially expressed TFs and constructed a regulatory network linking TFs and IRGs.

Construction of a prognostic risk model

Gene-weighted values were calculated using the regression coefficients of the multivariate Cox regression model. The equation used in this analysis was as follows:

Risk score =∑i=1nCoe (gene i)∗Exp (gene i)

“Coe (gene i)” represents the regression coefficient of gene i estimated from the multivariate Cox analysis, and “Exp (gene i)” is the expression of gene i. The prognostic risk model was used to calculate the risk score of each patient. Furthermore, KM curve was performed using the “survival” and “survminer” R packages.

External validation of the proposed signature in the GSE13507 and GSE32894 cohorts

We utilized the same risk score formula and cut-off value in two external validation cohorts, GSE13507 and GSE32894, to validate the IRG-based prognostic risk model. The prognostic model was presented in each dataset that contained the differentially expressed genes, a risk plot, and the distribution of risk score. Additionally, we evaluated the area under the ROC curve with the “survival ROC” package to assess the survival differences in the external datasets.

Validation of IRGs in the risk model by western blot

Bladder cancer tissues and normal adjacent tissues were obtained from six patients admitted to Shandong Provincial Hospital Affiliated to Shandong First Medical University. T24 cells were cultured in RPMI 1640, and SV-HUC-1 cells were maintained in F-12K medium. The medium was refreshed every other day. Bladder tissues and urothelial cell lines were lysed with RIPA buffer. 25 μg of protein quantified by BCA kit in the samples were subjected to 6–10% SDS-PAGE and transferred to a polyvinylidene fluoride membrane. The membrane was blocked with 5% skim milk and incubated with primary antibodies at 4° C overnight (Supplementary Table 3). After hybridization of secondary antibodies, the protein expression level was detected by the chemiluminescence method (Amersham Imager 600, GE, USA).

Inference of immune cell infiltration in BLCA tissues

Here we estimated immune infiltration between high and low risk score groups with the LM22 (22 types of immune cells) signature file via the CIBERSORT algorithm. After 1000 permutations of CIBERSORT, the distribution of 22 subtypes of immune cells in each patient was exhibited, along with the p-values, correlation coefficients, and root mean squared error (RMSE), which evaluates the accuracy of the results of each sample.

Statistical analysis

Statistical analysis was conducted using R software, and p < 0.05 indicated a significant difference. To evaluate the accuracy of the prognostic risk model, the “survivalROC” package was used to calculate the AUC. The t test was used to evaluate continuous variables, and the χ2 test was used to compare categorical variables. Multivariate analysis was performed to determine the independent prognostic significance of the immune-related risk signature. The relationships between the risk score and the clinical characteristics were evaluated using the “beeswarm” package.

Abbreviations

BLCA: bladder cancer; IRGs: immune-related genes; DEIRGs: differentially expressed immune-related genes; NMIBC: non muscle invasive bladder cancer; MIBC: muscle-invasive bladder cancer; TCGA: The Cancer Genome Atlas; FC: fold-change; TF: transcription factor; AUC: area under the curve; ROC: receiver operating characteristic; KM: Kaplan Meier.

Author Contributions

Conceptualization, Liqiang Guo; Methodology, Shengtian Zhao; Data curation, Zhen Ma; Writing Original Draft Preparation, Qiong Wu; Writing and Editing, Liqiang Guo; Project Administration, Mingzhen Yuan; Funding Acquisition, Shengtian Zhao. All authors have read and agreed to the submitted version of the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest.

Funding

This study was funded by Clinical Research Center of Shandong University (2020SDUCRCC001), Shandong Province Natural Science Fund (ZR2020QH067), the Jinan Clinical Medical Technology Innovation Program (201907063), National Natural Science Funds of China (81670686, 81970578), and Engineering laboratory of urinary organ and functional reconstruction of Shandong Province.

Editorial Note

&

This corresponding author has a verified history of publications using a personal email address for correspondence

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. Babjuk M, Böhle A, Burger M, Capoun O, Cohen D, Compérat EM, Hernández V, Kaasinen E, Palou J, Rouprêt M, van Rhijn BW, Shariat SF, Soukup V, et al. EAU Guidelines on Non-Muscle-invasive Urothelial Carcinoma of the Bladder: Update 2016. Eur Urol. 2017; 71:447–61. https://doi.org/10.1016/j.eururo.2016.05.041 [PubMed]
  • 3. 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. 2017; 171:540–56.e25. https://doi.org/10.1016/j.cell.2017.09.007 [PubMed]
  • 4. 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]
  • 5. von der Maase H, Hansen SW, Roberts JT, Dogliotti L, Oliver T, Moore MJ, Bodrogi I, Albers P, Knuth A, Lippert CM, Kerbrat P, Sanchez Rovira P, Wersall P, et al. Gemcitabine and cisplatin versus methotrexate, vinblastine, doxorubicin, and cisplatin in advanced or metastatic bladder cancer: results of a large, randomized, multinational, multicenter, phase III study. J Clin Oncol. 2000; 18:3068–77. https://doi.org/10.1200/JCO.2000.18.17.3068 [PubMed]
  • 6. Aggen DH, Drake CG. Biomarkers for immunotherapy in bladder cancer: a moving target. J Immunother Cancer. 2017; 5:94. https://doi.org/10.1186/s40425-017-0299-1 [PubMed]
  • 7. Poli G, Egidi MG, Cochetti G, Brancorsini S, Mearini E. Relationship between cellular and exosomal miRNAs targeting NOD-like receptors in bladder cancer: preliminary results. Minerva Urol Nefrol. 2020; 72:207–13. https://doi.org/10.23736/S0393-2249.19.03297-1 [PubMed]
  • 8. Schneider AK, Chevalier MF, Derré L. The multifaceted immune regulation of bladder cancer. Nat Rev Urol. 2019; 16:613–30. https://doi.org/10.1038/s41585-019-0226-y [PubMed]
  • 9. Balar AV, Galsky MD, Rosenberg JE, Powles T, Petrylak DP, Bellmunt J, Loriot Y, Necchi A, Hoffman-Censits J, Perez-Gracia JL, Dawson NA, van der Heijden MS, Dreicer R, et al, and IMvigor210 Study Group. Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. Lancet. 2017; 389:67–76. https://doi.org/10.1016/S0140-6736(16)32455-2 [PubMed]
  • 10. Szabados B, van Dijk N, Tang YZ, van der Heijden MS, Wimalasingham A, Gomez de Liano A, Chowdhury S, Hughes S, Rudman S, Linch M, Powles T. Response Rate to Chemotherapy After Immune Checkpoint Inhibition in Metastatic Urothelial Cancer. Eur Urol. 2018; 73:149–52. https://doi.org/10.1016/j.eururo.2017.08.022 [PubMed]
  • 11. Rosenberg JE, Hoffman-Censits J, Powles T, van der Heijden MS, Balar AV, Necchi A, Dawson N, O’Donnell PH, Balmanoukian A, Loriot Y, Srinivas S, Retz MM, Grivas P, et al. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial. Lancet. 2016; 387:1909–20. https://doi.org/10.1016/S0140-6736(16)00561-4 [PubMed]
  • 12. Balar AV, Castellano D, O’Donnell PH, Grivas P, Vuky J, Powles T, Plimack ER, Hahn NM, de Wit R, Pang L, Savage MJ, Perini RF, Keefe SM, et al. First-line pembrolizumab in cisplatin-ineligible patients with locally advanced and unresectable or metastatic urothelial cancer (KEYNOTE-052): a multicentre, single-arm, phase 2 study. Lancet Oncol. 2017; 18:1483–92. https://doi.org/10.1016/S1470-2045(17)30616-2 [PubMed]
  • 13. Sharma P, Retz M, Siefker-Radtke A, Baron A, Necchi A, Bedke J, Plimack ER, Vaena D, Grimm MO, Bracarda S, Arranz JÁ, Pal S, Ohyama C, et al. Nivolumab in metastatic urothelial carcinoma after platinum therapy (CheckMate 275): a multicentre, single-arm, phase 2 trial. Lancet Oncol. 2017; 18:312–22. https://doi.org/10.1016/S1470-2045(17)30065-7 [PubMed]
  • 14. 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]
  • 15. 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]
  • 16. 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]
  • 17. Qiu H, Hu X, He C, Yu B, Li Y, Li J. Identification and Validation of an Individualized Prognostic Signature of Bladder Cancer Based on Seven Immune Related Genes. Front Genet. 2020; 11:12. https://doi.org/10.3389/fgene.2020.00012 [PubMed]
  • 18. Beck S, Kelly A, Radley E, Khurshid F, Alderton RP, Trowsdale J. DNA sequence analysis of 66 kb of the human MHC class II region encoding a cluster of genes for antigen processing. J Mol Biol. 1992; 228:433–41. https://doi.org/10.1016/0022-2836(92)90832-5 [PubMed]
  • 19. Sun C, Deng F, Meng L, Chen G. Correlation between TAP detection and common digestive tract precancerous lesions. Oncol Lett. 2018; 15:1616–20. https://doi.org/10.3892/ol.2017.7496 [PubMed]
  • 20. Ma A, Fan D, Yan F. A study of the application of TAP combined with transvaginal ultrasound in the diagnosis of early-stage endometrial cancer. Oncol Lett. 2018; 16:5186–90. https://doi.org/10.3892/ol.2018.9250 [PubMed]
  • 21. Liu Z, Cai J, Yu Y, Fang H, Si Y, Jankee JJ, Shen M. Tumor Abnormal Protein as a Novel Biomarker in Papillary Thyroid Carcinoma. Clin Lab. 2017; 63:479–85. https://doi.org/10.7754/Clin.Lab.2016.160903 [PubMed]
  • 22. Wu XY, Huang XE. Clinical Application of Serum Tumor Abnormal Protein (TAP) in Colorectal Cancer Patients. Asian Pac J Cancer Prev. 2015; 16:3425–28. https://doi.org/10.7314/apjcp.2015.16.8.3425 [PubMed]
  • 23. Cheng Y, Chen Y, Zang G, Chen B, Yao J, Zhang W, Wang H, Yu L, He P, Zhang Y, Wu H. Increased Expression of TAP Is Predictive of Poor Prognosis in Patients with Non-Small Cell Lung Cancer. Cancer Manag Res. 2020; 12:1941–46. https://doi.org/10.2147/CMAR.S239593 [PubMed]
  • 24. Sun Y, Cheng MK, Griffiths TR, Mellon JK, Kai B, Kriajevska M, Manson MM. Inhibition of STAT signalling in bladder cancer by diindolylmethane: relevance to cell adhesion, migration and proliferation. Curr Cancer Drug Targets. 2013; 13:57–68. [PubMed]
  • 25. Qin Z, Wang Y, Tang J, Zhang L, Li R, Xue J, Han P, Wang W, Qin C, Xing Q, Yang J, Zhang W. High LINC01605 expression predicts poor prognosis and promotes tumor progression via up-regulation of MMP9 in bladder cancer. Biosci Rep. 2018; 38:BSR20180562. https://doi.org/10.1042/BSR20180562 [PubMed]
  • 26. Chang Y, Jin H, Li H, Ma J, Zheng Z, Sun B, Lyu Y, Lin M, Zhao H, Shen L, Zhang R, Wu S, Lin W, et al. MiRNA-516a promotes bladder cancer metastasis by inhibiting MMP9 protein degradation via the AKT/FOXO3A/SMURF1 axis. Clin Transl Med. 2020; 10:e263. https://doi.org/10.1002/ctm2.263 [PubMed]
  • 27. Chou SH, Tseleni-Balafouta S, Moon HS, Chamberland JP, Liu X, Kavantzas N, Mantzoros CS. Adiponectin receptor expression in human malignant tissues. Horm Cancer. 2010; 1:136–45. https://doi.org/10.1007/s12672-010-0017-7 [PubMed]
  • 28. Ohashi K, Ouchi N, Kihara S, Funahashi T, Nakamura T, Sumitsuji S, Kawamoto T, Matsumoto S, Nagaretani H, Kumada M, Okamoto Y, Nishizawa H, Kishida K, et al. Adiponectin I164T mutation is associated with the metabolic syndrome and coronary artery disease. J Am Coll Cardiol. 2004; 43:1195–200. https://doi.org/10.1016/j.jacc.2003.10.049 [PubMed]
  • 29. Zhang K, Zhou B, Zhang P, Zhang Z, Chen P, Pu Y, Song Y, Zhang L. Prognostic value of serum nicotinamide phosphoribosyltransferase in patients with bladder cancer. Croat Med J. 2014; 55:507–13. https://doi.org/10.3325/cmj.2014.55.507 [PubMed]
  • 30. 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]
  • 31. Murdocca M, De Masi C, Pucci S, Mango R, Novelli G, Di Natale C, Sangiuolo F. LOX-1 and cancer: an indissoluble liaison. Cancer Gene Ther. 2021. [Epub ahead of print]. https://doi.org/10.1038/s41417-020-00279-0 [PubMed]
  • 32. Hirsch HA, Iliopoulos D, Joshi A, Zhang Y, Jaeger SA, Bulyk M, Tsichlis PN, Shirley Liu X, Struhl K. A transcriptional signature and common gene networks link cancer with lipid metabolism and diverse human diseases. Cancer Cell. 2010; 17:348–61. https://doi.org/10.1016/j.ccr.2010.01.022 [PubMed]
  • 33. Terada T. Autopsy case of primary small cell carcinoma of the urinary bladder: KIT and PDGFRA expression and mutations. Pathol Int. 2009; 59:247–50. https://doi.org/10.1111/j.1440-1827.2009.02358.x [PubMed]
  • 34. Terada T. Urinary bladder urothelial carcinoma with expression of KIT and PDGFRA and showing diverse differentiations into plasmacytoid, clear cell, acantholytic, nested, and spindle variants, and into adenocarcinoma, signet-ring cell carcinoma, small cell carcinoma, large cell carcinoma, and pleomorphic carcinoma. Int J Clin Exp Pathol. 2013; 6:1150–56. [PubMed]
  • 35. Cheng C, Song D, Wu Y, Liu B. RAC3 Promotes Proliferation, Migration and Invasion via PYCR1/JAK/STAT Signaling in Bladder Cancer. Front Mol Biosci. 2020; 7:218. https://doi.org/10.3389/fmolb.2020.00218 [PubMed]
  • 36. Ou-Yang S, Liu JH, Wang QZ. Upregulation of RAC3 in bladder cancer predicts adverse clinical outcome and increased tumor immune response. Int J Clin Exp Pathol. 2020; 13:2937–49. [PubMed]
  • 37. Eruslanov E, Neuberger M, Daurkin I, Perrin GQ, Algood C, Dahm P, Rosser C, Vieweg J, Gilbert SM, Kusmartsev S. Circulating and tumor-infiltrating myeloid cell subsets in patients with bladder cancer. Int J Cancer. 2012; 130:1109–19. https://doi.org/10.1002/ijc.26123 [PubMed]
  • 38. Feng CC, Wang PH, Ding Q, Guan M, Zhang YF, Jiang HW, Wen H, Wu Z. Expression of pigment epithelium-derived factor and tumor necrosis factor-α is correlated in bladder tumor and is related to tumor angiogenesis. Urol Oncol. 2013; 31:241–46. https://doi.org/10.1016/j.urolonc.2010.12.001 [PubMed]
  • 39. Larsen ES, Joensen UN, Poulsen AM, Goletti D, Johansen IS. Bacillus Calmette-Guérin immunotherapy for bladder cancer: a review of immunological aspects, clinical effects and BCG infections. APMIS. 2020; 128:92–103. https://doi.org/10.1111/apm.13011 [PubMed]
  • 40. Furuya H, Chan OT, Hokutan K, Tsukikawa Y, Chee K, Kozai L, Chan KS, Dai Y, Wong RS, Rosser CJ. Prognostic Significance of Lymphocyte Infiltration and a Stromal Immunostaining of a Bladder Cancer Associated Diagnostic Panel in Urothelial Carcinoma. Diagnostics (Basel). 2019; 10:14. https://doi.org/10.3390/diagnostics10010014 [PubMed]
  • 41. Luo Y, Zeng G, Wu S. Identification of Microenvironment-Related Prognostic Genes in Bladder Cancer Based on Gene Expression Profile. Front Genet. 2019; 10:1187. https://doi.org/10.3389/fgene.2019.01187 [PubMed]
  • 42. Bhattacharya S, Dunn P, Thomas CG, Smith B, Schaefer H, Chen J, Hu Z, Zalocusky KA, Shankar RD, Shen-Orr SS, Thomson E, Wiser J, Butte AJ. ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. Sci Data. 2018; 5:180015. https://doi.org/10.1038/sdata.2018.15 [PubMed]
  • 43. Bhattacharya S, Andorf S, Gomes L, Dunn P, Schaefer H, Pontius J, Berger P, Desborough V, Smith T, Campbell J, Thomson E, Monteiro R, Guimaraes P, et al. ImmPort: disseminating data to the public for the future of immunology. Immunol Res. 2014; 58:234–39. https://doi.org/10.1007/s12026-014-8516-1 [PubMed]
  • 44. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
  • 45. Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics. 2011; 27:431–32. https://doi.org/10.1093/bioinformatics/btq675 [PubMed]
  • 46. Bandettini WP, Kellman P, Mancini C, Booker OJ, Vasu S, Leung SW, Wilson JR, Shanbhag SM, Chen MY, Arai AE. MultiContrast Delayed Enhancement (MCODE) improves detection of subendocardial myocardial infarction by late gadolinium enhancement cardiovascular magnetic resonance: a clinical validation study. J Cardiovasc Magn Reson. 2012; 14:83. https://doi.org/10.1186/1532-429X-14-83 [PubMed]
  • 47. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, Fridman WH, Pagès F, Trajanoski Z, Galon J. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009; 25:1091–93. https://doi.org/10.1093/bioinformatics/btp101 [PubMed]
  • 48. Bindea G, Galon J, Mlecnik B. CluePedia Cytoscape plugin: pathway insights using integrated experimental and in silico data. Bioinformatics. 2013; 29:661–63. https://doi.org/10.1093/bioinformatics/btt019 [PubMed]
  • 49. Mei S, Meyer CA, Zheng R, Qin Q, Wu Q, Jiang P, Li B, Shi X, Wang B, Fan J, Shih C, Brown M, Zang C, Liu XS. Cistrome Cancer: A Web Resource for Integrative Gene Regulation Modeling in Cancer. Cancer Res. 2017; 77:e19–22. https://doi.org/10.1158/0008-5472.CAN-17-0327 [PubMed]
Download PDF
Corresponding Authors

Mingzhen Yuan
yuanmingzhen2005@126.com

Shengtian Zhao
zhaoshengtian@sdu.edu.cn

Keywords
immune-related genes prognosis bladder cancer bioinformatics analysis The Cancer Genome Atlas
Paper Sections