Single-nucleotide variations, clinical information, and RNA sequencing data of ccRCC patients, along with normal samples, were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). For the training cohort, data from 517 ccRCC patients from TCGA were included after excluding individuals who had a survival time of zero or missing essential information. A validation cohort comprising 446 patients diagnosed with ccRCC was sourced from the cBioPortal database. After excluding individuals with incomplete or missing data, 439 patients were retained for further analysis. Additionally, for further validation, the cohort of GSE53757 dataset composed of 72 individuals utilizing the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array), were enrolled. To assess the translational expression of relevant genes between tumor and normal tissue, data from the Human Protein Atlas (HPA) online database were also utilized (https://www.proteinatlas.org/) (Uhlén et al., 2015).

We collected RNA-sequencing data comprising 537 tumor samples and 72 normal samples from the TCGA database. Utilizing the FerrDb website (http://www.zhounan.org/ferrdb/current/operations/download.html), we identified 844 ferroptosis-related genes. Differential expression and volcano plots of these genes between the tumor and normal tissues were conducted with the “limma” and “ggplot” package. Prognostic ferroptosis-related genes associated with disease progression were identified through univariate and Lasso Cox regression analyses. To normalize the RNA sequencing data, we applied a log2 transformation, after which we constructed gene signatures using the “glmnet” R package. The risk score for this signature was calculated using the following formula: Risk score = ∑ (Expi × βi). The “exp” represents the expression level of each gene and “β” denotes the regression coefficient of each gene.

Incorporating the gene signature score along with main clinical characteristics such as age, gender, stage, and grade, we conducted Cox regression analysis to predict 1-, 3-, and 5-year survival probabilities of this disease. The corresponding score for each factor in model were calculated with R package “nomogramEx.” Subsequently, patients in the training cohort were categorized into two groups based on the median risk score derived from the model.

The OSRC2 and CAKI2 cell lines, commonly employed in renal cancer research (Akifumi et al., 2024; Guodong et al., 2013; Hiroshi et al., 2021; Yasutoshi et al., 2012), were obtained from the American Type Culture Collection (ATCC). These cell lines were cultured in RPMI-1640 medium (Gibco, United States) supplemented with 10% fetal bovine serum (FBS) (Gibco, United States) and 1% penicillin/streptomycin (Gibco, United States). All cell cultures were maintained at 37°C in a 5% CO₂ atmosphere. To knock down specific target genes, small interfering RNA (siRNA) was employed in conjunction with Lipofectamine^® 3000 according to the manufacturer’s protocols (Xiwen et al., 2023). Sequences of siRNAs used in this research is shown in Supplementary Table 1. Corresponding primer sequences for the targeted genes are provided in Supplementary Table 2. The gene silencing efficiency of each targeted gene are presented in Supplementary Figure 1 and the raw data is shown in Supplementary Table 3.

Cells were plated in 96-well plates at a concentration of 2,000 cells per well in 200 μL of growth medium. After incubation for the specified duration, the medium was aspirated, and 20 μL of diluted CCK-8 solution (Dojindo, Japan) was added. The plates were further incubated at 37°C for 1 h, followed by measurement of absorbance at 450 nm (Xiwen et al., 2023). Experiments were performed in three independent groups to ensure consistency and reliability of the results.

All cells were seeded in six-well plates at a density of 1,000 cells per well. The culture medium was refreshed every 3 days. After a ten-day incubation period, the medium was aspirated, and the cells were washed twice with PBS. The colonies were then fixed in methanol for 20 min and subsequently stained with 0.1% crystal violet solution for an additional 30 min at room temperature. Digital images of the plates were captured for permanent recording, and colony counting was performed using ImageJ software (version 1.8.0).

Gene set enrichment analysis (GSEA) version 4.1.0 from the Broad Institute (United States) was employed to identify gene sets that exhibited significant differences between high- and low-risk groups. Enrichment was considered significant if the false discovery rate (FDR q-value) was <0.25, |NSE| was <2, and the normal P value was <0.05. Multiple GSEA plots were generated using R packages “plyr,” “grid,” “ggplot2,” and “gridExtra.”

The cellular component, immune microenvironment and cell immune response data were calculated by different algorithms including TIMER (Li et al., 2016), CIBERSORT (Newman et al., 2015), EPIC (Racle et al., 2017), IPS (Charoentong et al., 2017), MCPcounter (Becht et al., 2016) and quantiseq (Finotello et al., 2019) with the R package “IOBR” (Zeng et al., 2021). A heatmap was created to visualize immune differences between the low and high-risk groups. Additionally, data calculated by ssGSEA algorithm were analyzed by multiple packages “GSVA,” “limma” and “GSEABase.” Potential immune check-points (ICBs) were reported in published literature, and the expression levels of these ICBs between the two groups were illustrated using box plots.

Single-nucleotide variation (SNV) data were obtained from the TCGA database. The genes exhibiting the highest rates of SNV between the two groups were calculated using Microsoft Excel. Waterfall plot was constructed by R package “GenVisR” to display the SNVs associated with our signature.

We employed several statistical tools for our analyses, including R software (version 4.0.2.), GraphPad Prism 5 Software and SPSS 22.0. To assess gene intensity via immunohistochemistry based on different staining types (none = 0, light = 1, moderate = 2, deep = 3), we calculated the scores using the formula: ∑ (percentage* intensity). The R package “rms” was utilized to generate a nomogram to predict 1-, 3-, and 5-year survival probability. R packages “survminer” and “survival” were used to draw Kaplan-Meier plotters to evaluate the discrepancy of overall survival (OS) between different subgroups with chi-squared test. Additionally, a Cox proportional hazards analysis presented hazard ratios (HRs) and 95% confidence intervals (CIs) for various clinical characteristics between groups. Time-dependent receiver operating characteristic (ROC) curves for 1-, 3-, and 5-year survival probabilities were analyzed using the R packages “survivalROC” and “timeROC.” The proportions of immune cells, immune microenvironment, immune checkpoint expression levels and several important genes expressions were compared between groups with Wilcoxon test. A significance level of P < 0.05 was applied, with all p values being two-tailed.

In Supplementary Figure 2, we present the study’s flowchart. A total of 844 ferroptosis-related genes were obtained from FerrDb website. Following differential expression analysis between normal and tumor tissue from TCGA database, 164 genes were selected. Through univariate and lasso cox analysis, eight ferroptosis-related genes were selected (EZH2, AURKA, BID, PLA2G6, EPAS1, SCP2, PRKAA2, ALDH3A2), which significantly influence the prognosis of patients with ccRCC (Figure 1A). The filter variable dynamic process diagram (Figure 1B) and cross verification curve (Figure 1C) were performed with the lambda.1 se criteria to show the process of the LASSO regression. The risk score of the ferroptosis-related gene signature was calculated using the following formula: (−0.000045 × EPAS1) + (−0.002464 × SCP2) + (−0.014891 × ALDH3A2) + (−0.010727 × PRKAA2) + (0.061863 × BID) + (0.026166 × AURKA) + (0.012338× EZH2) + (0.029049× PLA2G6).

Initially, the translated data for these ferroptosis-associated genes was obtained from the Human Protein Atlas database (Figure 2A), in which all genes were identified through immunohistochemistry. The genes EZH2 was significantly higher expressed in tumor tissue while SCP2, ALDH3A2, and PRKAA2 were significantly upregulated in normal tissue. Figure 2B illustrates the transcriptional profiles of these eight genes in the TCGA dataset, revealing expression patterns that were consistent with those observed in Figure 2A. In ccRCC tissue, compared to normal kidney tissue, EPAS1, BID, AURKA, EZH2, and PLA2G6 exhibited significantly elevated expression levels. Conversely, the expression of SCP2, ALDH3A2, and PRKAA2 significantly reduced in ccRCC tissue. Supplementary Figure 3 illustrates the transcriptional data for various predictive ferroptosis-related genes from the GSE53757 dataset. The expression results were generally in agreement with those in the TCGA dataset except PLA2G6 gene, which exhibited higher expression in normal kidney tissue.

We identified the biological functions of eight ferroptosis-related genes in clear cell renal cell carcinoma (ccRCC). The results from the CCK-8 assay demonstrated that silencing EPAS1 and ALDH3A2 resulted in increased cell viability in OSRC2 and CAKI2 cell lines, while silencing BID, AURKA, and PLA2G6 led to a decrease in cell viability (Figure 3A). Additionally, wound healing assays revealed that the knockdown of EPAS1 and ALDH3A2 enhanced the invasive capabilities of OSRC2 and CAKI2 cells, whereas silencing BID, AURKA, and PLA2G6 significantly reduced their invasiveness (Figure 3B). Plate Colony formation assays further indicated that the colony formation ability of OSRC2 and CAKI2 cells increased following the silencing of EPAS1 and ALDH3A2, while a decrease in colony formation was observed after silencing BID, AURKA, and PLA2G6 (Figure 3C; Supplementary Figure 4). These results suggest that ferroptosis-related genes play a significant role in regulating renal cell carcinoma metastasis.

The risk score plots derived from the ferroptosis-related gene signature in ccRCC patients and normal individuals were analyzed and presented in Figures 4A, B. These analysis reveals a significant elevation in risk scores for ccRCC when compared to normal kidney tissue, highlighting the potential role of ferroptosis-related genes in the pathogenesis of ccRCC.

The ferroptosis-related gene signature exhibited potential clinical value in diagnosing ccRCC with AUC = 0.937 (95%CI: 0.910–0.964), as illustrated in Figure 4C, which not only highlights the effectiveness of this gene signature in distinguishing ccRCC patients from healthy individuals but also suggests its potential as a valuable biomarker for early diagnosis and treatment planning.

A ferroptosis-related gene model and nomogram were constructed by multivariate cox regression analysis (Figure 5A), which were composed with ferroptosis-related gene signature score, age, sex, stage and grade of tumor. Age points were calculated with formula: age points = 0.68 * age −16.92. In this model, the point for male patients was set at 0.56, while female patients received 0 points. The points for stages “1, 2, 3, and 4” are as follows: “0.26, 0, 14.85, and 33.29,” respectively, while the points for grades “1, 2, 3, and 4” are “0, 92.07, 96.98, and 100”, respectively. The points for the ferroptosis-related gene signature were computed using the formula: points = 37.94 * signature score + 15.18. These cumulative points were utilized to predict the survival probabilities of patients at 1, 3, and 5 years. Based on the total score, patients were stratified into two groups, with a median total score of 151.84 points serving as the cutoff. This classification was applied consistently in both the training and validation cohorts.

The Kaplan–Meier analysis revealed that high-risk ccRCC patients which was identified by ferroptosis-related gene model, were more likely to experience poor overall survival compared to low-risk individuals in the training cohort (TCGA dataset) (HR 4.243, 95% CI = 3.119–4.243, p < 0.0001) (Figure 5B). Similar findings were observed in the external testing cohort (cBioPortal database) (Figure 5C). The survival status plots, and risk score plots of ferroptosis-related gene model in two groups were displayed in Figures 5D, E, respectively, indicating more death cases as the risk score increase in ccRCC patients, showing a good distinguished power of the model. The area under curve (AUC) values for predicting survival probability at 1, 3, and 5 years were 0.875, 0.818, and 0.790, respectively, in the training cohort (Figure 5F). Notably, the ferroptosis-related gene signature exhibited greater accuracy in predicting the 5-year survival rate of ccRCC patients compared to individual ferroptosis-related genes (Figure 5G). Furthermore, it exhibited superior predictive efficacy for 5-year survival when compared to conventional clinical characteristics such as age, gender, grade, stage and TMN classification (Figure 5H).

Subgroup analyses of the ferroptosis-related gene signature were conducted among patients diagnosed with ccRCC to evaluate the model’s performance across various clinical characteristics (Figure 6). These analyses revealed that high-risk patients consistently exhibited poorer overall survival rates compared to low-risk patients across multiple subgroups. This trend was evident in patients categorized by tumor stages T1 to T4 (Figures 6A–C), as well as those classified with N0 (no regional lymph node involvement) and M0 (no distant metastasis) stages (Figures 6D–H). Notably, the model’s predictive accuracy extended to all stages of disease, encompassing stages 1–4 (Figures 6I–K). Moreover, both male and female patients classified as high-risk demonstrated significantly worse survival outcomes compared to low-risk individuals (Figures 6L, M). Additionally, the analysis included different age groups, and the results indicated that high-risk patient experienced diminished overall survival compared to their low-risk peers (Figures 6N, O).

Gene set enrichment analysis (GSEA) was employed to investigate the differences in biological profiles between high-risk and low-risk patients with ccRCC. Pathways such as the intestinal immune network for mTOR signaling, fatty acid metabolism, the citrate cycle (TCA cycle), renal cell carcinoma, and peroxisome pathways were enriched in high-risk patients, suggesting the alterations in this pathway may contribute to the aggressive behavior of tumors in high-risk patients. Conversely, the low-risk group displayed enrichment in pathways related to immune response and tissue rejection, including the intestinal immune network for IgA production, allograft rejection, graft versus host disease, and proteasome pathways (Figure 7), which suggest that the low-risk patients may have a better immune response, potentially facilitating better control of tumor growth and progression.

A heatmap was generated using CIBERSORT, EPIC, IPS, MCPcounter, quantiseq and TIMER algorithms to illustrate the immune landscape in high- and low-risk groups (Supplementary Figure 5). Differences in infiltration of 16 types of immune cells between two groups were analyzed using ssGSEA (Figure 8A). Cancer-inhibitory immune cells, including activated dendritic cells (aDCs), CD8⁺ T cells, macrophages, tumor cell infiltrating lymphocytes (TILs) and T regulatory cells (Tregs), were significantly higher in the high-risk group, suggesting a stronger anti-tumor immune response in these patients. Additionally, thirteen immune functional scores, such as CC chemokines (CCR), checkpoints, cytolytic activity, human leukocyte antigen (HLA), and inflammation-promoting factors, showed significant differences between these two groups. Specifically, the expression levels of immune checkpoints such as PD1, CTLA4, LAG3, TIGIT, LILRB, and BTLA were significantly elevated in the high-risk group (Figure 8B). These findings suggest a potential therapeutic advantage from immunotherapy for high-risk patients, highlighting the role of these pathways in cancer progression.

We analyzed all genes exhibiting single nucleotide variant (SNV) mutations, ultimately identifying the 15 genes with the highest mutation rates, as illustrated in Figure 9A. Among these, VHL displayed the highest mutation rate (60.6%), followed by PBRM1 (50.2%), TTN (20.8%), SETD2 (15.2%), BAP1 (12.3%), etc. The prevalence of mutations in these genes showed their potential significance in the pathogenesis of ccRCC. Interestingly, our analysis revealed that patients with mutations in PBRM1 had significant lower risk scores compared to those with a wild-type background. In contrast, mutations in SETD2 and BAP1 were correlated with increased risk scores when compared to their wild-type counterparts (Figure 9B), indicating a potential association with more aggressive disease characteristics and poorer outcomes. Kaplan-Meier analysis revealed that patients with a BAP1 mutation exhibited poorer overall survival (OS) compared to those with wild-type BAP1 (Figure 9C), suggesting that BAP1 mutation may serve as an independent prognostic factor for predicting the outcome of ccRCC patients.