Differentially expressed FRGs in KIRC
Using transcriptome data from the TCGA database, we analyzed the mRNA levels of 61 FRGs in KIRC and normal tissues. A total of 19 differentially expressed FRGs were eventually identified using the criteria of |log2FC| > 1 and FDR < 0.05, including 13 downregulated genes (MT1G, ACSF2, CHAC1, ACSL4, AKR1C2, PEBP1, PTGS2, AKR1C1, CBS, GOT1, ACO1, FDFT1, and HMGCR) and 6 upregulated genes (ALOX12, CD44, SLC7A11, ALOX5, HMOX1, and ALOX15B) in KIRC tissues compared to the normal kidney tissues (Fig. 1a, b). The box diagram was utilized to exhibit the expression patterns, median values, and data ranges of the differentially expressed FRGs in tumor and normal samples (Fig. 1c). The interaction network among these genes was presented in Fig. 1d, and the result indicated that PTGS2 and HMOX1 seemed to be the hub genes in this network.
Two ferroptosis subgroups were different in clinical phenotypes and OS via consensus clustering analysis
According to expression levels of 61 FRGs, the total TCGA-KIRC cohort was clustered into 2 subgroups (cluster 1 and cluster 2) with k = 2 as the optimal value because the grouping was suboptimal when they were divided into more than 2 clusters (Fig. 2a-c). Moreover, principal component analysis (PCA) was performed to compare the transcriptional profile between cluster 1 and cluster 2. The result demonstrated that there was a significant distinction between the two subgroups (Fig. 2d). To better understand the clustering result and its relationships with survival outcomes, we compared the OS between cluster 1 and cluster 2 and observed cluster 2 had shorter OS than cluster 1 for the KIRC patients (Fig. 2e). We then evaluated associations between the clustering and the clinicopathological parameters of KIRC patients. The result showed that these 2 clusters were different in grade (P < 0.05), AJCC stage (P < 0.001), T stage (P <0.01), N stage (P <0.01), and survival status (P <0.05), but did not show any significant differences in age and gender (Fig. 2f). Therefore, these results suggested that ferroptosis was closely related to clinical phenotypes and the progression of KIRC.
Functional enrichment analyses of the differentially expressed FRGs
To elucidate the biological functions and pathways of the 19 differentially expressed FRGs, GO functional annotation and KEGG pathway enrichment analyses were conducted. The results indicated that these differently expressed FRGs were significantly enriched in the biological process (BP) related to several metabolic process, such as cofactor metabolic process, fatty acid metabolic process, and fatty acid derivative metabolic process (Fig. 3a, b). Furthermore, lipoxygenase pathway was also involved. In terms of molecular function (CC), we found that the differently expressed FRGs were significantly enriched in peroxisomal membrane, microbody membrane, and caveola. Through the molecular function (MF), the differently expressed FRGs were notably enriched in oxidoreductase activity, dioxygenase activity, and lyase activity. In the KEGG pathway enrichment analysis, these genes were shown to be mostly associated with pathways in arachidonic acid metabolism, serotonergic synapse, and ferroptosis (Fig. 3c, d). In addition, we explored the effect of the differentially expressed FRGs in multiple classical signaling pathways on KIRC using the GSLA database. The results revealed that some FRGs was associated with the activation or inhibition of oncogenic pathways. For example, SLC7A11 expression was related to the activation of apoptotic, cell cycle, and EMT pathways; PTGS2 expression was related to the inhibition of cell cycle, DNA damage response, hormone AR, and hormone ER (Additional file 2: Figure S1).
Prognosis-related FRDs selecting and construction of a prognostic signature based on seven FRGs in the TCGA training cohort of KIRC patients
We conducted a univariate Cox regression analysis on the transcriptome data from the TCGA-KIRC dataset and found that the expression of 11 FRGs (CBS, GOT1, FDFT1, HMOX1, CD44, ACO1, AKR1C2, PEBP1, CHAC1, HMGCR, and SLC7A11) was significantly related to the KIRC survival (P < 0.05; Fig. 4a). Genes (CBS, CD44, AKR1C2, CHAC1, and SLC7A11) with HR > 1 were considered as risk genes, while the remaining six genes (GOT1, FDFT1, HMOX1, ACO1, PEBP1, and HMGCR) with HR < 1 as protective genes. Based on the expression profile of the 11 genes mentioned above in the training cohort, a prognostic risk signature was constructed by the least absolute shrinkage and selection operator (LASSO) Cox regression analysis. As a result, a 7-gene signature (CBS, HMOX1, CD44, AKR1C2, CHAC1, HMGCR, and SLC7A11) was identified based on the optimal value of λ (Fig. 4b, 4c). Survival analyses based on the optimal cut-off expression value of each gene showed that high expression of risk genes (CBS, CD44, AKR1C2, CHAC1, and SLC7A11) was correlated with poor prognosis, while high expression of protective genes (HMOX1 and HMGCR) displayed the opposite patterns (Additional file 3: Figure S2).
Based on the 7 candidate FRGs, the risk score of each patient was calculated according to the following formula: Risk score = (3.8463 × CBS) + (-0.0021 × HMOX1 ) + (0.0101 × CD44) + (0.0208 × AKR1C2) + (0.0252 × CHAC1) + (-0.1354 × HMGCR) + (0.3774 × SLC7A11). We then used the median risk score as a cutoff point for classifying KIRC patients in the training cohort (n = 294) into high-risk group (n = 147) and low-risk group (n = 147). Kaplan-Meier survival curve analysis showed that OS was significantly different between the predicted two risk groups and the high-risk group had significantly shorter survival time compared to low-risk group (P = 4.709e-06; Fig. 4d). The survival rates at three and five years in the high-risk group were 61.9% and 49.2%, respectively, while the corresponding rates in the low-risk group were 86.3% and 76.7%, respectively. The ROC curve analysis indicated the risk signature had a promising predictive value for KIRC survival prediction (AUC = 0.734) (Fig. 4e). We then ranked the risk scores of the patients and then analyzed their distributions (Fig. 4f). The distributions of risk scores and the survival status suggested survival rate and time were significantly increased in the low-risk group compared to the high-risk group (Fig. 4g). The heatmap displayed that the expression levels of the 7 candidate FRGs in the high- risk and low-risk groups (Fig. 4h).
Validation of the prognostic signature in the TCGA testing cohort, the TCGA entire cohort, and the ICGC cohort
To verify the accuracy and robustness of the prognostic risk signature, the prognostic value of the signature was further validated in the testing cohort (n = 192), the entire TCGA-KIRC cohort (n = 486). The prognostic risk score was calculated for patients in each cohort based on the prognostic signature. Patients of the testing cohort (n = 192) were stratified into the high-risk group (n = 96) and low-risk group (n = 96) with the median risk score as the cutoff. Similar to the procedure in the entire cohort, the 489 patients of the entire cohort were divided into the high-risk group (n = 243) and low-risk group (n = 243). The detailed clinical features of all KIRC patients were listed in Table 1.
We observed that the results in the testing and entire cohort were consistent with the outcome in the training cohort. Kaplan-Meier survival curves showed that patients in the high-risk group had a significantly worse OS than their low-risk counterparts in both the testing cohort (P = 2.149e-08) (Fig. 5a) and the entire TCGA-KIRC cohort (P = 7.724e-13) (Fig. 5d). In the testing cohort, the survival rate at 3 years in the high-risk and low-risk group was 65.1% and 87.9%, respectively. The survival rate of high-risk group was 32.4% at 5 years, and that of the low-risk group was 84.8%. Similarly, in the entire TCGA-KIRC cohort, the survival rates at 3 and 5 years in high-risk group were lower than those in the low-risk group (63.4 vs. 87.7% at 3 years, 41.3 vs. 80.2% at 5 years). ROC curve analysis showed that the AUC values for in the testing cohort and the entire TCGA-KIRC cohort were 0.762 (Fig. 5b) and 0.749 (Fig. 5e), respectively. The risk score distribution, survival status, and the risk gene expression in the testing cohort and the entire TCGA-KIRC cohort were shown in Fig. 5c and 5f.
Then the prognostic signature was validated in the ICGC cohort (n = 90). Patients were divided into high‐risk group (n = 45) and low‐risk group (n = 45) with the median risk score. The OS was significantly poorer in the high‐risk group than in the low‐risk group (P = 1.592e-02) (Additional file 4: Figure S3A). The AUC value for the prognostic signature was 0.71, suggesting well-prediction performances (Additional file 4: Figure S3B). The distribution of risk score, survival status, and gene expression of KIRC patients in the ICGC cohort was presented in Additional file 4: Figure S3C, which were similar to the above cohorts. Taken together, these results revealed that the risk signature could accurately predict the prognosis of KIPP patients.
The prognostic signature is an independent prognostic factor for KIRC patients
To determine whether the risk signature can be used as an independent prognostic factor, univariate and multivariate Cox regression analyses were performed with the risk score and relevant clinical factors, including age, gender, grade, AJCC stage, T stage, and M stage. In the TCGA training cohort, univariate analyses showed that age and risk score were significantly associated with OS. Besides, subsequent multivariate analyses suggested that age and risk score were still significantly associated with OS (Table 2). Similar results were obtained for both the testing and the entire TGCA-KIRC cohorts (Table 2). Therefore, the risk score calculated based on the prognostic risk signature was an independent adverse prognostic factor for OS in KIRC patients.
Prognostic risk score indicated strong associations with clinical characteristics in KIRC
To investigate whether the risk signature could better predict KIPC clinicopathological features, an analysis was performed to explore the associations between the risk signature and clinical parameters. Significant differences were observed between two groups in grade (P = 6.264e-06) (Fig. 6a), AJCC stage (P = 1.6973-06) (Fig. 6b), T stage (P = 4.884e-06) (Fig. 6c), and M stage (P = 0.002) (Fig. 6d). Simultaneously, we observed that advanced-stage tumors were significantly associated with the high‐risk group, whereas, the early-stage tumors were correlated with the low‐risk group.
To detect the prognostic value of our risk score model in different subgroups, we further investigated stratified survival analysis using the following clinical variables: age (≤60 and >60), gender (female and male), tumor grade (G1-2 and G3-4), AJCC stage (I & II and III & IV), T stage (T1-2 and T3-4), and M stage (M0 and M1). Interestingly, survival analysis indicated the high-risk group had a significantly lower OS rate than those in the low-risk group for all hierarchical cohorts (Fig. 7). Thus, these results suggested that the classification of the risk signature can be used to precisely identify patients with poor prognosis, without the consideration of clinical parameters.
Development of a personalized prognostic nomogram
A nomogram is a powerful tool that has been extensively applied to quantitatively determine individuals’ risk in clinical decision making by incorporating multiple clinical factors [36, 37]. To establish a viable method for predicting survival in patients with KIRC, we developed a prognostic nomogram based on the constructed prognostic risk signature and serval clinical features, including age, gender, grade, AJCC stage, T stage, M stage. The nomogram was devoted to evaluate the probability of 1-, 3-, and 5-year survival (Fig. 8a). Each factor was assigned a score in proportion to its risk contribution to survival in the nomogram. The C-index that was used to evaluate the OS of the nomogram was 0.771. Calibration curves showed optimal agreement when compared with an ideal model (Fig. 8b), particularly for 3- and 5-year survival predicted probabilities. Decision curve analysis (DCA) of the nomogram indicated that the nomogram had a wide and practical range of threshold probability for the TCGA-KIRC cohort for predicting survival rates (Fig. 8c).
The expression patterns, SNVs, CNVs of the seven candidate genes in the signature for KIRC patients
Subsequently, we analyzed the correlation between the expression of each of the seven prognostic risk signature genes and the clinicopathological features of KIRC patients in the TCGA cohort. In terms of grade alone, CBS, CD44, and CHAC1 increased with tumor grade, while HMGCR was decreased. No significantly difference in the expression values of HMOX1, AKR1C2, and SLC7A11 was detected between different tumor grades (Fig. 9a). As for different AJCC stage, CBS, HMGCR, CHAC1, and HMOX1 were significant differentially expressed, with higher expression levels of CBS and CHAC1 indicating higher advanced AJCC stage, while HMGCR and HMOX1 showed the opposite trend (Fig. 9b). Regarding T stage, it was noted that CBS, CD44, CHAC1, and SLC7A1 were significantly up-regulated in advanced T grade, whereas HMGCR was significantly down-regulated (Fig. 9c). CBS, CD44, and HMGCR also represented similar trends in N stage as T stage (Fig. 9d). In general, the expression levels of CBS, CD44, CHAC1, and SLC7A1 were positively correlated with tumor progression, while HMOX1 and HMGCR were negatively correlated with tumor progression, which was consistent with the above study. In addition, AKR1C2 expression appeared to be independent of tumor progression.
We then used the GSCA database to determine the single nucleotide variations (SNVs) of the seven candidate genes in the signature for KIRC patients. The results indicated that the most frequent mutation type was single nucleotide polymorphism (SNP) (Additional file 5: Figure S4A), and missense mutation was the most fraction among these mutations (Additional file 5: Figure S4B). In addition, C > T transversion accounted for the most common of SNV (Additional file 5: Figure S4C). The characteristic of the frequently mutated genes was showed in Additional file 4: Figure S4D. We also analyzed the copy number variations (CNVs) of the seven candidate genes in the signature and found heterozygous mutations (amplification and deletion) in all genes (Additional file 5: Figure S4E).
Drug sensitivity of the seven candidate genes in the signature for KIRC Patients
We next utilized the CellMiner database to explore the expression of seven candidate genes in diverse kidney cancer cell lines, including 786-O, A498, ACHN, CAKI-1, RXF 393, SN12C, UO-31, and TK-10. We observed that the expression levels of these genes showed great heterogeneity in different cell lines (Additional file 6: Figure S5A). In addition, we also investigated the drug sensitivity of the seven candidate genes for the KIRC patients using the GSCALite database. Four of them (CD44, SLC7A11, AKR1C2, and HMOX1) were highly related to drug sensitivity to a number of chemotherapy drugs (Additional file 6: Figure S5B), which provided direct support for drug targeted therapy.