Development of prognostic signature based on immune-related genes in muscle-invasive bladder cancer: bioinformatics analysis of TCGA database

Background: Muscle-invasive bladder cancer (MIBC) with high tumor stages accounts for most bladder cancer patient mortality. Platinum-based chemotherapy provides insufficient survival benefits; however, immunotherapy is a promising option for MIBC. Results: There were 31 differentially expressed IRGs that significantly correlated with the clinical outcomes of MIBC patients. A prognostic signature based on 12 IRGs (MMP9, RBP7, ADIPOQ, AHNAK, OAS1, RAC3, SLIT2, EDNRA, IL34, PDGFD, PPY, IL17RD) performed moderately in prognostic predictions with area under the curve (AUC) equal to 0.76. The high-risk patient group presented worse survival outcomes (hazard ratio 1.197, 95% confidence interval 1.103–1.299, p < 0.001). Furthermore, immune cell infiltration analysis showed increased tumor infiltration of macrophages in the high-risk group. Conclusion: This novel prognostic signature can effectively divide MIBC patients into different risk groups, allowing for intensive treatment of high-risk individuals who have worse predicted survival outcomes. Methods: Bioinformatics analyses were conducted using the Cancer Genome Atlas (TCGA) database. Differentially expressed genes and survival-associated immune-related genes (IRGs) were analyzed through a computational algorithm and Cox regression. The potential mechanisms of IRG expression were explored with transcription factors, and a prognosis classification based on IRG expression was developed to stratify patients into distinct risk groups.

AGING progression using specific and sensitive biomarkers could reduce the likelihood of advanced disease.
Immunotherapy has been a major driver of personalized medicine and has demonstrated aggressive anti-tumor effects by stimulating the immune system [5,6]. Since immune components in the tumor microenvironment play an important role in tumor cell gene expression and clinical outcomes [7][8][9], measurement of these immune components provides a method of predicting the longterm prognosis of cancer patients. The differential expression of immune-related genes (IRGs) is considered an effective biomarker in many types of cancer.
This study explores IRG expression's potential clinical utility for prognostic stratification and targeted MIBC immunotherapy. We analyzed IRG expression with corresponding clinical information to develop an individualized prognostic model for MIBC patients. Bioinformatics analyses were conducted to explore the underlying regulatory mechanisms of IRGs. Our results provide a foundation for developing personalized treatment for MIBC patients.

Identification of differentially expressed IRGs
The edgeR algorithm identified 4,876 differentially expressed genes, of which 3,453 were upregulated and 1,423 downregulated ( Figure 1A, 1C). From this set of genes, we extracted 260 differentially expressed IRGs, including 120 upregulated and 140 downregulated ( Figure 1B, 1D). The list of differentially expressed genes and IRGs is presented in Supplementary Tables 1, 2, respectively.

AGING
According to the GO enrichment analysis, inflammatory pathways were the most frequently enriched. "T cell activation," "collagen-containing extracellular matrix," and "receptor-ligand activity" were the most common terms among biological processes, cellular components, and molecular functions, respectively ( Figure 2A). As for KEGG pathways, cytokine-cytokine receptor interactions were most often enriched by differentially expressed IRGs ( Figure 2B).

Identification of survival-associated IRGs
MIBC-specific IRGs related to survival outcomes (p < 0.01) were selected for the survival analysis to make the prognosis index more precise (Table 1). There were 31 IRGs that met this criterion, four of which were risk factors and the others protective factors.

TF regulatory network
The regulatory relationships of the survival-associated IRGs and TFs were analyzed to explore the potential mechanisms corresponding to clinical significance ( Figure  3). We examined the expression profiles of 318 TFs and found 77 that were expressed differently in MIBC than in normal tissues. Setting the correlation score cutoff value to > 0.5 resulted in 5 TFs that were screened to form the regulatory network in association with the IRGs.

