Establishment of a novel risk score model by comprehensively analyzing the immunogen database of bladder cancer to indicate clinical significance and predict prognosis

IFRGs with differential expression

According to the Limma algorithm, we first recognized 4881 differentially expressed BCa genes, of which 1423 were downregulated and 3458 were upregulated (Figure 1A and 1C). To investigate the immune correlations of these genes, we further confirmed 260 IFRGs, including 140 downregulated and 120 upregulated IFRGs (Figure 1B and 1D). We next discovered the functional enrichment of these identified IFRGs by gene ontology (GO). Interestingly, “receptor ligand activity”, “positive regulation of response to external stimulus” and “extracellular matrix” were the most common entries among molecular functions (MF), biological processes (BP) and cellular components (CC), respectively (Figure 2A). Meanwhile, we found, based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, that “cytokine-cytokine receptor interaction” was the most frequently enriched pathway of these IFRGs (Figure 2B). According to the GO and KEGG analysis results, extracellular components are probably the potential terms triggering cell-cell interactions, which inspired us to focus more attention on the elements and requirements of achieving intercellular communication, such as receptors, ligands and extracellular signals, in further studies of BCa.

Figure 1. Differentially expressed BCa genes and IFRGs. Heatmap (A) and volcano plot (C) show the differentially expressed genes between BCa tissues and nontumor tissues. The black dots represent genes without differential expression; the red dots represent the significantly upregulated genes, and the green dots represent downregulated genes. FDR < 0.05, log₂ | FC | >1 and P < 0.05. Immune function-related genes (IFRGs) with various expression levels are illustrated in the heatmap (B) and volcano plot (D). The black dots represent IFRGs without differential expression; the red dots represent upregulated IFRGs, and the green dots represent downregulated IFRGs. FDR < 0.05, log₂ | FC | >1 and P < 0.05.

Figure 2. Functional enrichment analysis of differentially expressed IFRGs. The top pathways of IFRGs are shown in biological process, cellular component, molecular function (A), and KEGG pathway (B).

The characteristics of sIFRGs

With the consideration of predicting prognosis, we explored valid IFRGs as biomarkers to guide clinical treatments. Together with clinical outcomes, we analyzed the survival correlation of these IFRGs, of which 57 were closely involved in the progression of BCa and significantly related to OS (sIFRGs) (P<0.05). To investigate the potency of predicting prognosis and the characteristics of BCa, a forest plot of the 57 sIFRGs was employed to detect the expression profiles in non-BCa and BCa samples (Figure 3 and Table 1). As illustrated in the forest plot, the majority of these genes were upregulated in BCa samples (Figure 3). Consistently, the hazard ratio forest plot revealed that the largest number of sIFRGs played negative roles. The aforementioned results suggest that these upregulated sIFRGs corelating with poor prognosis seem to be the underlying genetic markers.

Figure 3. Survival-related values of sIFRGs. Forest plot of the hazard ratios showing the survival-related values of sIFRGs. Red parts represent upregulated sIFRGs, and green parts represent downregulated sIFRGs.

Through protein-protein interaction (PPI) network analysis, we further found that AGTR1, STAT1, ISG15 and CXCL12 were the core genes among the 57 sIFRGs (Figure 4A). Furthermore, the functional enrichment analysis of the 57 sIFRGs revealed similar results to all 260 IFRGs; for example, “receptor ligand activity”, “regulation of smooth muscle cell proliferation” and “extracellular matrix” were the most enriched terms in MF, BP and CC, respectively (Figure 4B). Unlike the KEGG analysis of the 260 IFRGs, “human cytomegalovirus infection” was identified to be the most enriched among the pathways of sIFRGs, and “cytokine-cytokine receptor interaction” and “neuroactive ligand-receptor interaction” were also enriched in the KEGG analysis of the 57 sIFRGs (Figure 4C), which indicates the significant role of intercellular activity once again.

Figure 4. Functional enrichment analysis of differentially expressed sIFRGs. PPI network (A) of sIFRGs and the top pathways of sIFRGs are shown in biological process, cellular component, molecular function (B), and KEGG pathway (C).

TF-related network

