Identification of a Novel Ferroptosis-Related Gene Signature for Prediction of Prognosis in Bladder Urothelial Carcinoma

Data acquisition

RNA sequencing data from 411 BLCA tissues were extracted from the Cancer Genome Atlas (TCGA) and corresponding clinicopathological data were obtained as a training set. The RNA sequencing transcriptome data were estimated as log2(X + 1) transformed FPKM normalized counts. Similarly, the microarray data from 165 primary urothelial tissue samples were extracted from the GSE13507 dataset published in the Gene Expressed Omnibus (GEO) and used as a validation set. Data on 60 ferroptosis-related genes were retained mainly from the findings of previous literature [7, 16–19]. The clinical information of both TCGA and GEO cohort was presented in Table 1.

Visualization of interaction among FRGs

The protein-protein interaction network (PPI) was constructed using the STRING database (https://string-db.org) and the functional network was determined with a medium confidence score greater than 0.4, along with other default parameters.

Establishment of a ferroptosis related-gene risk model

Samples without complete clinical data were removed, and Univariate Cox analysis was performed on the remaining samples to identify potential prognostic factors for patients with BLCA. We screened for ferroptosis-related genes associated with prognosis based on P < 0.05. LASSO regression was then performed to eliminate FRGs that might overfit our ferroptosis risk model. The LASSO algorithm was used for variable selection and shrinkage with the ‘glmnet’ R package. Additionally, we adopted a multivariate Cox analysis to determine favorable prognostic FRGs according to the linear combination of the Cox coefficient and gene expression. The following calculation formula was used for the analysis: Risk score = ∑i=1N(Exp×Coef) . N, Coef, and Exp represent gene number, coefficient value, and level of gene expression, respectively. All bladder cancer patients were divided into high-risk and low-risk groups by using the median as the boundary value. We generated a receiver operating characteristic (ROC) curve to verify the accuracy of the risk model for OS prediction using the ‘survivalROC’ R package. ‘Survival’ and ‘SurvMiner’ packages were used to analyze the overall survival status of BLCA patients, and Kaplan-Meier (KM) analysis was performed to compare the differences between high-risk and low-risk ferroptosis groups.

Characterization of genetic alteration and DNA methylation among key FRGs

The cBioPortal platform (http://www.cbioportal.org/) was adopted to gather information on genetic alteration and DNA methylation of hub FRGs in BLCA patients. Based on the cBioPortal platform, the OncoPrinter module was used to analyze the frequency of gene alteration. Thereafter, the Plot module was used to calculate and visualize the correlation between relative linear copy number values (CNV), DNA methylation, and mRNA expression.

Gene set enrichment analysis (GSEA) for high ferroptosis risk subgroup

GSEA was performed to detect whether there was a significant difference in gene expression between high-risk and low-risk subgroups in the enrichment of ‘MSigDB’ collection [h.all.v7.2.symbols.gmt (hallmarks)]. Gene set permutations were performed 1,000 times for this analysis. The phenotype label was represented by a high-risk or low-risk score.

Prediction of chemotherapy response between high and low risk groups

The chemotherapy response for cisplatin, Gemcitabine, docetaxel and sunitinib of each bladder cancer patient in the TCGA cohort was calculated based on the Genomics of Drug sensitivity in Cancer (GDSC, https://www.cancerrxgene.org) by using “pRRophetic” Rpackage.

Statistics

The R language platform was used to analysis these data (v.4.0.3, https://cran.r-project.org/). The wilcoxon rank-sum test was performed to distinguish the difference of mRNA levels of FRGs between different risk groups. The univariate cox and lasso regression were used to screen hub FRGs and the multivariate cox regression was used to establish our model. The KM curve was generated using R package ‘survival’ and assessed by the log-rank test. The corresponding coefficients of hub FRGs weighted and obtained from multivariate cox regression were used to calculate the risk score of each BLCA patient.

Identification of ferroptosis signature for prediction of BLCA prognosis

In this study, 60 ferroptosis-related hub genes were adopted based mainly on the previous study of Liang JY, et al. [19]. We performed a protein–protein interaction network analysis using the STRING database to better understand the interaction between the 60 ferroptosis-related genes (Fig. 1A). Generally, proteins that are more instrumental exhibit a considerable number of adjacent nodes. Thus, we used a ‘barplot’ R package to identify the number of adjacent nodes in each protein, and TP53 presented the maximum (Fig. 1B).

Fig. 1

The relationship between 60 ferroptosis-related genes and Establishment of a ferroptosis-related gene risk model. (A) The protein-protein interaction network of 60 ferroptosis-related genes. (B) The number of adjacent nodes that each ferroptosis-related gene possessed. (C) Relationship between the expression of FRGs and overall survival of bladder urothelial carcinoma patients by univariate Cox analysis. (D) The Lasso coefficient profiles of FRGs in BLCA. (E) Partial likelihood deviance was plotted against log (lambda). The vertical dotted lines represent the lambda value with minimum error. The largest lambda value is where the deviation is within one standard error (SE) of the minimum. (F) Establishment of a ferroptosis risk signature to predict BLCA prognosis by multivariate Cox regression. (G) Spearman correlation analysis of six ferroptosis-related genes in the TCGA cohort. The red circles represent positive correlation,the green circles represent negative correlation. The area of the circles represent the value of correlation. (H) Spearman correlation analysis of six ferroptosis-related genes in the GEO cohort.

To establish a ferroptosis risk signature for prediction of prognoses in patients with BLCA, univariate and multivariate regression analyses were performed using the 60 FRGs in TCGA training cohort. Accord-ing to univariate COX analysis (Fig. 1C), 12 ferro-ptosis-related genes were significantly correlated with overall survival (OS). We performed LASSO regression to eliminate FRGs that might overfit our model (Fig. 1D, E). Finally, nine genes were left to be subordinated into the multivariate Cox model. As visualized by multivariate COX analysis, six ferroptosis-related genes were selected to construct the predictive model (Fig. 1F). The risk-score formula was generated as follows:

Risk score =

(0.14×FADS2) + (0.18×GCLM) +

(-0.11×ALOX5) + (-0.31×ACSL4) +

(0.33×ACACA) + (0.19×CRYAB)

All six FRGs were found to be not significantly correlated with each other in both TCGA and GEO cohorts (Fig. 1G, H), indicting our risk model based on these six FRGs avoided the overfitting caused by collinearity.

Identification of genetic alteration status and methylation among six FRGs

Genetic alteration status was explored using the cBioportal platform. The results demonstrated that CRYAB, ACLS4, GCLM, FADS2, ALOX5, and ACACA exhibited 1.0%, 1.5%, 1.5%, 1.5%, 3.0%, and 9.0%genetic alteration, respectively (Fig. 2A), and most of the alterations pertained to variations in the copy number. The correlation between relative copy number values and mRNA expression levels of ACACA reached a value of 0.42 (Fig. 2G), as determined by Pearson’s correlation analysis, while that of other genes was less than 0.30 (Fig. 2B-F). We further investigated the association between mRNA expression and methylation status of six FRGs, although methylation data for only four of the six FRGs was available to perform the comparison. The results showed that three FRGs, namely GCLM, FADS2, and ALOX5, had robust negative correlations be-tween gene expression level and DNA methylation (Fig. 3A-D).

Fig. 2

The genetic changes of six hub FRGs. (A) The genetic alterations of six hub genes. (B-G) The correlation between mRNA expression and relative linear copy number values.

Fig. 3

The methylation of hub FRGs.(A-D) The correlation between mRNA expression and DNA methylation of hub FRGs.

Identification of prognostic value of ferroptosis-risk model in BLCA

In consideration of the significant biological functions of ferroptosis to facilitate tumor growth and progression, we conducted a systematic investigation of the prognostic value of the ferroptosis model. The patients were divided into a high ferroptosis-risk subgroup and a low ferroptosis-risk subgroup according to the median value in both TCGA (Fig. 4A) and GEO (Fig. 4E) databases. As the expression levels of CRYAB, FADS2, ACACA, and GCLM were enhanced with higher ferroptosis risk scores, while the expression level of the two remaining FRGs decreased with higher risk scores in both TCGA (Fig. 4B) and GEO (Fig. 4F) databases, it was proposed that ALOX5 and ACSL4 served as protective factors. Our data also indicated that patients with higher risk scores had a higher mortality risk than those with lower scores in TCGA database (Fig. 4C). The results were validated in the GEO database (Fig. 4G). Additionally, we performed KM analysis to assess the prognostic potential of the ferroptosis risk model in patients with BLCA. As shown in Fig. 4D, a high ferroptosis risk score was associated with poor overall survival in TCGA -BLCA training set, which was authenticated by the GEO testing set (Fig. 4H).

Fig. 4

Risk score analysis of ferroptosis-gene risk model in TCGA and GEO cohort.(A) separation of patients into high and low risk group by the median risk score.(B)Heatmap plot of six genes between high and low risk group.(C)Survival status of patients.(D)Survival analysis according to risk score. (E) separation of patients into high and low risk group by the median risk score.(F)Heatmap plot of six genes between high and low risk group.(G)Survival status of patients.(H)Survival analysis according to risk score.

FRG risk model shows the capability to predict prognosis

To estimate the predictive accuracy of the ferroptosis risk model in the one-to-five year survival rate group, a receiver operating (ROC) curve was performed using the information from TCGA and GEO cohorts. It is shown in Fig. 5A that the proportions under the ROC curve (AUC) were 0.630, 0.671, and 0.703, at one year, three years, and five years, respectively. This figure demonstrates the predictive value of our model. In parallel, we performed another ROC curve to assess the prognostic prediction manifestations of our ferroptosis risk signature in GEO datasets (Fig. 5B). Fortunately, we obtained a much larger AUC score, indicating the accuracy of our model to predict the overall survival rate of BLCA patients.

Fig. 5

Prognostic value of the ferroptosis risk signature in BLCA.(A,B)ROC curves indicating the predictive efficiency of the ferroptosis risk signature on the 1-, 3-, and 5-years survival rate;(C-F)Univariate and multivariate Cox analyses evaluating the independent prognostic value of the ferroptosis signature in terms of OS in BLCA patients.

Additionally, we performed univariate and multivariate Cox regression analyses to assess the independent prognostic value of the ferroptosis risk model in terms of overall survival. The univariate analysis indicated that higher ferroptosis risk scores significantly contributed to an unfavorable OS (Fig. 5C) and was further validated in GEO cohort (Fig. 5D). Other parameters were also embedded into this analysis, including age, sex, TNM stage, and papillary levels. Multivariate analysis demonstrated that high ferroptosis risk scores were independently correlated with unfavorable OS in these patients (Fig. 5E), indicating that ferroptosis risk score is a potential independent prognostic factor for patients with BLCA. This was further validated by the GEO datasets with the variables of grade and invasiveness (Fig. 5F).

Ferroptosis-related gene is correlated with clinicopathological characteristics and molecular subclassifications

In consideration of the robust biological functions of ferroptosis in tumor growth and development, we conducted a comprehensive investigation to evaluate the correction among six identified ferroptosis-rela-ted genes and the clinicopathological characteristics as well as the molecular subclassifications of bladder carcinomas, including the T-stage, N-stage, WHO grade along with PAM, TCGA, Lund and CMS molecular subclassifications. We classified the T-stages into T1, T2 and T3&4, and then generated a heatmap of gene expression levels and T-stages in TCGA datasets (Fig. 6A). The map showed that the expression level of ALOX5 was significantly decreased in the T3&4 stage tumor groups compared to the other three genes with statistical significance, implying that ALOX5 might be a protective factor in the low ferroptosis risk group. Similarly, Fig. 6C visualized the heatmap of gene expression levels and N-stages in TCGA cohort. Furthermore, we studied the relationship between gene expression and the WHO grade. The heatmap demonstrated that ALOX5 expression level was significantly decreased in high-grade BLCA (Fig. 6E). The correlations between T-stages, N-stages, WHO grade and gene expression levels were authenticated by quantitative analyses (Fig. 6B, D, F). With the development of second-generation sequencing, molecular subclassification of bladder cancer were established which contributed greatly to the individual treatment [20–22]. Hence, we explored the relationships between Six FRGs and molecular subclassifications. Various molecular subclassifications data for each bladder cancer patient were obtained from the previous research [23]. Figure 6G-J showed the heatmap of the expression of six FRGs in different molecular subclassifications. Results demonstrated that the expression levels CRYAB, FADS2 and GCLM were overexpressed in basal group. In contrast, patients with high expression of ALOX5 seems to generate a luminal subclassification.

Fig. 6

Ferroptosis-realted gene expression is correlated with clinicopathological features and molecular subclassifications of BLCA. (A)Heatmap showing six ferroptosis gene expression profiles in different T-stage from TCGA cohort. (B) The expression level of six FRGs in BLCA with different T-stage. (C) Heatmap showing six ferroptosis gene expression profiles in different N-stage from TCGA cohort. (D)The expression level of six FRGs in BLCA with different N-stage. (E) Heatmap showing six ferroptosis gene expression profiles in different WHO grade from GEO cohort. (F) The expression level of six FRGs in BLCA with different WHO grade. (G-J) Heatmap showing six ferroptosis gene expression profiles in different molecular subclassifications in TCGA cohort. ^*P < 0.05,^**P < 0.01,and ^***P < 0.001

GSEA demonstrates ferroptosis-related signaling pathways

GSEA was performed to compare the high-risk ferroptosis and low-risk ferroptosis groups, to verify related signaling pathways activated in the high-risk subgroup. Gene sets were differentially enriched in the high ferroptosis groups of TCGA database, as they were related to tumor proliferation such as G2M checkpoint, E2F targets, and hedgehog signaling (Fig. 7A). These were further validated in the GEO database (Fig. 7B).

Fig. 7

GSEA enrichment between high and low ferroptosis risk groups. (A)GSEA revealing that genes in the high ferroptosis risk group were enriched for hallmarks of malignant tumors in TCGA-BLCA cohort. (B)The results were further validated by GEO data.

The relationship between ferroptosis-related signature and molecular subclassification, histological type and chemotherapy response

Previous work demonstrated that the expression of six FRGs were strongly correlated with various molecular subclassification. Hence, we investigated the relationship between our risk signature and various molecular subclassifications. Results demonstrated that patients in basal group tend to possess higher risk score indicating unfavorable prognosis (Fig. 8A). The relationship between our risk signature and TCGA, Lund and CMS molecular subclassifications were showed in Fig. 8B-D. Additionally, we also investigated the association between angiolymphatic invasion and ferroptosis-related gene signature. Results showed that patients with angiolymphatic invasion tend to own higher risk scores (Fig. 8E). Moreover, we subsequently detected the correlation between our risk signature and different histological types. Figure 8F revealed that patients with MIBC were correlated with higher risk scores, indicating that our risk signature had the ability to distinguish MIBC patients. Chemotherapy is the common method for treating patients diagnosed with bladder cancer. Therefore, the difference responses for cisplatin, gemcitabine, docetaxel and sunitinib between high and low risk groups were calculated based on the GDSC database. The estimated IC50 values demonstrated that patients in high-risk group had a better response for cisplatin, docetaxel and sunitinib (Fig. 8G).

Fig. 8

Association between ferroptosis-related signature and molecular subclassifications and chemotherapy response. (A-D) the association between ferroptosis-related signature and PAM, TCGA, Lund and CMS subclassification. (E) the association between angiolymphatic invasion and ferroptosis-related signature. (F) the relationship between ferroptosis-related signature and histological type (NMIBC and MIBC). (G) Various drug response between high and low risk groups.