Bioinformatics Analysis Finds Immune Gene Markers Related to the Prognosis of Bladder Cancer

Bladder cancer is one of the most common malignant tumors of the urinary system that seriously threatens the health of a population. In recent years, the application of immunotherapy has significantly changed the treatment of bladder cancer, but only some patients can benefit from the treatment with immune-checkpoint inhibitors. Many problems are unsolved in the field of bladder cancer immunotherapy, especially in the search for genes that are critical to the level of immune cell infiltration and new effective therapeutic targets. We attempted to use bioinformatics analysis to identify immune gene markers related to the prognosis of bladder cancer and to establish a prognostic signature for bladder cancer patients based on their immune gene expression profiles. We used univariate Cox proportional hazards regression analysis, the least absolute shrinkage and selection operator (LASSO) Cox regression, and multivariate Cox proportional hazards regression analysis from The Cancer Genome Atlas bladder cancer cohort (TCGA-BLCA). Fifteen genes related to prognosis were screened using the survival analysis, correlation analysis, cancer and adjacent cancer differential expression analysis, and mutation analysis. The potential biological role of these genes was determined using survival analysis and principal component analysis (PCA). The receiver operating characteristic (ROC) curve assesses the prognostic value of the predictive signature. The gene ontology (GO), Kyoto Encyclopedia of Gene and Genome (KEGG), Gene set enrichment analysis (GSEA), and other methods were used to reveal the differential gene enrichment in the signaling pathways and cellular processes of high- and low-risk groups. The single-sample GSEA (ssGSEA) method was used to quantify the infiltration levels of 24 immune cells in the tumor immune microenvironment and these immune genes were found to be closely related to the tumor immune microenvironment. In summary, we screened 15 immune genes that were closely related to bladder cancer overall survival (OS) and may be potential prognostic indicators of bladder cancer. They may have research and clinical application value in bladder cancer immunotherapy. We used 15 immune genes to construct a new immune-related gene signature that was verified and could be helpful in improving individualized prognosis of patients with bladder cancer.


