Comprehensive analysis of new prognostic signature based on ferroptosis-related genes in clear cell renal cell carcinoma

The expression pattern and correlation of ferroptosis genes in ccRCC

A heatmap was generated to analyze the expression pattern of FRGs in ccRCC. Most of the genes (52/64, 81.3%) showed a significantly different expression level in the ccRCC tissues compared to the normal tissues (Figure 1A). Furthermore, univariate Cox analysis revealed that 35 out of the 52 DEGs were significantly associated with OS (Figure 1B). Next, a Venn diagram was constructed to screen out the prognostic ferroptosis-related DEGs, with results showing that 27 DEGs were associated with prognosis (Figure 1C). Finally, the 27 prognostic FRGs were preserved (p < 0.05, Figure 1D, 1E).

Moreover, correlation analysis was used to investigate the interactions among all selected genes, with results suggesting that most of the FRGs had a positive correlation (Figure 1F). The PPI network indicated that SLC7A11, GCLM, and NFE2L2 were the hub genes (Figure 1G).

Construction of prognostic signature in the TCGA cohort

The entire TCGA cohort was randomly divided into train set (n = 264) and test set (n = 261) using the “caret” package. The 27 DEGs mentioned above were then measured as predictive genes and exposed to LASSO analysis. From the results, six FRGs were screened out based on the optimal value of λ (Supplementary Figure 1A, 1B). Subsequently, multivariate Cox analysis was performed, and a four genes prognostic model was finally constructed (Figure 2A). The risk score was then determined using the coefficients and expression level of each gene: risk score = (0.015) × CD44 + (-0.005) × DPP4 + (-0.017) × NCOA4 + (0.432) × SLC7A11. According to the median risk score, patients in the train set were grouped into high- and low risk sets, and followed by comparison of the OS using the K-M curve. Results showed that ccRCC patients in the high-risk group had a poorer OS compared to patients in the low-risk group (Figure 2B). The distribution of the four FRGs signature based on risk scores is also displayed in Figure 2C, 2D, which was consistent with results of the K-M curve. Furthermore, time-dependent ROC analysis demonstrated that the prognostic signature had an advanced predictive performance for OS, with AUC values equal to 0.756, 0.753 and 0.769 at one, two and three years, respectively (Figure 2E). Finally, PCA and t-SNE analyses revealed that the two groups of patients were distributed in two different directions (Figure 2F, 2G).

Validation of the four-gene signature

According to the median risk score from train group, patients in test group were also categorized into high- and low-risk. Results showed that the survival outcome of high-risk patients was worse than that of low-risk patients, which was consistent with results of the train set (Figure 3A). The distribution of risk scores confirms that patients with high-risk scores have poorer survival outcomes compared to patients with low-risk scores (Figure 3B, 3C). The ROC curve demonstrated that the four-gene signature had a preferable predictive capacity for OS, with AUC values of 0.718, 0.644 and 0.670 at one, two and three years, respectively (Figure 3D). In addition, PCA and t-SNE analyses confirmed that the two groups of patients with ccRCC were dispensed in two directions (Figure 3E, 3F).

Furthermore, the robustness of the four-gene model was verified for its clinical application using KIRC samples obtained from our medical center. IHC images showed that the normal renal tissue had weak staining for CD44, DPP4 and SLC7A11 in the cytoplasm and cell membrane, while NCOA4 displayed the opposite trend in normal tissues (Figure 4A, 4C, 4E). On the other hand, moderate and strong staining patterns for CD44, DPP4, and SLC7A11 were observed in the cytoplasm and cell membrane of tumor tissues (Figure 4B, 4D, 4F), while NCOA4 only exhibited weak positive staining on the cell membrane of tumor tissues (Figure 4G, 4H). These unique IHC staining patterns illustrate that the four genes can be used to predict clinical outcome and can distinguish cancerous tissue from normal tissue.

Subsequent qRT-PCR analysis also showed that the expression level of CD44, DPP4 and SLC7A11 in ccRCC samples are significantly higher than that in paired non-tumor samples (Figure 4I). In addition, the expression level of NCOA4 in ccRCC samples is lower than that in paired non-tumor samples (p < 0.05, Figure 4I).

Independent prognostic value of the constructed signature in ccRCC

To determine whether the four-gene signature could serve as independent prognostic factor, univariate and multivariate Cox regression analyses were operated on both the train and test sets. Univariate Cox analysis indicated that age, tumor grade, tumor stage and risk score were associated with OS in both train and test sets (p < 0.05, Figure 5A, 5C). The variables with associated P value <0.1 were then enlisted into the multivariate Cox analysis, which showed that the predictors were risk score and tumor stage in both train and test sets (p < 0.05, Figure 5B, 5D).

Functional analyses

To illustrate the biological functions and pathways of genes that were differentially expressed in ccRCC tissues compared with non-cancerous tissues, GO enrichment and KEGG pathway analyses were conducted. GO analysis suggested that the DEGs were enriched in peptidase regulator activity, signaling receptor activator activity, receptor ligand activity, and plasma membrane (FDR < 0.05, Figure 6A), which is consistent with previous studies [18–21]. On the other hand, KEGG pathway analyses implied that the DEGs were enriched in the TNF signaling pathway, PPAR signaling pathway, p53 signaling pathway, folate biosynthesis, bile secretion, and cholesterol metabolism (Figure 6B), which is congruent with preceding studies [22–27]. However, the DEGs were also enriched in several immune functions and immune-associated pathways, such as cytokine-cytokine receptor interaction, lymphocyte mediated immunity, and complement activation, complement and coagulation cascades (FDR < 0.05, Figure 6A, 6B).

