Abstract

Background. Establishing a robust signature for prognostic prediction and precision treatment is necessary due to the heterogeneous prognosis and treatment response of clear cell renal cell carcinoma (ccRCC).

Objectives. This study set out to elucidate the biological functions and prognostic role of ferroptosis-related long non-coding RNAs (lncRNAs) based on a synthetic analysis of competing endogenous RNA networks in ccRCC.

Materials and methods. Ferroptosis-related genes were obtained from the FerrDb database. The expression data and matched clinical information of lncRNAs, miRNAs and mRNAs from The Cancer Genome Atlas (TCGA) database were obtained to identify differentially expressed RNAs. The lncRNA-miRNA-mRNA ceRNA network was established utilizing the common miRNAs that were predicted in the RNAHybrid, StarBase and TargetScan databases. Then, using progressive univariate Cox regression, least absolute shrinkage and selection operator (LASSO) and multivariate Cox regression analysis of gene expression data and clinical information, a ferroptosis-related lncRNA prognosis signature was constructed based on the lncRNAs in ceRNA. Finally, the influence of independent lncRNAs on ccRCC was explored.

Results. A total of 35 ferroptosis-related mRNAs, 356 lncRNAs and 132 miRNAs were sorted out after differential expression analysis in the TCGA-KIRC. Subsequently, overlapping lncRNA-miRNA and miRNA-mRNA interactions among the RNAHybrid, StarBase and TargetScan databases were constructed and identified; then a ceRNA network with 77 axes related to ferroptosis was established utilizing mutual miRNAs in 2 interaction networks as nodes. Next, a 6-ferroptosis-lncRNA signature including PVT1, CYTOR, MIAT, SNHG17, LINC00265, and LINC00894 was identified in the training set. Kaplan–Meier analysis, PCA, t-SNE analysis, risk score curve, and receiver operating characteristic (ROC) curve were performed to confirm the validity of the signature in the training set and verified in the validation set. Finally, single-sample gene set enrichment analysis (ssGSEA) and ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) analysis showed that the signature was related to immune cell infiltration.

Conclusions. Our research underlines the role of the 6-ferroptosis-lncRNA signature as a predictor of prognosis and a therapeutic alternative for ccRCC.

Key words: bioinformatics analysis, ferroptosis-related lncRNA, competing endogenous RNA network, prognostic signature, clear cell renal cell carcinoma

Background

Based on the data from the American Cancer Society, an estimated 79,000 new cases of kidney cancer occurred in the USA in 2022, with nearly 13,920 deaths attributed to this disease.¹ Histologically, renal cell carcinoma (RCC) constitutes the vast majority (90%) of kidney cancer cases, primarily comprising clear cell RCC (ccRCC; 70%), papillary RCC (10–15%) and chromophobe RCC (5%).² The incidence of RCC has been gradually increasing due to population aging, obesity and environmental pollution.¹ According to the latest global data, about 4.4 per 100,000 people are diagnosed with RCC every year. The incidence of RCC is twice as common in men as in women, which poses a significant health risk to men.^{3, 4}

Early surgical intervention is currently the main treatment method for clear cell RCC (ccRCC), with an overall 5-year survival rate of >90% for most patients who undergo early resection surgery.^{3, 4} However, about 30% of patients experience metastatic recurrence after nephrectomy, which significantly affects postoperative survival.⁵ Additionally, ccRCC demonstrates high resistance to chemotherapy and radiotherapy, leading to a poor prognosis.⁶ According to reports, the 5-year survival rate for patients with advanced ccRCC is only 11.7%.⁷ The survival outcomes of ccRCC patients vary greatly, making it difficult for clinicians to accurately predict patient prognosis. Traditionally, the clinical and pathological features of patients have been used to evaluate the risk of recurrence and predict prognosis. Recently, researchers have explored molecular biomarkers to reliably predict the prognosis of ccRCC.⁸ In 2018, a validated prognostic molecular signature, ClearCode34, was introduced. Based on the expression of 34 genes, ClearCode34 demonstrated satisfactory prediction performance.⁹ Moreover, a study published in 2022 by Helmink et al. utilized deep RNA sequencing to construct B-cell-related gene signature, aiming to predict the impact of anti-PD-1 therapy on the efficacy and prognosis of ccRCC patients.¹⁰ It can be seen that by integrating RNA sequencing data with survival information, molecular-based characteristics provide new options for assessing the prognosis of ccRCC patients. Therefore, considering the high incidence and mortality rate of ccRCC, it is crucial to explore molecular features with prognostic value that affect ccRCC patients.

