Identification of the hypoxia-related risk signature to predict bladder cancer prognosis.
Figure 1 showed the flow chart of this study process. We downloaded the RNA-seq and clinical prognostic information in the TCGA-BLCA cohort. The hypoxia-related gene set was obtained from GSEA (hallmark-hypoxia). Next, we developed PPI network analysis by the STRING online database and Cytoscape software to further find out the interactions between these genes (Figure 2A). A total of 50 genes with the most significant interaction degrees were screened (Figure 2B).
To construct the hypoxia-related risk model, univariate and multivariate Cox analysis were analyzed with the top 50 genes in the TCGA-BLCA cohort. In the univariate Cox analysis, 14 HRGs were significantly associated with patients’ OS (Figure 2C). After multivariate Cox regression analysis, 7 HRGs were identified and chosen to construct the signature for OS (Figure 2D). The risk score signature was developed as follows: risk score = (0.119 × SLC2A3 expression level) + (-0.387 × ALDOB expression level) + (0.320 × FOXO3 expression level) + (-0.216 × SDC4 expression level) + (-0.227 × VEGFA expression level) + (0.165 × EGFR expression level) + (0.232 × GPC1 expression level).
Prognostic significance of the hypoxia-related risk signature in bladder cancer patients.
For investigating the prognostic significance of the hypoxia-related risk signature. As shown in the Figure 3A, the expressions of the 7 HRGs were upregulated in high risk score in both TCGA and GSE32894 set. The risk score of patients in the low- and high- risk groups were visualized (Figure 3B). As the hypoxia risk score increased, an increasing rate of mortality in the patients (Figure 3C,D). Moreover, KaplanMeier analysis was used to evaluate the prognostic significance of the hypoxia risk signature. We found that high hypoxia risk score was correlated to poor OS in the TCGA and GSE32894 set compare to low risk score (Figure 3E,F).
The hypoxia-related risk signature for OS is an independent prognostic factor of bladder cancer patients.
To determined the predictive accuracy of the hypoxia risk signature, we used the ROC curve to assess the model. The AUC of the signature for prediction of 1-, 3-, and 5-year OS were 0.661, 0.676 and 0.710, respectively, in the TCGA set and 0.600, 0.594, 0.636 respectively in the GSE32894 set (Figure 4A,B). Then, the univariate and multivariate Cox analysis were used to evaluate whether the independent prognostic value of hypoxia risk signature for OS. The univariate Cox analyses showed that the risk score was associated with OS like other variables including age, gender and WHO grade (Figure 4C). Next, multivariate Cox analysis indicated that the risk score was independently correlated with the OS (Figure 4E). The results were validated in GSE32894 set (Figure 4D,F).
Relationships between the prognostic signature with clinicopathological variables
To investigate whether the prognostic signature correlated with clinicopathological variables in TCGA-BLCA. Bladder cancer patients were stratified according to age, gender, satge, T stage and N stage. The result showed that high-risk patients in those clinical parameters had significantly shorter OS time than low-risk patients (Figure 5), which suggest that the hypoxia-related signature could be applicable to clinical factors.
Correlation between the hypoxia risk signature and immune cell infiltration
For investigating the utility of the risk signature in reflecting the immune cell environment. We used the CIBERSORT analysis to estimate expression level of the 22 immune cell types infiltration between different risk level bladder cancer patients. The composition of the immune cell population in the patients was summarized (Figure 6A). The low hypoxia risk patients performs notably higher proportions of follicular helper T cells (p=0.011), CD8 T cells (p=0.0032) and plasma cells (p=0.0077) compared to the high-risk group (Figure 6E-G). However, a higher proportion of mast resting cells (p=0.023), neutrophils cells (p=0.0061) and CD4 memory resting T cells (p=0.0023) were enriched in high-risk group (Figure 6B-D).
Functional analysis of the prognostic signature
We further verify the underlying mechanism involved in the low- and high-risk groups. GSEA analysis showed that the signaling pathways such as hypoxia, epithelial-mesenchymal transition, inflammatory response and complement were the most significantly enriched pathways in the high risk groups (Figure 7).
Potential of the hypoxia risk signature associated with immunosuppressive microenvironment
Immunotherapy has been a promising treatment for advanced urothelial carcinoma. It was confirmed that Cancer-Immunity Cycle regulating cancer cells and immune response, which affects the utility to immune therapies. The process of Cancer-Immunity Cycle is initiated by the release of cancer-associated antigen. Then the related antigens are identified with dendritic cells and transfer to lymph nodes, accompanied by activating T cells. Those effector cells next migrate and infiltrate the tumor stroma, specifically recognize and eliminate cancer cells. Every step of the cycle needs the coordination of stimulatory and inhibitory factors. Here, we focused on the genes negatively mediating the process in low- and high- risk groups. The related genes signatures were obtained from Tracking Tumor Immunophenotype website (http://biocc.hrbmu.edu.cn/TIP/index.jsp). Genes enriched in the negative regulation of the cycle were significantly increased in the high risk score group (Figure 8A), which indicates that high hypoxia risk patients may associated with poor immunotherapy efficacy.
Moreover, the association between the hypoxia-related signature and the expression levels of important immune checkpoint genes (i.e., PD-L1, PD-1, CTLA-4 and LAG-3) were investigated. As shown in Figure 8B-E, the four immune checkpoints were correlated with hypoxia risk score and upregulated in the high hypoxia risk group. Furthermore, we evaluated the expression of some immunosuppressive cytokines in the low- and high- risk groups. High hypoxia risk group significantly showed a high expression level of immunosuppressive cytokines (Figure 8F).
These results demonstrated that high hypoxia risk patients may develop an immunosuppressive microenvironment and insensitive to immunotherapy.
Validation of 7 HRGs expression results using qRT-PCR analysis
We analyzed the expression of 7 HRGs mRNAs by qRT-PCR in 45 paired paracancerous and cancer tissues. As show in Figure 9 A-G, expression level of SLC2A3, FOXO3, EGFR and GPC1 were higher in tumor samples. No difference was found in ALDOB, SDC4 and VEGFA expression. Moreover, Gene Expression Profling Interactive Analysis (GEPIA) database was used to analyze the HRGs with patients’ OS in TCGA-BLCA. The results showed that high expression of FOXO3 and EGFR and low expression of VEGFA were closely correlated with poorer survival of bladder cancer (Figure 9 H-K). Immunohistochemistry data obtained from the Human Protein Atlas were used to verify the expression of HRGs in normal and tumor tissues was showed in supplemented Figure S1.