To further detect the underlying mechanisms involved in the clinical differences caused by these sIFRGs, we attempted to discover their regulatory mechanisms. Eventually, 318 transcription factors (TFs) were examined, and we found 77 TFs with differential expression, of which 41 TFs were upregulated and another 37 TFs were downregulated (Figure 5A and 5B). According to the relevant TFs (Supplementary Table 1) and 57 sIFRGs, we then implemented a regulatory network illustrating the relationships among these sIFRGs (Figure 5C). A series of TFs were found to be related to sIFRGs, of which “STAT1” and “NFATC1” showed a closer correlation with their corresponding sIFRGs. The results inspired us to further explore whether these TFs are truly regulated by sIFRGs and whether the variations in TFs play significant roles in the progression and prognosis of BCa in future studies.

Figure 5. The mediated relationships between sIFRGs and differentially expressed transcription factors. Volcano plot (A) and heatmap (B) showing the differentially expressed transcription factors (TFs). The black dots represent TFs without differential expression; the red dots represent upregulated TFs, and the green dots represent downregulated TFs. FDR < 0.05, log₂ | FC | >1 and P < 0.05. The regulatory network (C) was composed of relevant TFs and sIFRGs. The red circles represent upregulated sIFRGs, the green rectangles represent downregulated sIFRGs, and the blue rhombuses represent relevant TFs. The red lines represent positive regulation, and the green lines represent negative regulation.

Clinical consequences assessment