INTRODUCTION
Bladder cancer is the ninth most common cancer worldwide. It has the characteristics of a difficult early diagnosis, rapid metastasis, and unsatisfactory treatment. In the past 10 years, the current treatment plan has not remarkably improved the 5-year survival rate, which needs to be addressed urgently. Additionally, finding effective biomarkers that assess and promote the diagnosis, treatment, and prognosis of bladder cancer (Dobruch et al., 2016) is important. To date, the prognosis of bladder cancer mainly depends on histopathological diagnosis and the tumor staging system. However, traditional methods are not sufficient to accurately evaluate the prognosis of bladder cancer patients and meet the needs of clinicians (Jordan and Meeks, 2019). Therefore, it is imperative to develop reliable and precise prognostic biomarkers to help clinicians optimize the treatment strategies.
The ICT is a treatment against CTLA-4, PD-1, or PD-L1. Recently, ICT has been applied to many aggressive cancers and it has changed the interventions for urinary cancers including advanced bladder cancer. The inhibition of the interactions between PD-1 and PD-L1 can restore the anti-tumor activity of the T cells and enhance the immune attack on the antigen (Ghasemzadeh et al., 2016;Kim, 2016). Bladder cancer is currently a highly immunogenic malignant tumor. In recent years, the ICT has achieved very good results in bladder cancer. In particular, the application of PD-1/PD-L1 inhibitors has greatly improved the incidence of benefit from survival for some patients. However, it is undeniable that only some of these patients can benefit from the treatment of immune checkpoint inhibitors as some patients do not respond to ICT or become resistant (Cheng et al., 2018;Felsenstein and Theodorescu, 2018;Schneider et al., 2019). Many problems remain to be solved in the bladder urothelial carcinoma (BUC) immunotherapy, especially in predicting immunotherapeutic biomarkers and finding new effective therapeutic targets.
Although several studies have proposed numerous biomarkers for predicting the efficacy of a treatment such as the expression of PD-L1, TMB, and microsatellite instability (MSI) biomarkers, most of these markers focus on the tumor invasion of the lymphocyte or TME. Disturbance in the immune response in a TME plays a decisive role in the development of bladder cancer. The constituent immune cells of a TME are an important part of the tumor tissue (Neal et al., 2018;Lemos et al., 2019). Many recent studies have shown that the effect of immune checkpoint inhibitors is affected by the tumor immune microenvironment that consists of effector CD8 + , CD4 + cells, regulatory T cells, and dendritic cells (DCs) (Barry et al., 2018;Lambrechts et al., 2018;Sun et al., 2018). Improving the clinical response to an immune checkpoint blockade will require a deeper understanding of the factors affecting the local immune balance in the TME. Therefore, to explore the regulatory mechanism of the tumor immune microenvironment and the factors influencing the immune checkpoint inhibitors, we need to look for genes that are critical in affecting the level of the immune cell infiltration. Thus, targeted research and developmental interventions are of great significance for the diagnosis and treatment of bladder cancer.
In recent years, gene expression databases have been used to mine valuable therapeutic genes, identify promising prognostic factors, and analyze the molecular mechanisms of various cancers (Wei et al., 2019). Unlike the traditional individual molecularprognostic-prediction indicators, a signature combining multiple genes can significantly improve the accuracy of prognosis prediction. Based on the importance of immune regulation in the diagnosis and treatment of bladder cancer and the general role of many immune genes in the prognosis of bladder cancer, we used single-factor Cox proportional hazards regression analysis to screen prognostic genes from the 314 immune-related genes in TCGA bladder cancer cohort (TCGA-BLCA). These genes were then subjected to the LASSO Cox regression and multi-factor Cox proportional hazard regression analyses to obtain 15 genes that would help establish the optimal risk signature. Survival analysis, correlation analysis, cancer and adjacent cancer differential expression analysis, and mutation analysis were carried out to explore the potential biological role of these genes. According to the median risk score, the patients were divided into the high-risk and low-risk groups. Survival analysis, PCA analysis, and ROC curve assessed the prognostic value of the risk scores. The GO and the KEGG databases were screened by the GSEA to explore the key signal pathway differences between the high-risk and low-risk populations. Finally, the ssGSEA method was used to quantify the infiltration levels of 24 immune cells in the tumor immune microenvironment and to explore the correlation between the 15 immune genes and the tumor immune microenvironment.
In conclusion, we screened 15 immune genes that were closely related to the overall survival (OS) of bladder cancer. They may prove to be potential prognostic indicators of bladder cancer, powerful predictors of immune checkpoint inhibitor responses, and/or new targets for immunotherapy. Simultaneously, we used 15 immune genes to construct and verify a new immune-related gene signature that could improve the individualized prognosis prediction of bladder cancer patients.

Database
The RNA-seq data and data on the clinical characteristics (including patient age, sex, stage, smoking status, and TNM staging) of the bladder cancer cohort were obtained from TCGA 1 database.

Selection of Immune-Associated Genes
The IMMUNE_RESPONSE and IMMUNE_SYSTEM_PROCESS 2 immune gene sets were obtained from the Molecular Signatures Database (MsigDB 2 ) and 314 duplicate immune-related genes were removed. The mRNA expression data of these 314 genes were obtained from TCGA-BLCA.

Identification and Validation of the Prognostic Gene Signature
The "survival" package of R language was used to perform the univariate Cox proportional hazards regression analysis to screen the immune genes that were significantly related to the OS of the TCGA-BLCA cohort. Using the "glmnet" package of R for the LASSO Cox regression analysis, a reduction in the dimensionality of high-dimensional data was achieved by limiting the sum of absolute values of the coefficients to less than a predetermined value. Variables with relatively small contributions were given coefficients of zero; only the genes with non-zero coefficients in the LASSO regression analysis were selected for further analysis. Finally, the obtained genes were subjected to multi-factor Cox proportional hazard regression analysis and screening to obtain 15 immune genes that would determine the best prediction signature are shown in Figure 3. These genes were selected to further calculate the risk score of each patient (Bao et al., 2014;Cheng et al., 2016).
According to the median value of the risk coefficient, the patients were divided into the high-risk and low-risk groups. The univariate Cox proportional hazard regression analysis and multivariate Cox proportional hazard regression analysis were performed on the risk value by using the "survival" package of the R language. Cox proportional hazard regression signature includes the risk score, age, gender, grade, T-, N-, and M-phase, and the smoking status. The Kaplan-Meier survival analysis was subsequently performed using the R "survival" package. The sensitivity and specificity of the ROC curve were used to evaluate the prognostic performance of the signature and the PCA was used to analyze the expression pattern of the grouped samples. A correlation analysis of the 15 immune genes was performed using the R "corrplot" package in the Pearson method and the results were displayed in the form of a Circos diagram are shown in Figure 4. The expression of these 15 immune genes was compared between cancer tissues and normal tissues. The mutations of the 15 immune genes in the TCGA-BLCA cohort were downloaded from the cBioPortal website 3 . The Kaplan-Meier survival analysis was performed on these 15 genes.

