DE-FRGs identification and functional enrichment analysis
As listed in TableS1, a total of 293 ferroptosis-related genes were obtained from FerrDb database (http://www.zhounan.org/ferrdb/operations/download.html), GeneCards database (https://www.genecards.org/Search/Keyword?queryString=ferroptosis) and the previous literature. Then we obtained 71 DE-FRGs with the screening conditions mentioned above, including 34 down-regulated FRGs and 37 up-regulated FRGs (Figure 2A). Figure 2B showed the top 10 up- and down-regulated FRGs. Subsequently, we conducted the GO and KEGG analyses using the “clusterProfiler” R package to investigate the potential biological action mechanisms and signaling pathways of these DE-FRGs. As illustrated in Figure 3A, in the BP group, these DE-FRGs were mainly involved in carboxylic, organic, and unsaturated fatty acid biosynthetic process, response to hypoxia, and smooth muscle cell proliferation; in the CC group, these genes were significantly involved in membrane raft, membrane microdomain and apical plasma membrane; in the MF group, the DE-FRGs were enriched in oxidoreductase activity, and iron ion, coenzyme, heme and tetrapyrrole binding. Additionally, KEGG analysis revealed that several carcinogenic signaling pathways, such as Ferroptosis, HIF−1 signaling pathway, central carbon metabolism in cancer and MicroRNAs in cancer, were enriched among the DE-FRGs (Figure 3B)
Development and validation of a prognostic ferroptosis-related gene signature
Using univariate Cox analysis regression, totally 32 DE-FRGs significantly associated with patients’ OS were identified (TableS2). Regulators of the 12 genes were remained using the LASSO regression (Figures 4A-4B), which minimizes the risk of overfitting. Subsequently, we divided all patients randomly into training set (n=374) and test set (n=156) according to the ratio of 7:3. Then, we identify 7 hub FRGs (HMGCR, MT1G, BID, EIF4A1, FOXM1, TFAP2C and CHAC1) and construct a prediction risk signature by stepwise multivariate Cox regression analysis (Figure 4C and Table 2). The HMGCR was considered as protective factors with HR < 1, whereas MT1G, BID, EIF4A1, FOXM1, TFAP2C and CHAC1 were risk factors with HR > 1. The risk score of each patient is calculated as the following formula: Risk score = (-0.430*HMGCR) + (0.077 * MT1G) + (0.471* BID) + (0.318* EIF4A1) + (0.319 * FOXM1) + (0.348* TFAP2C) + (0.364* CHAC1). In addition, we evaluated the prognostic value of in ccRCC, and Kaplan-Meier survival analysis showed that overexpression of MT1G, BID, EIF4A1, FOXM1, TFAP2C, CHAC1 and low expression of HMGCR were associated with the poor OS of patients with ccRCC (Figure 4D-4J). Then, 374 ccRCC patients in the training set were divided into the low-risk group and the high-risk group according to the median value risk score. Kaplan–Meier curve analysis showed that high-risk predicted poor overall survival (Figure 5A). Figure 5B showed the distribution of the risk score, survival status of patients and heatmap plot of the seven hub FRGs between high and low risk group. ROC curve analysis showed that the area under the ROC curve for 1-year, 3-year, and 5-year OS were 0.773, 0.723, 0.731 (Figure 5C). Subsequently, patients in the testing set and the entire group were dichotomized into high-risk and low-risk group based on the median value risk score as the cutoff in the training set. We observed that high risk score indicated a shorted overall survival rate and increased deaths both in the testing set and the entire group (Figure 6A, 6B, 7A and 7B). Additionally, in the testing set, AUC was 0.659 for 1-year OS, 0.689 for 3-year OS, and 0.783 for 5-year OS, respectively (Figure 6C), in the whole group, AUC was 0.736 for 1- year OS, 0.712 for 3-year OS, and 0.749 for 5-year OS, respectively (Fig. 7C).
Identification of independent prognostic factors
In Univariate analysis, we found that 7-FRG based prognostic signature (HR = 1.396, p <0.001), age (HR =1.029; p = 0.007), histological grade (HR = 2.603; p <0.001), clinical stage (HR = 2.091; p <0.001), T status (HR=2.640; p <0.001), N status (HR=4.345; P=0.002) and M status (HR = 4.333; p <0.001) were associated with poorer overall survival of patients with ccRCC (Figure 8A). Furthermore, multivariate survival model after variable selection showed that age (HR = 1.066; p <0.001), M status (HR = 4.548; p =0.005), and the 7-FRG based prognostic signature (HR = 1.284; p <0.001) were independently associated with unfavorable OS of ccRCC patients (Figure 8B).
Relationship between the hub FRGs with clinicopathological parameters
Based on the Wilcoxon-signed rank test, we explored the association between clinicopathological parameters and the mRNA expression levels of seven hub FRGs. We found that patients with lower HMGCR expression level showed a higher histological grade (Figure 9A). As for BID, its gene expression level was correlated with clinical stage, -histological grade, and TMN classifications (Figure 9B-9F). Patients with T3–4 status had higher EIF4A1 expression than patients with T1–2 status (Figure 9G). Increased expression of FOXM1 was significantly related to a high histological grade, advanced clinical stage and high TMN status (Figure 9H-9L). Upregulation of TFAP2C expression was significantly associated with an advanced clinical stage and high T status (Figure 9M-9N). Patients with T3–4 status had higher CHAC1 expression than patients with T1–2 status (Figure 9O). In addition, we observed that a high-risk score associated significantly with an advanced clinical stage, high histological grade and TMN classifications (Figure 9P-9T). However, no significant correlation was found between the MT1G expression level and clinical features of ccRCC.
Functional analyses related to the 7-FRG based prognostic signature
To explore the potential biological mechanisms and signaling pathways that were associated with the 7-FRG based prognostic signature, the DEGs between the high-risk and low-risk groups were used to conduct GO and KEGG analyses. As illustrated in Figure 10A, the main biological processes (BP) involved include extracellular structure organization, epidermis development, extracellular matrix organization, negative regulation of proteolysis, hormone metabolic process, negative regulation of peptidase activity, negative regulation of endopeptidase activity, cellular hormone metabolic process, cornification and acute−phase response. The most abundant cellular component (CC) terminology was collagen−containing extracellular matrix, presynapse, endoplasmic reticulum lumen, blood microparticle, anchored component of membrane, collagen trimer, Golgi lumen, high−density lipoprotein particle, plasma lipoprotein particle and lipoprotein particle. The most abundant molecule function (MF) term was receptor ligand activity, metal ion transmembrane transporter activity, G protein−coupled receptor binding, serine hydrolase activity, serine−type endopeptidase activity, serine−type peptidase activity, cytokine activity, peptidase inhibitor activity, sodium ion transmembrane transporter activity and hormone activity. Results of the KEGG pathway analysis showed that the most abundant pathways were Neuroactive ligand−receptor interaction, Cytokine−cytokine receptor interaction, Protein digestion and absorption, IL−17 signaling pathway, Complement and coagulation cascades, Hypertrophic cardiomyopathy, Arachidonic acid metabolism, Mineral absorption, Linoleic acid metabolism and Phototransduction (Figure 10B).
Association of the 7-FRG based prognostic signature with tumor-infiltrating immune cells and immune checkpoints
To uncover the relationship between the FRG-based prognostic model and tumor immune microenvironment (TME), we explored the difference of tumor-infiltrating immune cells between two groups. The results showed the abundance of Dendritic cells resting, Macrophages M2, Mast cells resting and Monocytes were significantly enriched in the low-risk group (Figure 11A-11D). However, Macrophages M0, Plasma cells, T cells follicular helper and T cells regulatory (Tregs) enriched more in the high-risk group (Figure 11E-11H). Tumor cells could escape from immune response by activating the immune checkpoint pathway. Thus, we estimated the association between the expression level of immune checkpoints and the FRG-based prognostic model, and results showed that PDCD1, CTLA4, LAG3, CD276 and TIGIT were highly expressed in patients with high-risk score (Figure 12A-12E). These results collectively suggest that the changes in the proportion of immune infiltrating cells and the expression of crucial immune checkpoints might be associate with the risk stratification of ccRCC.
External verification of the hub FRGs
Using the ONCOMINE database, we found that HMGCR, MT1G, TFAP2C and CHAC1 were significantly down-regulated in ccRCC samples by comparison with normal renal samples, whereas BID and EIF4A1 were highly expressed in ccRCC specimens compared to noncancerous specimens (Figure 13A and Table 3). However, no difference of FOXM1 relative expression was found between ccRCC tissues and healthy control tissues. We also evaluated the expression levels of the hub FRGs in ccRCC samples and normal samples with GEPIA. Results showed that MT1G and CHAC1 were significantly lower in ccRCC samples than normal renal samples, BID and FOXM1 were significantly higher in ccRCC specimens than noncancerous specimens, however, there was no significant difference in HMGCR, EIF4A1 and TFAP2C expression between ccRCC and normal renal tissues (Figure 13B-13H). In addition, the representative protein expression of HMGCR, BID, EIF4A1, FOXM1, TFAP2C and CHAC1 was investigate in the Human Protein Profiles and shown in Figure 13I-13N. However, MT1G was not found on the website. The value of the hub FRGs in the overall survival of ccRCC patients was also assessed using the Kaplan-Meier plotter database, and results showed that overexpression of MT1G, BID, EIF4A1, FOXM1, TFAP2C, CHAC1 and low expression of HMGCR were related to the poor OS of patients with ccRCC (13O-13U), which are consistent with our above results. Furthermore, we analyzed the genetic alterations of the hub FRGs in ccRCC using the cBioportal database and results were exhibited in Figure 13V.