Cell death induction, including apoptosis, necroptosis, pyroptosis, and autophagy, is the core mechanism of anti-tumor drugs. Approximately 10 years ago, the concept of ferroptosis was introduced as a new form of programmed cell death.¹¹ Ferroptosis is an iron-dependent form of regulatory cell death caused by excessive lipid peroxidation and is involved in the development and progression of different tumors, including ccRCC.¹² Currently, inducing ferroptosis in ccRCC is a promising strategy. Miess et al. has proved that inhibiting the synthesis of glutathione can make ccRCC sensitive to ferroptosis and finally prevent tumor growth.¹³ In addition, a study by Yang et al. also confirmed that the Hippo pathway effector TAZ can regulate the ferroptosis sensitivity of RCC.¹⁴ Therefore, modulating ferroptosis may have a number of important implications for future RCC therapeutic practices.

Long non-coding RNAs (lncRNAs) are non-coding transcripts containing more than 200 nucleotides.¹⁵ Accumulating studies have shown that lncRNAs play a significant role in the occurrence, development and metastasis of ccRCC, and can serve as reliable prognostic factors.^{16, 17} Instead of using a single lncRNA to analyze its prediction of disease, it is more effective to comprehensively analyze the expression profile of certain pathway-related lncRNAs. However, only a few studies have screened ferroptosis-related lncRNAs in ccRCC. Moreover, very little information is available on lncRNA signatures to explain the relationship between lncRNAs and ferroptosis-related genes using the competing endogenous RNA (ceRNA) network.

Objectives

Our objective was to establish a ferroptosis-related lncRNA signature for predicting the prognosis of ccRCC patients. Furthermore, by constructing ceRNA networks and analyzing the immune microenvironment, we sought to elucidate the underlying mechanisms of ferroptosis-related lncRNAs in ccRCC.

Materials and methods

Data acquisition and differentially expressed gene analysis