Pathway Analysis
The "edgeR" package calculation in R language was used to perform a differential analysis of the mRNAs of the low-risk and the high-risk groups. To perform functional annotation from the GO 4 for mRNAs with FDR values less than 0.05, the biological functions of differential genes, including biological processes (BPs), cellular components (CCs), and molecular functions (MFs) were analyzed. The KEGG 5 database analyzes the metabolic pathways and signal transduction pathways in which differential genes are significantly enriched. A GSEA 6 was then performed to reveal the signaling pathways and BPs in which differentially expressed genes were enriched between the high-risk and low-risk subgroups.

Tumor Immune Microenvironment Analysis
Inference of Infiltrating Cells in the TME Based on the immune cell marker genes provided by Bindea et al. (2013), a ssGSEA was used to quantify the infiltration levels of the 24 types of immune cells, including the T lymphocytes, DCs, and natural killer cells; the TCGA-BLCA (Finotello and Trajanoski, 2018) database was also used for this purpose. According to the level of immune cell infiltration, the patients were divided into the high-infiltration group and low-infiltration group. Heat maps were plotted to observe the relationship between risk value, age, T, N, and M stages, gender, and the immune infiltration levels of various immune cells in the high-and low-risk groups. Finally, the correlation between at-risk cells and immune cells was calculated using the Pearson method.

Statistical Analysis
All analyses were performed using the R programming language 7 . Univariate and multivariate Cox proportional hazard regression analyses were also used to assess the relationship between the risk score and OS. The ROC analysis was used to detect the sensitivity and specificity of the genetic signature risk scores to predict survival. The area under the ROC curve (AUC) was used as an indicator of prognostic accuracy. In all analyses, P-values less than 0.05 were considered statistically significant.

Construction of a Prognostic Signature for TCGA-BLCA
A total of 314 immune-related genes were obtained from the MSigDB. Univariate Cox regression analysis was performed on these genes. Seventy-six immune genes from the TCGA-BLCA database were found to be significantly related to the OS. These genes were subjected to the LASSO regression analysis to calculate the correlation coefficients. The coefficients of each gene are shown in Figures 1A,B. Twenty-nine immune genes were obtained for the multi-factor Cox regression analysis. The signature performed best when only 15 genes were included. The results of the multi-factor Cox regression analysis for 15 genes   (Figures 2A,B). The ROC analysis detects the sensitivity and specificity of the genetic risk scores in predicting survival. The AUC of riskScore was 0.751 ( Figure 2C), indicating that the signature displayed good sensitivity and specificity for predicting survival. The Kaplan-Meier survival analysis showed that the high-risk group had a lower prognosis than the low-risk group (P < 0.01, Figure 2D). The PCA results showed that our 15 immune genes could better divide the high-and low-risk patients into two groups compared with all other genes ( Figure 2E) and the immune gene set 304 genes ( Figure 2F). These results confirm the sensitivity and specificity of the signature.