Evaluation of clinical outcomes
Based on the LASSO regression analysis of the expression of IRGs, the patients were divided into two risk groups (Figure 4). The formula for the prognostic index was as follows: With the AUC of the ROC curve equal to 0.760, the prognostic index showed moderate potential for prognosis prediction ( Figure 5A). According to the results from the multivariate Cox regression, the prognostic risk score could act as an independent factor after adjusting for clinical parameters (hazard ratio 1.197, 95% confidence interval 1.103-1.299, p < 0.001] ( Figure  5B). Overall survival of the low-risk group was significantly greater than the high-risk group (p < 0.001) based on the Kaplan-Meier survival curves ( Figure 5C).
To further explore the impact of IRGs on clinical outcomes, we conducted multivariate regression of the IRGs to elucidate the relationships between IRG expression and the clinicopathological factors in MIBC ( Table 2). Results showed patients in the high risk group presented higher T stage, M stage and clinical stage.
The immune cell infiltration analysis was performed to identify changes in the tumor microenvironment   AGING produced by the differential expression of IRGs ( Figure 6). Among the six types of immune cells studied (B cells, CD4 + T cells, CD8 + T cells, macrophages, neutrophils, and dendritic cells), only macrophages showed a significant difference in infiltration abundance between the two risk groups (β = 0.18, p < 0.001).

DISCUSSION
According to our bioinformatics analysis, 31 differentially expressed IRGs significantly correlated with the clinical outcome of MIBC patients. A prognostic signature based on 12 IRGs (MMP9, RBP7, ADIPOQ, AGING AHNAK, OAS1, RAC3, SLIT2, EDNRA, IL34, PDGFD, PPY, IL17RD) performed moderately in prognostic predictions with AUC equal to 0.76. Patients in the high-risk group experienced worse survival outcomes compared with those in the low-risk group. Furthermore, immune cell infiltration analysis showed increased tumor infiltration of macrophages in the highrisk group, suggesting that an abnormal immunoreaction may have led to the disparate prognoses of this group.
Immune regulation plays an essential role in bladder cancer through intrinsic and extrinsic control of immune cell activation [10]. Tumor cells can exacerbate immunosuppressive pathways, fostering a tolerant microenvironment. The tumor microenvironment comprises many different immune cell subpopulations endowed with anti-tumor or pro-tumor activity [11]. Furthermore, different levels of tumor-infiltrating immune cells in the tumor microenvironment affect the prognosis of patients [12]. As for MIBC, molecular biomarkers and histopathology are predictive tools to assess possible clinical outcomes and guide further therapies [13,14]. Based on these findings, IRG signatures could be used to identify potentially high-risk patients and thus allow for individualized treatment strategies.
Many of the 12 IRGs in our prognostic signature have previously been associated with cancer survival outcomes. Researchers have demonstrated that MMP9 is significantly associated with poor survival of cancer patients by way of promoting the invasiveness of cancer cells [15], and high expression of MMP9 has been observed in bladder cancer patients with poor prognosis and rapid tumor progression [16]. According to Elmasry et al., overexpression of RBP7 was an independent biomarker of poor cancer-specific survival in colon cancer [17], although no studies of RBP7 in MIBC have been reported. ADIPOQ is relevant in obesity and may potentially lead to metabolic syndrome [18], thus possibly worsening the prognosis of MIBC patients. As for AHNAK, two studies have conducted prognosis model analyses with this and other bladder cancer genes [19,20], and both showed poor outcomes with high levels of AHNAK expression. One study has reported that the methylation status of SLIT2 is associated with low stage and histological grade for the initial diagnosis of non-MIBC [21]. Still, the gene expression level was unmentioned in the study. EDNRA overexpression was shown to be associated with metastasis and poor outcome in advanced bladder cancer patients [22]. Further research is needed to elucidate the impact of these IRGs on MIBC outcomes more fully.
Although several previous studies are focusing on distinct survival outcomes for differentially expressed genes [23][24][25], no genome-wide profiling studies of MIBC have been conducted. Until recently, the  Platinum-based chemotherapy is still widely used for metastatic bladder cancer without significant survival improvement [26]. However, in recent years, cancer immunotherapy has played an important role in clinical cancer management [27]. With the discovery of immune checkpoint inhibitors, MIBC patients may receive better prognoses. Programmed cell death 1 (PD-1, also known as PDCD1) on the surface of CD8 + T cells binds to programmed cell death ligand 1 (PD-L1, also known as CD274) produced by tumor tissue, leading to a limited host immune response. By increasing the level of AGING infiltrated CD8 + T cells, PD-L1 inhibitors demonstrate an effective anti-tumor immune response [28]. T cells play an important role in the "immune surveillance" theories of cancer, with complex effects, including expression of polymorphic antigen receptors for specific antigen recognition, possession of effector functions, and development of memory characteristics [29,30]. However, studies have only observed a significant increase in the level of CD4 + cells (specifically Th17 cells) in the blood of bladder cancer patients [31]; thus, the role of tumor-infiltrated CD4 + cells has not been elucidated. Our results showed no distinct expression differences of CD4 + cells in the two risk groups, indicating that the level of CD4 + cells probably does not independently reflect the risk classification for MIBC. Previous research demonstrated that macrophages are divided into different subtypes, with anti-tumor (M1) and pro-tumor activities (M2), but the holistic phenotype depended on the cytokine microenvironment in the tumor tissues [32]. For MIBC, high levels of tumorinfiltrated macrophages indicate a dismal prognosis and adjuvant chemotherapy resistance [33]. Our current analyses showed an upregulated level of macrophages in the high-risk group of MIBC patients. Moreover, Zhou et al. reported tumor-infiltrating neutrophils' involvement in tumorigenesis and a negative correlation between neutrophils and CD8 + T cells [34], which indicates a state of immunosuppression in MIBC caused by a high level of tumor-infiltrating neutrophils.