To distinguish the heterogeneous clinical prognostic outcomes, based on the differential expression of sIFRGs, we performed a risk scoring model to divide patients with BCa into a high-risk group and a low-risk group. (Figure 6A–6C). The formula was as follows: [Expression level of CALR * (0.000964)] + [Expression level of MMP9 * (0.000387)] + [Expression level of PAEP* (0.035026)] + [Expression level of RBP7 *(0.014213)] + [Expression level of STAT1 * (-0.00786)] +[Expression level of CACYBP * (0.019418)+[Expression level of AHNAK * (0.016702)]+[Expression level of RAC3 * (0.015394)]+[Expression level of SLIT2* (-0.20438)]+[Expression level of EDNRA * (0.098941)]+[Expression level of IGF1 * (0.264303)]+[Expression level of NAMPT * (0.013366)]+ [Expression level of NTF3 * (-0.66003)]+[Expression level of PPY * (0.014214)]+[Expression level of ADRB2 * (-0.13911)]+[Expression level of SH3BP2 * (-0.0922)]. Based on the divisions of the IRRS, we further detected the survival probabilities of high-risk and low-risk patients with BCa. As illustrated in Figure 7A, low-risk BCa patients obtained higher survival probabilities than high-risk patients, which indicated that the IRRS could be regarded as a significant tool for the potential outcomes of patients with BCa. The area under the ROC curve was 0.792, suggesting moderate potential for the prognostic signature based on IRGs in survival monitoring (Figure 7B). Univariate Cox regression analysis indicated that the risk score, T-stage and N-stage were closely related to OS (p<0.05) (Table 2). In the multivariate Cox regression analysis, only the risk score was ascertained to have potential for predicting OS (Table 3). Therefore, we think it will be reasonable to believe that the “risk score”, which is both significant in univariate and multivariate Cox regression analysis, is promising for use in predicting the prognosis of BCa.

Figure 6. An immunogen-related risk score model (IRRS) was established according to sIFRGs. Heatmap of expression profiles of contained sIFRGs (A). The red parts represent upregulation, the green parts represent downregulation, and the black parts represent no difference. FDR < 0.05, log₂ | FC | >1 and P < 0.05. Distribution of groups and risk score rank (B). Survival status between the low-risk group and the high-risk group (C).

Figure 7. Survival curve and receiver operating characteristic (ROC) curve. The survival curve between the low-risk group and the high-risk group (A). Survival-related receiver operating characteristic (ROC) curve validation of the survival value of the risk score (B).

Clinical use of the risk score signature

To further investigate the relevance of the sIFRGs and clinicopathological features of BCa, we analyzed the relationship between the IRRS and clinical and demographic characteristics, including age, sex, grade, staging, T-stage, N-stage, and M-stage. Under the IRRS, the scores of female patients (Figure 8A), high-grade patients (Figure 8B), patients with advanced N-stage (Figure 8C), advanced staging (Figure 8D) and advanced T-stage (Figure 8E) significantly increased. However, no differences were observed with respect to age or in terms of M-stage. The above results elucidate some clinicopathological and demographic characteristics that are sensitive to the IRRS and further corroborate the clinical application value of the model.

Figure 8. The relationships between the IRRS and different clinical features. Relationships between the IRRS and sex (A), tumor grade (B), N-stage (C), tumor staging (D) and T-stage (E).

We also analyzed the relationships between the compositions of the risk scores and the aforementioned tumor characteristics. As illustrated in Supplementary Figure 1, some sIFRGs, including the higher expression of AHNAK and CACYBP, were closely correlated with advanced grade, staging, T-stage and N-stage. Meanwhile, strong positive relevance was observed between EDNRA expression and age, grade, stage and T-stage. Furthermore, the increased expression of IGF1 obviously positively influenced grade, staging and T-stage. The lower expression of MMP9 showed a negative correlation with grade and T-stage. The expression of NAMPT and SLIT2 was positively related to grade, staging and T-stage. STAT1 is a remarkable IFRG that is differentially expressed in various grades, stages, T-stages and M-stages. The high expression of CALR was correlated with advanced grade and early M-stage. The overexpression of RBP7 has a stronger relationship with more advanced staging and grade. The increased expression of PAEP, PPY and RAC3 was proven to be closely correlated with a higher grade. SH3BP2 overexpression indicated the early N-stage. Detection of the roles of different sIFRGs in indicating tumor characteristics provides insight into the further discovery of biomarkers.

In order to verify whether the immune genome accurately reflects the state of the tumor’s immune microenvironment, we analyzed the relationships between sIFRGs and immune cell infiltration (Figure 9A–9F). We found that macrophages showed the most significant relationship with sIFRGs (Figure 9E), which is probably due to the remarkable effects of macrophage polarization on the tumorigenesis, progression, prognosis and tumor behaviors of BCa. The results motivated us to further discover the underlying functions and mechanisms in future studies.

Figure 9. Relationship between the IRRS and infiltration abundances of six types of immune cells. The relationships were examined using Pearson correlation analysis. B cells (A); CD4 T cells (B); CD8 T cells (C); dendritic cells (D); macrophages (E); neutrophils (F).

CACYBP was overexpressed in patients with BCa, especially with advanced grades and T-stages

To further verify the relationships between sIFRGs and the clinicopathologic features, as well as discover the roles of sIFRGs in indicating clinical prognosis, CACYBP was selected. Consistently, the correlations of CACYBP expression levels and clinicopathological features of BCa as well as OS were similar to those of the integrated sIFRGs (Supplementary Figure 2). We further detected the expression levels of CACYBP in bladder cancer tissues and normal adjacent tissues of 50 patients with various grades, stages and T-stages in vivo. In addition, the expression of CACYBP at the genetic and protein levels was analyzed among normal urothelial cells and various BCa cell lines. As illustrated in Figure 10B, compared with the urothelial cell line SVHUC-1, the higher expression of CACYPB was examined in BCa cell lines including T24, UMUC3, BIU-87 and 5637 at the protein level. Accordingly, consistent with the western blot results, qRT-PCR demonstrated higher expression in BCa cell lines at the genetic level (Figure 10A). Next, we utilized pathological tissue slices to detect CACYPB expression in vivo. We found that the expression of CACYBP was more significant in T3 and 4 stages, grade 3 compared with that in earlier T-stages (Figure 10C) and grades (Figure 10D). The aforementioned results suggest that CACYBP, one of the sIFRGs, has the potential to indicate the clinical prognosis of BCa.

Figure 10. The overexpression of CACYBP in patients with BCa and its correlation with grades and T-stages. The results of RT-qPCR of CACYBP expression in bladder cancer cells and normal urothelial cells (A). Western blot analysis showed that the expression of CACYBP was higher in BCa cell lines than in normal urothelial cells (B). Increasing expression levels of CACYBP were detected by immunohistochemistry (C). CACYBP was also increased in various grades incrementally (D).

Human bladder cell lines

Bladder cancer tissues and normal adjacent tissues were collected from 50 patients admitted to the First Affiliated Hospital of Chongqing Medical University and diagnosed with bladder cancer. Human normal bladder epithelial cell lines SV-HUC-1 and bladder cancer cell lines (BIU-87, TCC-sup, T24, 5637, and UMUC-3) were purchased from the American Type Culture Collection (Manassas, Virginia, USA). Cells were cultured in 1640 and DMEM supplemented with 10% fetal bovine serum (FBS), 100 U/ml penicillin and 100 mg/ml streptomycin (Gibco, Gaithersburg, MD, USA). Cells were incubated at 37°C in 5% CO2. The medium was changed every 1-2 days.

Data download and preprocessing

Transcriptome RNA-sequencing data of BLCA samples were downloaded from the TCGA data portal (https://portal.gdc.cancer.gov/), which contained data from 19 nontumor tissues and 411 primary BCa samples. Clinical data about these patients were downloaded and extracted. Raw data were prepared for further analyses. These data were updated on September 17, 2019. RNA-seq results were combined into a matrix file using a merge script in the Perl language (http://www.perl.org/). Next, the Ensembl database (http://asia.ensembl.org/index.html) was used to convert the Ensembl IDs of genes into a matrix of gene symbols. IFRGs were obtained from the Immunology Database and Analysis Portal (ImmPort) database (https://www.immport.org/). ImmPort is a database that can accurately update immunological data in a timely manner. The data from ImmPort are a significant basis for immunological research. More importantly, a list of IFRGs provided by the database can be used for cancer research. These genes were identified as having an important role in the immune process.

Differential gene analysis

To identify IFRGs participating in the pathogenesis of BCa, the R software limma package (https://bioconductor.org/packages/release/bioc/html/limma.html) was used to screen for differentially expressed genes in cancer and adjacent normal tissues. We present all transcribed differential gene analysis data with the screening value of “FDR < 0.05, log₂ | FC | >1 and P < 0.05”. The differential IFRGs were all extracted from the differential genes. Functional enrichment analysis was performed through the GO and KEGG pathways to explore the underlying molecular mechanisms of differential genes and differential IFRGs.

Survival-related IFRGs

Differentially expressed IFRGs associated with clinical outcomes in patients with BCa were identified as sIFRGs. sIFRGs were selected by univariate Cox analysis using R software survival packages. Functional enrichment analysis was also performed on IFRGs that are closely related to overall survival (OS). Because these sIFRGs may have clinical applications, their clinical value is also worthy of systematic exploration. In order to explore the interaction between these genes, a PPI network based on the data was constructed on the STRING online database (https://string-db.org/). PPI networks can show relationships between many interacting genes. The criterion for a core gene is no less than five node degrees. Cytoscape software version 3.7.2 was used to display PPI results. We further explored their control mechanisms. TFs directly control the degree of gene expression. Therefore, it is necessary to explore the relationship between potentially competent TFs and clinically relevant IFRGs. Cistrome Cancer (http://cistrome.org/) is a data source that integrates cancer genomics data from TCGA with more than 23,000 ChIP-seq and chromatin accessibility analyses to provide regulatory links between TFs and transcriptomes. The Cistrome Cancer database contains a total of 318 TFs for experimental and computational cancer biology research. Then, we constructed an interaction network of prognosis-related TFs and sIFRGs.

Creation of the immunogen-related risk score model (IRRS)

The sIFRGs were submitted for multivariate analysis, while the integrated IFRGs were still used as an independent prognostic indicator to develop the IRRS. The IRRS was based on expression data multiplied by Cox regression coefficients. Patients were divided into low-risk groups and high-risk groups according to the median risk score. The value of IRRS was employed to evaluate various subtypes of patients with BCa. The tumor infiltrating immunocytes were evaluated through the TIMMER database. The level of immune infiltration in patients with BCa was downloaded, and the relationship between IRRS and immune cell infiltration was calculated.

To further investigate the relevance of the sIFRGs and clinicopathological features of BCa, we analyzed the relationship between the IRRS and clinicopathologic characteristics, of which the “TNM staging method” is the most common way to describe the tumor status. The division of “T-stage” BCa was based on the depth and extent of tumor invasion, with a deeper and more extensive invasive status in the more advanced T stages. “N-stage” reflects the lymph node metastasis conditions with more and larger metastatic lymph nodes in the more advanced N stages. “M-stage” is distinguished according to whether the tumor exhibits distant metastasis, and advanced M stages usually represent poor tumor conditions. In addition, “staging” is a comprehensive method combining T-stage, N-stage and M-stage to divide patients with BCa into I, II, III and IV staging.

Western blot

BCa and urothelial cell lines were lysed with RIPA lysis buffer (Beyotime, Haimen, China) and 1% PMSF (Beyotime) at low temperature for 30 min. The protein concentration was quantified using the BCA kit (Beyotime). An equal amount of protein sample was separated by 12% SDS-PAGE and transferred to a polyvinylidene fluoride membrane (Millipore Sigma). The membrane was blocked with 5% skimmed milk at room temperature for 2 hours and incubated at 4°C overnight with appropriate primary antibodies. Then, the membrane was washed with Tween buffer 3 times, and the secondary antibody was incubated for 1 hour at room temperature. Finally, the protein expression level was detected by the chemiluminescence method.

Immunohistochemical staining

Paraffin-embedded BCa tissue sections were dewaxed with xylene, followed by rehydration with reduced concentrations of ethanol (100%, 95%, 85%, and 75%). To obtain the antigen, the sections were immersed in sodium nitrate buffer and heated in sodium citrate for 30 minutes. Then, the endogenous peroxidase activity was removed with 3% hydrogen peroxide and blocked with normal goat serum at room temperature for 15 minutes. The sections were incubated with the primary antibody at 4°C overnight. Following rewarming at room temperature, ascorbic acid was applied at 37°C for 30 minutes. Next, streptavidin-horseradish peroxidase conjugate was added at 37°C for 30 minutes, followed by staining with DAB for 1 min and hematoxylin for 10 s. Finally, PROPLUS (IPP, V.6; MydiyByNeNICs, SilverSpring, MD, (US)) was used to analyze the expression levels of various markers. In addition, the sustained intensity and response intensity were regarded as the parameters to quantitatively measure and evaluate the protein expression levels.

Real-time quantitative PCR

TRIzol (Invitrogen) was used to extract total RNA from BCa and urothelial cell lines under various experimental conditions according to the manufacturer's instructions. The PrimeScript RT kit (Osaka, Japan of TaKaRa) combined with RNA (1 μg) was utilized to reverse transcribe cDNA. The reaction steps were as follows: 37°C for 15 min and 85°C for 5 s. Quantitative polymerase chain reaction (qPCR) was performed on an ABI 7500 real-time PCR system (Applied Biosystems) using the SYBR-Green method (TaKaRa). The reaction cycle conditions were performed (95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 34 s); the primer sequences are shown in Table 4. Three replicate assays were performed for each cDNA sample.

Statistical analysis

Gene function enrichment analysis was based on the R software packages of cluster profiler, org.Hs.eg.db, and enrichplot. The predicted target of the ROC curve is the prognosis of patients with BCa. To verify the prognostic performance, the ROC curve was drawn by the survival ROC package of the R software. The survival time, survival status and risk score of patients with BCa were used to predict the prognosis over a 5-year period. Then, the ROC curve was drawn, and the value of the AUC was calculated. The abscissa is the false positive rate, and the ordinate represents the true positive rate. The AUC generally ranges between 0.5 and 1, and values close to 1 indicate greater accuracy. Variations in clinical parameters were determined using independent t-tests. P<0.05 was considered statistically significant.

Ethics statement

Informed consent forms were signed by all patients before the study. The research protocol has been approved by the Ethics Committee of The First Affiliated Hospital of Chongqing Medical University and is based on the ethical principles of medical research involving human subjects in the Helsinki Declaration.