Features of the Prognostic Signature
We compared the expression of the 15 immune gene cancer tissues and adjacent tissues in the TCGA-BLCA cohort ( Figure 5A) and found that the expression of the CCR9, IL7, PTGER4, IL10, CTSG, and ZBTB16 proteins in the adjacent tissues was significantly higher than that in the cancerous tissues (P < 0.05). The expression of CEBPG and RUNX1 proteins in the cancer tissues was significantly higher than that in the adjacent tissues (P < 0.05). Next, we examined the mutations in these genes in bladder cancer ( Figure 5B) and found that FIGURE 3 | The patients were divided into two groups: low-risk and high-risk. As the risk score increased, the survival time of patients decreased and the number of deaths increased. The heat map shows the expression profile of the 15 immune genes in the prognostic markers.
the mutation frequency of PTGER4 was the highest, reaching 9%, with amplification mutations as the main mutation; it was followed by RUNX1, IL7, and CIITA, with mutation frequencies of 6, 5, and 5%, respectively. By analyzing the relationship between the 15 immune genes and the OS, it was found that the group showing a higher expression of the CDK6, HDAC7, CTSG, EREG, and ZBTB16 mRNAs had a lower survival time and was smaller. The group showing a higher expression of the ZAP70, IL7, and PTGER4 proteins had a significantly longer survival time than the lower expression group (P < 0.05, Figure 6).

Identification of the Involved Signaling Pathways
The signaling pathway enrichment analysis of differential mRNA (FDR < 0.05) in the low-risk and high-risk groups and the GO analysis showed that the differential genes were related to the extracellular matrix organization, extracellular structure organization, and extracellular matrix. The collagen-containing extra cellular matrix, extra cellular matrix structural constituents, integral binding, and other BPs are closely related (Figures 7A-C). The KEGG analysis revealed that these differential genes were mainly enriched in the PI3K-Akt signaling pathway, in the proteoglycans in cancer, human papillomavirus infection, and other signaling pathways ( Figure 7D and Table 2). The GSEA analysis showed that the signaling pathways such as E2F targets, hypoxia, G2/M DNA damage checkpoint, apical junction complex, epithelial-mesenchymal transition (EMT), KRAS Signaling UP, mTORC1 signaling, mitotic spindle, and TNFα signaling via NF-Kβ were significantly increased in the high-risk group (Figure 7E).

The Prognostic Signature Is Related to the Tumor Immune Microenvironment
We quantified 24 types of immune cells including the B cells, T cells, natural killer cells, macrophages, DCs, and myeloid subpopulations to investigate the composition of the tumor immune microenvironment and draw a heat map to observe the risk values, age, stages T, N, and M, gender, and immune infiltration (Figure 8). At the same time, we found cytotoxic cells, DC, eosinophils, CD56 bright NK cells, CD8 T cells, T cells, and T helper cells in the high-risk group. The level of infiltration of TFH and Th17 cells was significantly lower than that in the low-risk group, while the levels of macrophages, neutrophils, CD56 dim NK cells, NK cells, Th1 cells, and Th2 cells in the high-risk group were significantly higher than that in the lowrisk group (Figure 9A). At the same time, we analyzed the correlation between the immune cells. Among them, we focused on the significant positive correlation between CD8 + T cells closely related to the immune checkpoint inhibitors and iDC, DC, pDC, TReg, T cells, and cytotoxic cells ( Figure 9B) Figure 9C).