AGING
Although a limitation of our study is that the effects of specific genes in tumor growth, invasion, and metastasis could not be evaluated, our analysis provided important information on MIBC survival outcomes.

CONCLUSIONS
Our prognostic signature effectively divides MIBC patients into different risk groups; thus, individuals with worse predicted survival outcomes can be identified before the advancement of their disease and treated intensively with appropriate therapies. The differential gene expression identified in our study provides insight into the MIBC tumor microenvironment, contributing to the development of IRG-targeted therapies to improve the survival of MIBC patients.

Clinical samples and data acquisition
We obtained transcriptome RNA-sequencing (RNA-seq) data of bladder cancer samples from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/ projects/TCGA-BLCA), which included data from 414 primary cancers and 19 non-tumor tissues. Raw count data and relevant clinical information for the patients were also downloaded and extracted. A list of IRGs was obtained from the Immunology Database and Analysis Portal (ImmPort) database [35]. A list of 318 transcription factors (TFs) was downloaded from the Cistrome Cancer database [36] to explore the potential mechanisms of IRG expression. Patients with MIBC (tumor stage ≥T2) were identified with this clinical data and included in the study. Patients with a follow-up duration of less than 60 days were excluded. The analysis included 375 MIBC patients. To verify the robustness of the results, we used the Gene Expression Omnibus database to perform a cross-validation, among which 165 samples were selected. The flowchart outlining the sample selection is provided in Supplementary Figure 1.

Differential gene analysis
We used the limma (linear models for microarray data) software package for the statistical R programming language to identify the differentially expressed genes between the tumor tissues and normal tissues [37], and the differential expressions were analyzed using the Wilcoxon test [38]. We set a false discovery rate <0.05 and a log2(fold change) >1 as the cutoff values. We then performed the Spearman correlation test to identify the association of IRGs and TFs, using a correlation filter of 0.5. Also, Gene Ontology (GO) enrichment was performed to investigate the functions of the gene products [39] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [40] pathway enrichment was conducted to explore the potential molecular mechanisms of the differentially expressed IRGs.

Survival analysis and tumor-infiltrating immune cell analysis
Survival-associated IRGs for MIBC were identified using univariate and multivariate Cox regression analyses with a threshold value of p < 0.05. With the use of least absolute shrinkage and selection operator (LASSO) regression analysis, we selected specific IRGs. We calculated the individualized risk score with coefficients to construct a prognostic signature separating the highrisk and low-risk groups. Subsequently, Kaplan-Meier survival curves were generated based on the different risk groups to explore the survival differences. A receiver operating characteristic (ROC) analysis was performed to validate the predictive value of our model, and the area under the curve (AUC) was calculated. The Cox regression analyses were conducted for both clinical characteristics and risk scores to investigate whether risk score could independently predict survival prognosis. The tumor-infiltrating immune cell analysis (including B cells, CD4 + T cells, CD8 + T cells, macrophages, AGING neutrophils, and dendritic cells) was performed using the TIMER database [41]. In the above analyses, p < 0.05 was considered statistically significant.

AUTHOR CONTRIBUTIONS
All authors contributed to the planning and design of the study. K.J. and S.Q. were involved in review of the raw data and directly involved in the analysis. J.K.L. and X.N.Z. provided analytical feedback based on aggregated results. D.J. and X.H.Z. drafted the manuscript, with input from all authors. X.Y.L. was responsible for the chart making. All authors provided substantive review and commentary on multiple drafts and approved the final version. K.J. and S.Q. had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.