Overall, 112 ferroptosis-related genes (Supplementary Table 1) were obtained from the FerrDb Database (http://www.zhounan.org/ferrdb) concerning the following screening conditions: validated in human and protein-coding. Gene expression data and matched clinical profiles (including lncRNAs, microRNAs (miRNAs) and messenger RNAs (mRNAs)) of the KIRC cohort were downloaded from The Cancer Genome Atlas (TCGA; https://www.cancer.gov/ccg/research/genome-sequencing/tcga) database using the “GDCRNATools” package of R software,¹⁸ duplicate samples were removed, and only samples with sample_type of “PrimaryTumor” and “SolidTissueNormal” were retained. Then, the “GDCRNATools” package of R software was applied to analyze the differential expression of lncRNA, miRNA and ferroptosis-related mRNA based on the conditions: method: DESeq2, Normalization: Voom, p = 0.05 and log| Fold Change |=2.

Construction of the ceRNA network

The miRNA can regulate gene expression by targeting the 3’UTR of mRNA. The ceRNA network refers to non-coding RNA, such as lncRNA, that can competitively bind to miRNAs and reduce their inhibition on mRNA. The interactions of differentially expressed lncRNA and miRNA, as well as the interactions of miRNA and ferroptosis-related mRNA, were predicted utilizing TargetScan (www.targetscan.org), StarBase (https://starbase.sysu.edu.cn) and RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid) databases. The intersection between the 3 databases was identified using the “vennR” package of R software. Then, the ceRNA network was established using the common miRNAs in the 3 databases that connect lncRNAs and mRNAs. Following the mechanism of the ceRNA network, we only retained the relationships as upregulated lncRNA\downregulated miRNA\upregulated mRNA and downregulated lncRNA\upregulated miRNA\downregulated mRNA, which were used to construct the ceRNA network, and then the results were imported into Cytoscape 3.7.1 (https://cytoscape.org) for visualization.

Identification and validation of a prognostic lncRNA signature based on the ceRNA network

Using the survival data in TCGA and excluding the incomplete survival time data, we constructed the prognostic signature using lncRNAs in ceRNA. KIRC patients were randomly assigned into the training and test cohorts in a 1:1 ratio, utilizing the “caret” package in R 4.1.1. In the training set, univariate Cox regression analysis of overall survival (OS) was performed to explore the ferroptosis-related lncRNAs with prognostic values (p < 0.01) in the ceRNA network. Then, to minimize the risk of overfitting, least absolute shrinkage and selection operator (LASSO) regression analysis was utilized to construct a prognostic signature in the R package “glmnet”. The following formula:

risk score = (β1×G1+β2×G2+β3×G3+...+βn×Gn)

was used to calculate the risk score for each patient and to predict the prognosis of the patient according to the score, where β stands for the coefficient of each lncRNA, G means each lncRNA expression value, and n represents the number of lncRNAs.

Based on the median risk score of the training set, we stratified patients into 2 groups. The “Pheatmap” package was applied to draw a scatter diagram to describe the distribution pattern of risk scores and the corresponding survival time of each patient in the training and validation set. The “Stats” package was used for principal component analysis (PCA) to describe the gene expression distribution in the signature. Then, the “Rtsne” package was used for t-distributed stochastic neighbor embedding (t-SNE) analysis to describe the distribution of survival status in different risk groups. Receiver operating characteristic (ROC) curves, drawn by the “timeROC” package, were used to verify whether the signature could be considered a useful biomarker in the training and validation set. Finally, the “survival” package was used for survival analysis for each lncRNA in the signature to verify its prognostic value in the KIRC cohort. The mRNAs competitively inhibited by lncRNAs in the ceRNA network were verified in the Gene Expression Profiling Interactive Analysis (GEPIA; http://gepia.cancer-pku.cn/) database (based on the TCGA and Genotype-Tissue Expression (GTEx) data) in ccRCC.

Independent prognostic analysis of the ferroptosis-related lncRNA signature

To further estimate the independent prognostic value of the ferroptosis-related lncRNA signature, univariate and multivariate Cox regression analyses were used to find out whether it was influenced by other clinical features. Available clinical characteristics including sex, age and tumor-node-metastasis (TNM)-based staging were converted into dichotomous variables and used to calculate OS-based hazard ratios (HRs) and 95% confidence intervals (95% CI). A p < 0.05 was considered statistical significance.

Functional enrichment analysis

The “limma” R package was performed to select the risk-related differentially expressed genes (DEGs) in the KIRC cohort based on | log2 Fold Change |≥1.2 and false discovery rate <0.05 as a standard between the high- and low-risk groups. The “clusterProfiler” R package was applied to conduct Gene Ontology (GO; https://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/) analyses for risk-related DEGs based on the criterion: gene count >5 and p-value < 0.05. The infiltrating scores of 16 immune cells and the activities of 13 immune-related pathways were quantified with single-sample gene set enrichment analysis (ssGSEA) using the “GSVA” R package. Using the annotated gene set provided in Supplementary Table 2, we quantified the immune infiltration enrichment scores of different immune cells and immune-related functions to further study the correlation between risk scores and immune status.

Statistical analyses

Statistical analysis was performed using R software (R v. 4.1.1; R Foundation for Statistical Computing, Vienna, Austria). The key R-scripts used in this study were presented in the Supplementary material – R-scripts. Continuous variables were analyzed with the Wilcoxon test (Mann–Whitney test), whereas categorical data was analyzed using the χ² test or Fisher’s exact test. Univariate and multivariate Cox regression analysis was used to identify the related factors affecting the OS of KIRC patients. The results of the proportional hazards assumption for Cox regression were conducted based on Schoenfeld’s global and individual test (Supplementary Table 3, Supplementary Fig. 2,3). Kaplan–Meier analysis and the log-rank test were performed to calculate the survival difference. For all statistical analyses, a p-value < 0.05 was considered statistically significant without special explanation.

Results

Identification of DEGs in the KIRC cohort

After filtering the sample types of Additional-NewPrimary (1 transcriptome and 1 miRNA) and duplicate samples (8 transcriptomes and 28 miRNAs), a total of 583 transcriptome samples, and 583 miRNA samples, including 512 KIRC patients and 71 normal controls were obtained from the KIRC cohort of TCGA database. Through differential expression analysis, we screened out 356 lncRNAs (286 upregulated and 70 downregulated), 132 miRNAs (62 upregulated and 70 downregulated) and 35 ferroptosis-related mRNAs (18 upregulated and 17 downregulated) (Figure 1A–C).

Establishment of a ferroptosis-related ceRNA network in the KIRC cohort

To understand the regulatory crosstalk between different RNA molecules, we constructed a ferroptosis-related ceRNA network in the KIRC cohort. RNAHybrid, StarBase and TargetScan databases were utilized to predict the network of lncRNAs, miRNAs and mRNAs. In total, 8658,1263,160426 lncRNA-miRNA interactions and 1594,913,10067 miRNA-ferroptosis-related mRNA interactions were predicted, respectively (Figure 1D,E). Taking miRNA as the connecting point, we constructed a ceRNA network with 182 axes, according to the expression of genes in the network. The intersection of target miRNAs for downregulated lncRNAs, downregulated mRNAs with upregulated miRNAs, and target miRNAs for upregulated lncRNAs, upregulated mRNAs with downregulated miRNAs were selected and finally, a ceRNA network with 77 axes was obtained (Figure 1F). Among them, there were 27 lncRNAs (26 upregulated and 1 downregulated), 13 miRNAs (12 downregulated and 1 upregulated), and 9 mRNAs (8 upregulated and 1 downregulated). The regulatory relationships disclosed by these networks may provide an idea for exploring the molecular mechanism of the ferroptosis-related lncRNAs.

Construction of lncRNA signature based on the ferroptosis-related ceRNA network

Given the potential of ferroptosis-related ceRNA in regulating patient prognosis, we focused on 27 lncRNAs within the ceRNA network. A total of 508 eligible ccRCC patients (4 samples with missing survival time were removed) with integrated information were incorporated in the TCGA-KIRC dataset and randomly assigned to the train set (254 samples) and validation set (254 samples). Primarily, 27 lncRNAs in the ferroptosis-related ceRNA network were extracted to identify the prognostic risk model. Then, we performed the Cox proportional hazard assumption test based on Schoenfeld’s residuals. All the lncRNAs in the univariate Cox regression model satisfied the proportional hazards assumptions (p > 0.05) (Supplementary Table 3). The Schoenfeld’s residual plots did not suggest a violation of the proportional hazard assumption (Supplementary Fig. 2). Next, 14 lncRNAs significantly related to OS were selected in the univariate Cox regression analysis, which was considered as potential predictors. Finally, a LASSO regression algorithm was performed for feature selection. When the partial likelihood binomial deviation reaches the minimum value, the most appropriate tuning parameter λ for LASSO regression is 0.055. According to the penalty parameter (Lambda) in the model, we established a prognostic signature of KIRC patients consisting of 6-ferroptosis related lncRNAs including PVT1, CYTOR, MIAT, SNHG17, LINC00265, and LINC00894 (Figure 2A–C). A risk score for each patient was calculated according to the following risk formula:

risk score = 
β × expression value of PVT1
+ β × expression value of CYTOR
+ β × expression value of MIAT
+ β × expression value of SNHG17
+ β × expression value of LINC00265
+ β × expression value of LINC00894

We evenly categorized patients into a high-risk group (n = 127) and a low-risk group (n = 127) based on risk scores (Figure 3A). As depicted in Figure 3B,C, PCA analysis and t-SNE analysis were used to reduce the dimension of the data and observe the significant difference between the 2 groups. Our results indicated that patients in the high- and low-risk groups showed a significant 2-way distribution. As displayed in Figure 3D, there were significantly more deaths in high-risk patients than in low-risk ones. Additionally, the Kaplan–Meier curve revealed that patients in the high-risk group have corresponded to a worse OS in the training set (Figure 4A, logrank test, χ² = 32.6, p = 1.136e−08). Time-dependent ROC curves were performed to evaluate the predictive accuracy of the risk score for OS, and the area under the curve (AUC) was 0.791, 0.740 and 0.890 in predicting 1-, 5- and 10-years OS in KIRC, respectively (Figure 4B). More significantly, the signature shows superior prognostic prediction efficiency than tumor-node-metastasis (TNM) staging with the AUC of risk score of 0.78 on multi-factor ROC (Figure 4C).

Validation of the 6-lncRNA signature in the validation set

Subsequently, we performed the same analysis in the validation set. Once more, 254 patients in the validation set were divided into high-risk (n = 127) and low-risk (n = 127) groups based on the median risk values in the train set (Figure 5A). As predicted, the results of PCA analysis and t-SNE analysis indicated that the patients of the validation set in the high- and low-risk group also showed 2-way distribution, and patients in the high-risk group have corresponded to more death cases (Figure 5B–D). Moreover, like the training set, the Kaplan–Meier curve showed that the mortality of patients in the high-risk group was significantly higher than that in the low-risk group (Figure 6A, log-rank test, χ² = 17.9, p = 2.371e–5). The AUC reached 0.723, 0.714 and 0.780 in predicting 1-, 5- and 10-year OS, and the risk score in muti-ROC was 0.721 (Figure 6B,C). These results indicated that the prognostic signature can be regarded as a qualified prognostic evaluation.

Survival analysis, univariate Cox regression analysis and expression comparative analysis for the signature in the ceRNA network

Moreover, to further evaluate the signature, survival analysis was performed for each lncRNA in the signature (PVT1, CYTOR, MIAT, SNHG17, LINC00265, and LINC00894). As shown in Figure 7, the results showed that all 6 lncRNAs were significantly negatively correlated with OS. Then, to explore the regulatory relationship of the signature on ferroptosis-related genes, Cytoscape 3.7.1 software was used to show the specific regulation relationships of 6 lncRNAs in the signature and 7 ferroptosis-related genes (Figure 8A). By comparing ccRCC samples with normal kid­ney samples in TCGA and GTEx databases, we analyzed the expression of 7 lncRNA-regulated mRNAs (includ­ing CD44, PML, TAZ, CDKN2A, SCD, MYB, and CAV1) in the ceRNA network. The cutoff value was set to |log2 Fold Change | > 1, p < 0.05, and the results showed that all 7 ferrop­tosis-related genes were highly expressed in ccRCC, among which 5 were statistically significant, except TAZ and MYB (Figure 8B). Besides, in the FerrDb database, we discovered that these 7 genes all regulate ferroptosis with the help of lipid reactive oxygen species (ROS), except PML, and could be divided into 2 categories: ferroptosis suppressors (CAV1, CD44, PML, and SCD) and ferroptosis drivers (CDKN2A, TAZ and MYB) (Supplementary Fig. 1).

Prognostic independence analysis of ferroptosis-related lncRNA signature and clinical features

Univariate and multivariate Cox regression analyses were performed among the available variables (including age, gender and TNM staging) to determine whether the signature was an independent prognostic predictor for OS. All the covariates in the multivariate Cox regression model satisfied the proportional hazards assumptions and the global Schoenfeld’s test (Supplementary Fig. 3). Univariate Cox regression analysis showed that both in training set and validation set, TNM staging (training set: HR = 1.946, 95% CI = 1.604–2.361, Wald test, W = 45.650, p < 0.001; validation set: HR = 1.858, 95% CI = 1.533–2.252, Wald test, W = 39.900, p < 0.001) and risk signature (training set: HR = 3.208, 95% CI = 2.199–4.678, Wald test, W = 37.620, p < 0.001; validation set: HR = 3.600, 95% CI = 2.274–5.698, Wald test, W = 29.880, p < 0.001) were significantly correlated with OS (Figure 9A,C). After independent adjustment for other clinical features in multivariate Cox regression, risk signature remained a reliable prognostic predictor of OS in both cohorts (train set: HR = 2.612, 95% CI = 1.793–3.805, Wald test, W = 25.000, p < 0.001; validation set: HR = 2.416, 95% CI = 1.522–3.835, Wald test, W = 14.200, p < 0.001) (Figure 9B,D). The detailed Cox regression results can be found in Supplementary Table 4. Moreover, we examined the linear relationship between the log hazard and the continuous variable age. The results showed an overall linear trend in the relationship between the continuous covariate age and log hazard (Supplementary Fig. 4).

GO and KEGG enrichment analyses

To evaluate the biological functions and pathways associated with the prognostic signature, GO enrichment and KEGG pathway analyses were performed using DEGs between high- and low-risk groups. In total, 1,168 DEGs were significantly enriched in 367 biological processes, 51 cellular components, 52 molecular functions, and 10 KEGG entries. The top 10 enriched items of GO and KEGG were selected and demonstrated. As expected, GO analysis revealed that there were several molecular function entries related to lipid metabolism: these functions may be related to lipid peroxidation of ferroptosis. The KEGG analysis indicated that glycerophospholipid metabolism and phospholipase D signaling pathway were enriched. In addition, some immune-related items were also found to be significantly enriched in the 2 cohorts, including T cell activation, T cell differentiation, regulation of T cell activation, positive regulation of leukocyte cell adhesion, positive regulation of T cell activation, regulation of T cell differentiation, regulation of lymphocyte differentiation positive regulation of cell-cell adhesion, regulation of leucocyte cell adhesion, lymphocyte differentiation, cytokine receptor activity in molecular functions, immune receptor activity, cytokine binding, and so on (Figure 10A–D).

ssGSEA and ESTIMATE analysis of immune infiltration

The results from GO and KEGG suggest that the ferroptosis-related signature may be involved in the regulation of the immune microenvironment. Therefore, we explored its potential role in the TME. The ssGSEA approach was utilized to further evaluate the correlation between the risk score and immune cell infiltration status, and the results revealed that 9 of the 16 immune cells were found to have significantly high infiltration enrichment scores in the high-risk group, including CD8⁺ T cells, pDCs, T helper cells, Th1 cells, TIL, aDCs, B cells, and Th2 cells (Figure 11A). In terms of immune function, except for type II immune interferon response, the high-risk group was positively correlated with other immune-related functions (Figure 11B). Then, we performed the ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) algorithm to further analyze the differences in immune components of TME in the high- and low-risk groups, and our results indicated that the ESTIMATE score and immune score were increased significantly in the low-risk group (Wilcoxon rank sum test, W = 38063, p < 0.001; Wilcoxon rank sum test, W = 41373, p < 0.001) (Figure 11C), while stromal score was no significant difference between the 2 groups (Wilcoxon rank sum test, W = 31739, p > 0.05) (Figure 11C).

Discussion

Clear cell renal cell carcinoma displays considerable heterogeneity in both its molecular characteristics and biological behavior, thereby presenting substantial challenges for clinicians involved in the treatment of cancer patients.¹⁹ In clinical practice, grading and staging are commonly used clinical pathological parameters to assess the prognosis of ccRCC and guide treatment decisions.²⁰ However, patients with the same cancer staging often exhibit significant differences in prognosis and response to treatment.²¹ Specifically, conventional histopathological measures and similar treatment strategies may result in curing patients, while others may die prematurely. This indicates the limitations of clinical pathological parameters in describing the survival outcomes of ccRCC patients.²¹ Therefore, it is necessary to surpass the prognostic features of the current staging system in order to accurately identify patients who may develop aggressive diseases with poor survival rates and provide better guidance for clinical therapy. In the current study, we first identified lncRNAs, miRNAs and ferroptosis-related mRNAs through differential expression analysis of 583 samples downloaded from the KIRC cohort of the TCGA database. Then, through RNAHybrid, StarBase and TargetScan databases, 27 lncRNAs, 13 miRNAs and 9 mRNAs conforming to the ceRNA mechanism were predicted and selected. Through a series of analyses, including univariate Cox regression analysis, LASSO regression analysis and multivariate COX regression analysis, 6 lncRNAs significantly associated with OS were identified and used to establish a ferroptosis-related signature. To explore the distribution of each case with different risk values, PCA and t-SNE analyses were performed in the training and validation cohorts, and our results showed that KIRC patients with low-risk scores could be distinguished from those with high-risk scores. Moreover, the death probability distribution analysis also found that the number of deaths in KIRC patients with high-risk scores was significantly higher than that in KIRC patients with low-risk scores. Through Kaplan–Meier curve analysis and Cox regression analysis, we identified that the ferroptosis-related signature could stratify the risk for ccRCC patients based on OS and serve as an independent prognostic factor for clinical outcomes. Furthermore, the ROC curve analysis demonstrated that the ferroptosis-related signature could accurately predict short- and long-term survival in ccRCC patients, with significantly higher predictive accuracy compared to most clinical features. Taken together, these results indicate that the ferroptosis-related signature could be a stable and robust tool for clinicians to predict the prognosis of ccRCC patients.

Within the ferroptosis-related signature that we constructed, a higher risk indicated a poor prognosis. To further understand the underlying molecular mechanisms, we initially analyzed the effects of 6 lncRNAs on OS in different risk subgroups. Not surprisingly, the 6 lncRNAs (SNHG17, CYTOR, LINC00894, LINC00265, MIAT, and PVT1) included in the signature were significantly associated with the unfavorable prognosis of ccRCC. Among the 6 lncRNAs in the signature, some have previously been reported to be associated with ccRCC. Research has confirmed that lncRNA PVT1 is significantly overexpressed and can act as a ceRNA in the context of RCC.²² Additionally, Yan et al. identified that lncRNA MLAT is upregulated in ccRCC tissues and cell lines and can act as a sponge for miR-29c, increasing the expression of ZEB1.²³ Recently, a study identified that SNHG17 was significantly upregulated in ccRCC and could be used as a prognosis predictor.²⁴ Deng et al. uncovered that LINC00894 was upregulated and significantly related to the prognosis in ccRCC patients.²⁵ Also, the function of CYTOR has been identified to be an oncogene in many cancers. However, as well as the LINC00265, it has not been elucidated in the previous studies and our research uncovered their potential roles in KIRC. To further investigate the role of ferroptosis-related lncRNAs in gene regulation in ccRCC, we constructed a ceRNA network. According to the regulatory network in the FerrDb database, all 7 genes were somehow connected with lipid ROS to regulate ferroptosis, except that PML was a straight ferroptosis suppressor. The high lipid ROS connection was in accord with the lipid metabolism relevancy found in GO analysis, which revealed that the signature was highly associated with the lipid peroxidation process in ferroptosis. Moreover, we uncovered that CAV1, CD44, PML, and SCD could act as a suppressor of ferroptosis, while CDKN2A, MYB and TAZ were drivers of ferroptosis, which may explain the double-sided effect of ferroptosis on tumor inhibition and tumor growth. CAV1, CD44, PML, and SCD have been confirmed to be highly expressed in tumor tissues and promote tumor progression by inhibiting ferroptosis.^{26, 27, 28, 29} CDKN2A has been identified to play a vital role in tumorigenesis by enhancing p53-dependent transactivation and ferroptosis,³⁰ the hypermethylation of CDKN2A might be a predictor of poor prognosis for cancer.³¹ MYB, a vital transcription factor in solid tumors, can regulate the expression of CDO1; then, CDO1 can convert cysteine into taurine, reduce the utilization of cysteine, limit the synthesis of glutathione, inhibit the antioxidant capacity of cells, and eventually cause ferroptosis.³² A recent study showed that TAZ is highly expressed in RCC and can mediate cell density-regulated ferroptosis.¹⁴ These findings revealed potential regulatory mechanisms of ferroptosis-related molecules in the occurrence and development of ccRCC, providing potential biomarkers for the identification and personalized treatment of this cancer.

To better understand the transcriptional regulatory mechanisms of ferroptosis-related lncRNAs in ccRCC, we conducted a comprehensive exploration of their functions and pathways involved. The GO functional enrichment analysis revealed that these differentially expressed genes were involved in lipid metabolism-related functions, which may be associated with lipid peroxidation of ferroptosis, and the KEGG analysis revealed that several signaling pathways related to cancer and ferroptosis were enriched, such as glycerophospholipid metabolism, phospholipase D signaling pathway and lipase activity pathway. As previously reported, these pathways have been found to be abnormally overactivated in multiple types of cancers, which drives the proliferation, invasion and metastasis of cancer cells, and is often associated with poor clinical prognosis.^{33, 34, 35} Hence, these oncogenic pathways may explain the impact of the ferroptosis-related signature on the prognosis of ccRCC patients. Drugs targeting these abnormally overactivated pathways may have the potential to inhibit the progression of this cancer.

It has been reported that the mechanism by which ferroptosis cells trigger potent immune responses may have some similarities with traditional immunogenic cell death.³⁶ Considering the close relationship between ferroptosis and tumor immunity, we utilized 2 algorithms (ESTIMATE and ssGSEA) to investigate the immune microenvironment associated with the ferroptosis-related signature in ccRCC. In the TME component analysis, the high-risk group had higher ESTIMATE scores and immune scores. Some studies have demonstrated that the abundance of immune cell components serves as an independent prognostic factor that is crucial to the prognosis of ccRCC.^{37, 38} Therefore, we hypothesized that the influence of 6 lncRNAs (SNHG17, CYTOR, LINC00894, LINC00265, MIAT, and PVT1) on unfavorable survival among patients with ccRCC might be associated with the remodeling of immune components within the TME. This finding provides indirect evidence for the close association between the ferroptosis-related signature and the survival of ccRCC patients. The analysis of TME cell types from ssGSEA also fully supported this hypothesis. High-risk individuals had an active presence of immune cells, such as CD8⁺ T cells and Th1 cells, indicating enhanced anti-tumor immune activity. Overall, these discoveries provided novel insights into the mechanisms of signature lncRNAs regulating the immune microenvironment of ccRCC patients.

Limitations

Our study has some limitations. For instance, the ferroptosis-lncRNA signature was developed according to the TCGA public database, but it would be better to have strong external data to confirm its validity and practicality. Furthermore, our study is mainly based on integrated bioinformatics and lacks confirmation of these findings from valid clinical studies.

Conclusions

In the present study, we developed a ferroptosis-related lncRNA signature exhibiting high accuracy and stability, which shows potential as a prognostic prediction tool for ccRCC patients. Moreover, by constructing a ceRNA network and conducting immune infiltration analysis, we elucidated the potential mechanisms by which ferroptosis-related lncRNAs regulate ccRCC, and explored their role within the TME, offering novel insights for precise treatment based on ferroptosis signaling. Our study provides a potential option for prognosis prediction and personalized treatment in ccRCC patients.

Data availability

The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate

This article did not involve any studies with human or animal participants conducted by the authors, and therefore, ethical approval or consent was not required. The study did not require any administrative permission or licenses to access the original data used in the research.

Supplementary data

The Supplementary materials are available at https://doi.org/10.5281/zenodo.8032200. The package includes the following files:

Supplementary Fig. 1. Seven ferroptosis-related genes were divided into 2 categories in the FerrDb database. Ferroptosis suppressor: CAV1, CD44, PML and SCD. Ferroptosis driver: CDKN2A, TAZ and MYB.

Supplementary Fig. 2. The results of the proportional hazards assumption (lncRNAs selection using univariate Cox regression) based on Schoenfeld’s global and individual test.

Supplementary Fig. 3. The results of the proportional hazards assumption (selection of lncRNA signature and clinical features) based on Schoenfeld’s global and individual test.

Supplementary Fig. 4. The linear relationship between the log hazard and the continuous variable age.

Supplementary Table 1. Ferroptosis-related genes (n = 112) validated in humans and protein coding were obtained from the FerrDb database.

Supplementary Table 2. Annotated gene sets to quantify the immune infiltration enrichment scores of different immune cells and immune-related functions.

Supplementary Table 3. The results of the proportional hazards assumption for Cox regression based on Schoenfeld’s global and individual test.

Supplementary Table 4. The detailed Cox regression results, including HR, coefficient, 95% CI, Wald test, and likelihood ratio (LR) test