Considering the association between different DEGs and immune status, we further explored the interaction between risk score and immune status in the TCGA cohort. Major immune cells and functions were quantified with enrichment scores using ssGSEA. Results indicated that the scores of antigen-presenting cells (APC) and functions in the high-risk group, including DCs, iDCs and APC co-inhibition, were higher than in the low-risk group (p < 0.05, Figure 6C, 6D). However, the score of aDCs was lower in the high-risk group. Notably, cellular immunity, especially T cell-related immune cells and functions, including CD8⁺ T cell, Th1 cell, Th2 cell, Tfh cell and T cell co-stimulation, had higher scores in the high-risk group (p < 0.05, Figure 6C, 6D). Furthermore, consistent with the KEGG analysis, quantification of cytokine-cytokine receptor (CCR) interaction produced higher scores in the high-risk set (p < 0.05, Figure 6D). In addition, tumor infiltrating lymphocytes (TIL), the function of checkpoint regulators, and the ability of inflammation-promoting molecules displayed elevated scores in the high-risk group (p < 0.05, Figure 6C, 6D).

Finally, we performed co-expression analysis between four ferroptosis genes through Person’s correlation analysis. As shown in Supplementary Figure 2, based on the optimal correlation, we selected NCOA4 and SLC7A11 as the hub genes. Then, GSEA analysis was conducted to explore the underlying mechanisms based on the TCGA data. Notably, the results showed that overexpression of NCOA4 was associated with renal cell carcinoma, the JAK-STAT pathway, Toll-like pathway and the NOD-like receptor pathway (|NES|> 1, p < 0.05, FDR < 0.25, Supplementary Figure 3A–3D). In addition, overexpressed SLC7A11 enriched in cell cycle and NOD-like receptor pathway as well (|NES|> 1, p < 0.05, FDR < 0.25, Supplementary Figure 3E, 3F).

Data collection

Expression data and clinical information of 539 ccRCC and 72 normal patients were obtained from the TCGA database (https://portal.gdc.cancer.gov/). Expression data were normalized by “limma” package. After performing a comprehensive literature review, 64 FRGs were identified [13, 19, 38–41] (Supplementary Table 1).

Selection of intersected genes between differential expression- and prognostic ferroptosis-related genes

DEGs in paired tumor and non-tumor tissues were assessed by “limma” package. Univariate Cox analysis of overall survival (OS) was then conducted to obtain prognostic genes. Finally, intersected FRGs were screen out using the “venn” R package. The differential expression patterns, co-expression correlation analysis, and prognostic values were visualized by “survival”, “igraph” and “pheatmap” R package.

Construction and validation of the prognostic signatures

The TCGA cohort was randomly divided into two groups (train and test groups) by using the “caret” package. Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was conducted in train set with the "glmnet" R package. Multivariate cox analysis was performed to build prognostic models. The median value of risk scores was calculated to stratify ccRCC patients into high-risk and low-risk sets by the following formula: $\text{Risk}\text{score}={\sum }_1^{n}Coe{f}_n\times x_n$ (where Coef_n is the coefficient and x_n is the expression level of each intersected gene). Next, a Kaplan-Meier (K–M) curve was used to analyze the OS in train and test groups, respectively. Receiver operating characteristic (ROC) analysis was also executed to assess the prediction efficiency of prognosis. Based on the expression profile of four FRGs, PCA and t-SNE analysis were performed to demonstrate the expression patterns of various FRGs in diverse groups. Finally, IHC pattern was certified by utilizing KIRC clinical samples from our hospital. Univariate and multivariate Cox analyses were utilized to investigate independent prognostic factors for ccRCC patients.

Functional analyses

A protein–protein interaction (PPI) network was built for prognostic DEGs by the STRING database (http://string-db.org/). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on the DEGs were then operated by utilizing "clusterProfiler" package. In addition, the score of major immune cells and immune-related pathways were determined with single-sample gene set enrichment analysis (ssGSEA). GSEA was performed using GSEA v4.1.0 (http://www.broadinstitute.org/gsea/).

Immunohistochemistry

Ten pairs of ccRCC and adjacent normal tissues were collected from Shandong provincial Hospital affiliated to Shandong First University during January 2021 to February 2021. The experiment was approved by the Shandong Provincial Hospital Ethics Committee (Approval number: SWYX: NO.2021-118) and written informed consents were signed by each patient before the study began. IHC was performed according to previously described standard procedures [42]. All samples were incubated with rabbit polyclonal antiCD44 (ab189524), antiNCOA4 (ab62495), antiDPP4 (ab187048) and antiSLC7A11 (ab175186) antibodies overnight at 4° C and were then washed. Two pathologists independently assessed the IHC slides.

RNA analysis, extraction, and quantitative real-time PCR

Total RNA of frozen tissue was extracted by TRIzol reagent (Tiangen Biotech (Beijing)), and 1 μg of total RNA was reverse transcribed using the FastKing RT Kit (Tiangen Biotech (Beijing)) according to the manufacturer’s instructions. Followed by measuring with a real-time quantitative PCR system. The primers used in this study were provided by Beijing Dingguo Changsheng Company and are shown in Supplementary Table 2. The expression data was log2 transformed: log2(exp + 0.01).

Statistical analysis

Mann-Whitney tests were utilized to measure gene expression between ccRCC and non-tumor tissues. Associated samples with incomplete clinical information were eliminated. All data analyses were conducted using the R statistical package (R version 4.0.1). A two-tailed p < 0.05 was considered statistically significant.

Data availability statement

The data used to support the findings of this study are available from the corresponding author upon request.