Construction and Analysis of an Immune-Related Gene Prognostic Index for Bladder Cancer ()
1. Introduction
Bladder cancer is a common malignant tumor with the highest incidence and mortality rate among urological tumors [1] [2]. The traditional treatment of solid tumors is mainly surgical resection and chemotherapy [3] [4], but with the rise of immunotherapy, immune checkpoint inhibitors represented by PD-1 and PD-L1 have achieved very good therapeutic results in tumor treatment [5] [6].
The advent of immunotherapy has fundamentally changed the oncology treatment strategy, giving hope to some oncology patients, especially those with advanced cancer who are not sensitive to chemotherapy [7] [8]. Immunotherapy may be a better option for a small group of patients with bladder cancer [9]. Several biomarkers are available to help us better tap into these potentially beneficial patients, but their credibility needs to be validated. Clinical limitations similar to TMB and PD-L1 still exist in predicting tumor immunosuppressive site therapy [10] [11]. Also, there is limited research on the immune microenvironment (TME) of bladder cancer [12]. Therefore, establishing validated prognostic and therapeutic indicators is important for clinical decision-makers in the choice of diagnosis and treatment.
In this study, we used bladder cancer data from The Cancer Genome Atlas (TCGA) project to analyze differentially expressed genes in bladder cancer samples and normal samples, obtained immune-related differentially expressed genes after comparing them with immune-related genes in ImmPort, and screened them by Weighted Gene Co-expression Network Analysis (WGCNA) for prognosis The immune pivotal genes directly related to prognosis were screened by Weighted Gene Co-expression Network Analysis (WGCNA), and then a new validated Immune-Related Gene Prognostic Index (IRGPI) identifying bladder cancer mortality risk was constructed. Meanwhile, the molecular and immune characteristics of the constructed IRGPI were analyzed to determine its prognostic role in immunotherapy, and then compared with other tumor biomarkers, such as TIDE and TIS. Overall, IRGPI is a promising prognostic biomarker.
2. Methods
2.1. Patients and Datasets
Clinical characteristics, RNA-seq expression data, and somatic mutation information of bladder cancer patients were collected from the Tumor Genome Atlas database (https://portal.gdc.cancer.gov/). A total of 433 samples of RNA-seq expression data were collected, including 19 normal samples and 414 tumor samples. RNA sequence data and survival information of 256 bladder cancer samples (GSE13507) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The list of immune-associated genes was obtained from ImmPort (https://www.immport.org/shared/home).
2.2. Identification of Immune-Related Core Genes
A list of differentially expressed genes (p-value < 0.05, |log2FC| > 0.585) was identified by using the limma package of R based on RNA-seq data of bladder cancer samples obtained from TCGA. After considering the list of immune-related genes obtained from ImmPort in context, immune-related differentially expressed genes were obtained and analyzed by using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.
The core genes were then analyzed by WGCNA. First, a similarity matrix was constructed using the expression data by calculating the Pearson correlation coefficient between the two genes. Then, the similarity matrix was converted into an adjacency matrix with a network type of sign and a soft threshold of β = 5, and then into a Topological Overlap Metric (TOM) describing the degree of association between genes. The genes are clustered with 1-TOM as the distance, and then dynamic pruning trees are constructed for module identification. Finally, five modules were identified by setting the merging threshold function to 0.25. Based on the genes of significantly related modules (yellow module and blue-green module), the network was constructed using the edge between two genes with weight > 0.2. The top 50 genes in the network in terms of ranking are the central genes. The best cut-off value of Overall Survival (OS) was calculated for each core gene using R language, and immune-related core genes significantly associated with survival were screened for further analysis (p < 0.05, log-rank test). To reveal the associated genetic alterations, somatic mutations in 29 immune-associated HUB genes were analyzed using R.
2.3. Construction and Validation of Immune-Related Gene Prognostic Index (IRGPI)
Among the immune-related core genes that had been screened above, genes with significant effects on OS were screened and IRGPI was constructed by multifactorial Cox regression analysis. The IRGPI of each sample was obtained by multiplying the expression values of specific genes by their weights in the Cox model and then summing them. The prognostic power of the IRGPI was assessed by performing log-rank tests on the K-M survival curves of the TCGA and GEO cohorts. To verify the independent prognostic value of IRGPI, one- and multifactor Cox regression analyses were performed. Independent immune-related genes associated with prognosis were identified by univariate and multifactorial Cox regression. Then, risk score models were constructed based on the expression levels and coefficients of these genes. The risk score was calculated for each BLCA patient using the following formula: risk score = β1 * Exp1 + β2 * Exp2 + βi * Expi, where β denotes the independent prognosis-related RBP, the coefficient value proficiency represents the expression level of the prognosis-independent gene, and i represents the ith gene.
2.4. Comprehensive Analysis of Molecular Immune Characteristics of Different IRGPI Subgroups and ICI Treatment
In the signaling pathway analysis, the samples from the high and low IRGPI scoring groups were first analyzed for differential expression using R’s Limma package to identify signaling pathways involving differentially expressed genes, and then Gene Set Enrichment Analysis (GSEA) was performed on the KEGG and HARMARK gene sets using R’s clusterProfiler package (p < 0.05 and fdr < 0.25). Single-sample GSEA (SsGSEA) analysis was performed with R’s GSVA package for several representative gene sets, and survival differences were analyzed with Kaplan-Meier survival curves. Mutation data from two IRGPI subgroups were analyzed using R’s Maftools package, and IRGPI scores were correlated with PD-L1 expression and Total Mutation Load (TMB).
To characterize the immune properties of 414 bladder cancer samples, their expression data were imported into CiberSort (https://cibersortx.stanford.edu/) and iterated 1000 times to estimate the relative proportions of 22 immune cell types. We then compared the relative proportions of the 22 immune cells and clinicopathological factors between the two IRGPI subgroups.
To further determine the immune and molecular functions between IRGPI subgroups, we performed ssGSEA on certain gene signatures and compared the scores between the two IRGPI subgroups.
To explore the value of IRGPI in the prognosis of patients after immunotherapy, we also performed time-related ROC curve analysis, obtained the Area under the Curve (AUC), and compared the prognostic value between IRGPI, TIDE and TIS with the time ROC package of R. The TIDE score was calculated online (http://tide.dfci.harvard.edu/login/), and the TIS score was calculated as the mean of log2-scale normalized expression of 18 marker genes.
2.5. Statistical Analysis
Independent t-test was used for comparison of continuous variables between two groups. Categorical data were tested by chi-square test. The Wilcoxon test was used for comparison of TIDE scores between groups. Univariate survival analysis was performed using Kaplan-Meier survival analysis and log-rank test. Multifactorial survival analysis was performed using Cox regression model. Two-sided p values < 0.05 were considered significant.
3. Results
3.1. Immune-Associated Central Genes
In differential expression analysis (414 tumors vs. 19 normal controls), a total of 407 differentially expressed immune-related genes were obtained by crossing these genes with the list of immune-related genes obtained from ImmPort. To obtain immune-related central genes, WGCNA analysis was performed on candidate genes (n = 407). The log log(k) of nodes with connectivity k was negatively correlated with the log log(P(k)) of node probabilities and the correlation coefficient was greater than 0.85. Candidate genes were assigned to 5 modules. The yellow and blue-green modules that were closely associated with bladder tumors were selected for further analysis based on the Pierson correlation coefficient between module and sample characteristics for each module, and genes in these modules were selected for further analysis.
Differentially expressed RBPs were submitted to the STRING database (https://cn.string-db.org/) for PPI detection. PPI networks were then constructed and visualized using Cytoscape 3.7.0. The Molecular Complexity Detection (MCODE) plugin was used to screen key modules from the PPI network with both MCODE scores and node counts > 5. p < 0.05 was considered statistically significant. We obtained the top 50 immune-related hub genes with a threshold of degree > 20. As determined by K-M analysis, expression of 29 immune-related HUB genes was strongly associated with OS in bladder cancer patients (p < 0.05, logarithmic test, Figure 1(a)).
Then, we explored the characteristics of the 29 immune-associated HUB genes. As shown in Figure 1(b), most of the 29 immune-related core genes had missense mutations, among which AHNAK, PDGFRA, SLIT2 and NRP2 had high mutation rates, especially AHNAK mutation rate reached 9%.
Multifactorial Cox regression analysis was performed on 29 immune-related core genes; six independent genes associated with prognosis were identified by multifactorial Cox regression. We developed a prognosis-related risk score model based on the six independent prognosis-related genes. In the model, IRGPI = expression level of AHNAK * 0.49 + expression level of ILK * (−0.24) + expression level of OGN * (−0.18) + expression level of PDGFD * 0.29 + expression level of PPARGC1B * (−0.46) + expression level of JAM3 * 0.30.
Univariate COX regression analysis in the TCGA cohort showed that age, clinical stage and IRGPI risk coefficient were significantly associated with the prognosis of bladder tumors (Figure 1(c)). Further multivariate Cox regression analysis also showed that IRGPI and clinical stage were independent predictors of bladder cancer prognosis (Figure 1(d)). These results suggest that IRGPI and clinical stage are likely to be associated and could be used as bladder tumor a valid biomarker of prognosis. Patient PI values were calculated and patients were divided into high and low risk groups using the median as the cut-off value. Clinical staging of bladder cancer was divided into stage 1 - 4, relative to the two IRGPI subgroups, with stage II predominantly in the low IRGPI subgroup, and stage III, stage IV predominantly in the high IRGPI subgroup (p < 0.001, chi-square test, Figure 1(e)).
PD-L1 expression was higher in the subgroup with high IRGPI than in the subgroup with low IRGPI, with statistical value (p = 0.011), and IRGPI score was correlated with PD-L1 (R = 0.17, p = 0.00073, Figure 2(a)). TMB was slightly higher in the subgroup with low IRGPI, without statistical value (p = 0.074), and with IRGPI score was slightly correlated (R = 0.13, p = 0.0089, Figure 2(b)).
Kaplan-Meier survival analysis was performed. One-way cox regression analysis revealed that the high-risk group had lower OS than the low-risk group (p < 0.001, log test, Figure 3(a)), and then, using the GEO cohort as a validation cohort, the results were consistent with the TCGA cohort (p < 0.001, log test, Figure 3(b)).
3.2. Immunological Characteristics of the IRGPI Subgroup
The CIBERSORT algorithm was applied to calculate the proportion of 22 immune cells, and the Wilconxon test was used to compare the distribution of immune cells in different IRGPI subgroups. Resting memory CD4T cells, M0 macrophages, and M2 macrophages were more abundant in the subgroup with high IRGPI, while CD8T cells, and activated memory CD4T cells were more abundant in the subgroup with low IRGPI (Figures 4(a)-(c)).
![]()
(a)![]()
(b)
Figure 3. (a) K-M survival analysis of the IRGPI subgroup in the TCGA cohort. (b) K-M survival analysis of the IRGPI subgroup in the GEO cohort.
In the K-M survival analysis of 22 immune cells in both IRGPI subgroups, it was found that high infiltration of memory B cells, M0 macrophages, M2 macrophages, resting mast cells, monocytes, and neutrophils in the low IRGPI group had better OS compared with high infiltration of plasma cells, activated memory CD4 T cells, CD8 T cells, and follicular helper T cells in the high IRGPI group (p < 0.05, log test, Figure 4(d)).
![]()
(a)![]()
(b)![]()
(c)![]()
![]()
(d)
Figure 4. (a) IRGPI subgroups and proportion of 22 immune cells in bladder cancer patients in the TCGA cohort. (b) Proportions of 22 immune cell species in different IRGPI subgroups, with thick lines representing median values. The bottom and top of the box are the 25th and 75th percentile (interquartile range), respectively. Differences between the two subgroups were assessed by Wilcoxon test (ns, not significant; *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001). (c) Molecular and immune-related functions of different IRGPI subgroups. Gene sets for molecular and immune-related functions were analyzed by single simple gene set enrichment analysis (ssGSEA) and then compared between different IRGPI subgroups. The scattered points represent the ssGSEA scores of the two subgroups. The thick line represents the median. The bottom and top of the box are the 25th and 75th percentile (interquartile range), respectively. Statistically significant differences between the two subgroups were assessed with the Wilcoxon test (ns: not significant, *p < 0.05, ** p < 0.01, *** p < 0.001). (d) Kaplan-Meier survival curves of specific immune cells in different IRGPI subgroups.
3.3. ICI Treatment Efficiency in the IRGPI Subgroup
The subgroup with low IRGPI had lower TIDE scores than the subgroup with high IRGPI, meaning that patients with low IRGPI were more suitable for ICI therapy (Figure 5). The subgroup with high IRGPI had higher T-cell rejection scores and T-cell dysfunction. There was no difference in microsatellite instability
(MSI) scores between the two subgroups. In the analysis of the TCGA bladder cancer patient cohort, IRGPI had a significantly higher predictive value than TIDE and TIS, with the best at 3 years of follow-up.
4. Conclusions
In this paper, we found six immune-related prognostic genes to establish IRGPI with accuracy and validity in predicting the prognosis of bladder cancer and the efficacy of ICI treatment. Among the six core genes, the mutation rate of AHNAK reached 9%. As a tumor suppressor, AHNAK, also known as desmoglein, negatively regulates triple-negative breast cancer cell proliferation, triple-negative breast cancer xenograft growth, and metastasis through different signaling pathways [13]. In malignant tumors, AHNAK lacks uniform expression and has different intracellular localizations depending on the tumor cell type [14] [15] [16]. AHNAK in neuroblastoma cell lines is described as a nuclear protein whose expression is significantly suppressed [17], whereas, in a study on melanoma, an immune response was observed mainly in the cytoplasm of normal cells, in contrast to the disappearance of expression in melanoma cells [18].
IRGPI is associated with the infiltration of multiple immune cells [19], which contributes to further revealing the role and interrelationship of immune cells in the immune microenvironment of bladder cancer, and as the risk value of IRGPI increases, the later the clinical stage of the patient, the progressive impairment of the main immune function of the body leading to tumor development and metastasis, etc., these differences in the level of infiltration of immune-related cells may indicate the risk value potential immune mechanisms, which could have significant clinical application by detecting relevant infiltrating immune cells in bladder cancer patients to distinguish high and low-risk groups or as immune targets.
The tumor microenvironment consists of the tumor and its surrounding fibroblasts, adipocytes, immune cells, extracellular matrix, and tumor vasculature, which are the environment in which tumors grow, metastasize and evolve [20]. Among them, immune cells, inflammatory cells, and cytokines, which are closely related to immune response, constitute the tumor immune microenvironment and are mainly involved in tumor immune editing [21] [22]. The tumor microenvironment plays an important role in the regulation of PD-L1, and its composition is closely related to the efficacy of ICIs [23].
Tumor-associated macrophages (TAMs), which can be differentiated from CD14-positive monocytes or MDSC, similarly exert immunosuppressive effects and promote tumor metastasis. TAMs are classified into pro-inflammatory type M1 and pro-tumor type M2 according to their polarization status. In contrast, the M2-type TAM, which is subject to chemokine (CCL2), IL-10, and TGF-β-induced polarization, can suppress T-cell immune responses and promote angiogenesis by secreting IL-4, IL-5, IL-6, and ARG1 [24]. M1-macrophages suppress T-cell immune responses and promote angiogenesis by promoting IL-1β, TNF-α IL-6, IL-12, and other cytokine expressions as well as enhancing the function of helper T cell 1 (Th1) to exert anti-tumor effects, and M2-macrophages exert immunosuppressive functions by promoting the expression of immunosuppressive cytokines such as TGF-β1 and IL-10, thereby promoting tumor progression and metastasis [25]. TAMs were found to be abundant in bladder cancer in this study, and M2 macrophages were predominant, which is consistent with previous studies.
In recent years, there has been a lot of exploration of immunotherapeutic markers, and new markers have emerged, such as TIDE and TIS. The Tumor Immune Dysfunction and Exclusion (TIDE) predicted outcomes in melanoma patients treated with first-line anti-PD1 or anti-CTLA4 more accurately than other biomarkers such as PD-L1 levels and mutational load [26]. The results of the RATIONALE-307 phase III clinical trial results confirmed a strong association between TIS score and PD-1 inhibitor combination chemotherapy in delivering significant survival benefits to patients with sq-NSCLC [27]. In the present study, the predictive value of IRGPI was comparable to TIDE and TIS, and at longer follow-up, it proved to be a better predictor of OS, however, as an exploratory study the value of IRGPI in bladder cancer needs to be validated in a multicenter clinical study with a large sample.
In summary, IRGPI established by six immune-related prognostic genes can accurately predict the prognosis of clinical bladder cancer to a certain extent, which will provide new research ideas and directions for the etiology, pathogenesis, individualized immunotherapy, and prognostic judgment of bladder cancer.
Ethics Approval and Consent to Participate
The data used in our study were obtained from public databases TCGA and GEO, therefore, ethical approval was not required.
NOTES
*Corresponding author.