DISCUSSION
Bladder cancer is the ninth most common cancer in the world, affecting 430,000 people and causing 165,000 deaths each year. Although considerable time, effort, and expense have been invested in bladder cancer research, the overall morbidity and mortality has not improved significantly over the past 20 years (Antoni et al., 2017). Modern treatments for bladder cancer include surgery, chemotherapy, radiation therapy, and immunotherapy. Immunotherapy brings a new hope in the treatment of bladder cancer. A large number of basic and clinical experiments are currently underway, including allogeneic stem cell transplantation, antitumor vaccines, proinflammatory cytokines, chimeric antigen receptors, and adoptive T cell metastasis. One of the most promising methods considered is ICT (Yu et al., 2018;Zhu et al., 2019;Zhuang et al., 2020).
In the human immune cycle of tumors, the antigens produced by the tumor cells are captured by the DCs. The major histocompatibility complex (MHC)-I and MHC-II on the surface of the DCs present these antigens to the T cells for recognition and lead to the activation of effector T cells (Ramachandran et al., 2017;Pio et al., 2019;Ryan et al., 2019). The internal and external environment of the tumor cells during tumorigenesis and metastasis is called the TME. The TME contains tumor cells and surrounding immune cells, endothelial cells, fibroblasts, extracellular matrix, secreted cytokines, chemokines, etc. Tumors can create a series of favorable conditions for themselves through the TME and even escape the immune cycle. In cancer patients, the tumor's immune cycle does not perform well. The most important reason is that in the TME, there are some inhibitory signals that suppress the immune function of the effector T cells. The main role of the T cells is to distinguish healthy cells from pathogens or malignant cells by activating or deactivating various receptors on the surface of the T cells. The inhibitory signals include a class of signaling pathways called the immune checkpoints. These immune checkpoints are normally used to maintain the body's autoimmune tolerance and prevent these killer T cells from attacking their own cells because these molecules and their relevant receptors on the T cells "control" the immune system by blocking the immune function, so they are collectively referred to as checkpoint proteins. In order to escape the hunting by T cells, some tumor cells also generate some inhibitory signals on their surface. The immune function of the T cells is suppressed by the immune checkpoints. As a result, the immune system remains inactive against malignant cells, allowing their uncontrolled growth and proliferation. The immune checkpoint inhibitors interfere with the action of these checkpoint proteins to prevent tumors from suppressing the T cells and restart the tumor immune cycle (Kieler et al., 2018;Li et al., 2019). Immune checkpoint inhibitors, in principle, should have a wide range of killing capabilities against cancer, but a large number of cases of non-response have been found in clinical applications (Ellmark et al., 2015; Coffelt and de Visser, 2016). Recent research indicates that non-response of immune checkpoint inhibitors may be related to various factors such as tumor immune microenvironment regulation, single immune checkpoint inhibitor suppression, and blocked T cell infiltration during the immune cycle. Therefore, the subsequent objectives include exploring the regulatory mechanism of the TME, clarifying the mechanism of immune checkpoint inhibitor nonresponse, and finding promising new targets for immunotherapy.
Gene markers, often used to predict prognosis, have been reported to be more accurate than the TNM staging methods in multiple cancer species (Bao et al., 2014;Liu et al., 2019). In this study, we screened genes related to prognosis from 314 immune-related genes. Using univariate Cox proportional hazards regression analysis, Lasso regression analysis, and multivariate Cox proportional hazards regression analysis, we finally obtained 15 independent prognostic immune genes. The survival analysis, correlation analysis, cancer and adjacent-cancer differential expression analysis, and mutation frequency analysis were performed on these genes. Some of these genes were closely related to the occurrence and development of bladder cancer: CDK6, IL-10, and RUNX1, of which CDK4/6 inhibitors are a promising treatment strategy for the treatment of bladder cancer (Zhao et al., 2015;Sun et al., 2019;Tong et al., 2019). The primary bladder tumor cells secrete a large amount of IL-10. Removing IL-10 in co-cultures of monocytes and tumor cells can reduce the upregulation of PD-L1 in monocytes and affect the curative effect of the immune checkpoint inhibitors . RUNX1 is a novel direct target of miR-27a, which is involved in the regulation of sensitivity to bladder cancer chemotherapy (Deng et al., 2015). There are few reports on the role of IL-7  and HDAC7 in bladder cancer. It was found that IL7 might be involved in T cell activation and play a role in the anti-CTLA-4 immunotherapy. Niegisch et al. (2013) showed that in urothelial carcinoma, the upregulation of the mRNAs of HDAC2 and HDAC8 and the downregulation of the mRNAs of HDAC4, HDAC5, and HDAC7 are common findings. The role of genes other than those mentioned here has not been reported in bladder cancer: CCR9, ZAP70, PTGER4, CTSG, CEBPG, PF4, MAP3K7, ZBTB16, CIITA, and EREG. These genes may be new markers for predicting the prognosis of bladder cancer and new targets   for immunotherapy, but further basic and clinical laboratory identification is needed. Signal pathway enrichment analysis of the differential mRNA (FDR < 0.05) in the low-risk and high-risk groups and the GO analysis results consisted of the CC, BP, and MF. At the CC level, the differential genes are related to the inheritance junction, extracellular matrix, and plasma membrane protein complex. At the MF level, the differential genes are related to promoter promoter-specific DNA binding, RNA polymerase II proximal promoter sequence-specific DNA binding, etc. At the BP level, the differential genes are related to the regulation of body fluid levels, cell morphogenesis involved in neuron differentiation, etc. These pathways are closely related to the BPs of tumorigenesis and development. KEGG analysis found that these differential genes were mainly enriched in the classical cancer signaling pathways such as the PI3K-Akt signaling pathway, proteoglycans in cancer, human papillomavirus infection, and focal adhesion ( Figure 7D and Table 2).
For example, the PI3K-Akt signaling pathway has been studied to confirm that when its function is normal, the PI3K pathway regulates key cellular functions, including cell growth, movement, proliferation, and differentiation. However, excessive activation of the PI3K signaling pathway causes breast cancer and ovarian cancer (Shukla et al., 2007). New research shows that the PI3K-Akt signaling pathway plays an important role in the tumor immune microenvironment, affecting the efficacy of the immune checkpoint inhibitors Cristofoletti et al., 2019). The GSEA analysis found that the proteins in the high-risk group  Figure 7E and Table 3). Most of these pathways are closely related to tumorigenesis and the regulation of the TME. For example, hypoxia is a common feature of malignant tumors. It can regulate the tumor immune microenvironment by regulating a variety of immune cells. Hypoxia significantly reduces the T lymphocyte proliferation and activation, decreases the NKG2D receptor on the NK cells, and thereby inhibits the killing function of the NK cells, increases tumor-associated macrophages to induce angiogenesis, and reduces inflammation to promote tumor progression (Marques et al., 2012;Tian et al., 2017;Jiao et al., 2019). The EMT signaling pathway is an important BP for the epithelial-derived malignant tumor cells to acquire the ability to migrate and invade. It is of great significance in the occurrence, development, and metastasis of bladder cancer and participates in the TME regulation (Zhaojie et al., 2019). Interferons play a vital role in the immune response of the body toward malignant cells. Type I interferons (IFNα and IFNβ) directly regulate the transcription of more than 100 downstream genes, resulting in countless direct (via cancer cells) and indirect (via immune effector cells and vasculature) on tumors. The IFN-α/β receptor (IFNAR) signaling can promote the Treg function in autoimmunity. Activation of the IFNα signaling pathway leads to a more effective antiviral response and enhanced antitumor immunity (Gangaplara et al., 2018;Borden, 2019). In our study, the high-risk group was negatively correlated with the interferon alpha response signal pathway (NES = −1.427, P.adjust = 0.041, Table 3). Finally, we quantified the levels of 24 immune cell infiltrations in TCGA-BLCA tumor samples using the ssGSEA method. On comparing these levels in the high-and low-risk groups, we found that the levels of the cytotoxic cells, DC, CD8 + T cells, T cells, the T helper cells, TFH, Th17 cells, CD56 bright NK cells, and eosinophils were significantly lower than that in the low-risk group. It is well known that the higher the infiltration levels of the cytotoxic cells, DCs, CD8 + T cells, T cells, and T helper cells in a patient's tumors, the greater the survival benefit for the patients. The cell infiltration level groups such as DC, CD8 + T cells, T cells, and T helper cells affect the efficacy of ICT in a positive manner (Jiao et al., 2019). This may partly explain the phenomenon that the survival time of the high-risk group is significantly lower than that of the low-risk group ( Figure 9A). The above results further suggest the reliability of the prediction signature, its relevance to the TME, immunological examination, and the importance of the 15 immune genes. We then analyzed the correlation between these 15 immune genes and immune cells. We observed that most of the immune genes have a high correlation with the level of immune cell infiltration in the TME, but the specific mechanism needs further experimental investigation.

CONCLUSION
In summary, we screened 15 immune-related markers that have independent prognostic significance for bladder cancer. They may be used as potential prognostic indicators of bladder cancer and related to the level of tumor cell microenvironment immune cell infiltration. We hope to provide an additional feasible method for assessing the prognosis of bladder cancer and may provide valuable new targets for anti-tumor immunotherapy.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. These can be found in The Cancer Genome Atlas (https://portal.gdc.cancer. gov/).