The RNA-sequencing TPM data and corresponding clinical data of KIRC were retrieved from the TCGA database (https://portal.gdc.cancer.gov/), including 541 KIRC samples with complete survival data and 72 normal samples. The RNA-sequencing TPM data and corresponding clinical data of LGG were retrieved from the TCGA database (https://portal.gdc.cancer.gov/), including 479 LGG samples. The GSE22541 dataset was downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo) and included 68 CCRCC samples with survival data. The GSE53757 dataset was downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo) and included 72 CCRCC samples and 72 normal samples (von Roemeling et al., 2014). FDX1, LIPT1, LIAS, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A are thought to be cuproptosis-related genes and are involved in two structurally distinct copper-loaded ionophores (Tsvetkov et al., 2022).

RNA-sequencing TPM data from the TCGA database were used to compare the expression of cuproptosis-related genes in CCRCC specimens and normal specimens using the Wilcoxon rank-sum test. Statistical significance was set at p < 0.05. The GSE53757 dataset was used for further validation.

For significant prognosis-related cuproptosis-related genes, we performed gene co-expression analysis in TCGA CCRCC patients and set the absolute value of the correlation coefficient to greater than 0.4 with a p-value less than 0.001 to obtain the co-expression genes. To further understand the potential role of cuproptosis-related genes in CCRCC, GO and KEGG analyses were performed on co-expressed copper death-related genes.

We divided the training and testing cohorts into a ratio of 7:3 for patients with CCRCC. Clinical statistical analysis of the training and testing groups was performed using the chi-square test. In the training cohort, a univariate Cox regression analysis of cuproptosis-related genes was performed to identify the significant prognostically related genes. For significant prognosis-related cuproptosis-related genes, we used LASSO regression analysis to obtain independent prognostic genes in the training set. LASSO regression improves the accuracy and interpretability of the model and reduces the risk of overfitting (Tibshirani, 1997). Multivariate Cox regression analysis was conducted to obtain regression coefficients for independent prognostic genes. Finally, a 4-CRGs risk signature was established based on the multivariate Cox regression coefficient beta value, and the formula was as follows: risk score = EXPgene1∗ β1 + EXPgene2∗β2 + EXPgene3∗β3 + . . . + EXPgenen∗βn, where EXP is the expression level and β represents the regression coefficient from the multivariate Cox. In both cohort, by calculating the risk score for each sample, patients were divided into low- and high-risk groups using the median cut-off value. Furthermore, the KM curve was used to compare the overall survival (OS) between the two groups using the log-rank test. A time-dependent ROC curve analysis was used to assess the predictive power of the 4-CRGs risk signature. Finally, we perform external validation with the external validation cohort GSE22451 and TCGA-LGG. In each independent external validation cohort, based on the risk score, patients were classified into two groups. Maximally selected rank statistics was applied by using an R package “survival”, and “survminer” to identify the optimal cutting point to divide patients.

We screened for prognostic predictive factors including clinical characteristics and risk scores. Specifically, the univariate Cox proportional hazard model was employed to analyze the correlation between the risk score and OS, and multivariate Cox regression analysis was used to evaluate whether the established risk score could serve as an independent prognostic predictor. Further, to comprehensively assess patient survival, we constructed a nomogram integrating distinct clinicopathological information, including age, stage, and risk score, using the “rms” package. Additionally, the decision curve analysis (DCA) of 1, 3, and 5 years was calculated to evaluate whether the synthetic nomogram we established was suitable for clinical application.

We used the CIBERSORT algorithm to assess the degree of infiltration of 22 immune cells in different CCRCC samples (Newman et al., 2015). Single-sample gene set enrichment analysis (ssGSEA) was applied to explore the different infiltration degrees of immune-related functions in different CCRCC samples of the TCGA database using the R package “GSVA”.

TIMER is a website that can systematically analyze immune infiltration in various malignancies (https://timer.cistrome.org/) (Li et al., 2020). We investigated the relationship between gene expression and gene markers of TILs in CCRCC.

TISIDB (http://cis.hku.hk/TISIDB/index.php) was used to investigate the association of genes with immunostimulators in CCRCC.

The “pRRophetic” R package was used to predict the half-maximal inhibitory concentration (IC50) of some drugs in each sample regarding tumor treatment.

We used RNAseq data from the TCGA database and miRNAseq data, including 541 KIRC samples and 72 normal samples. The difference analysis was performed using the DESeq2 package, and |logFC|>1 and adj. p < 0.05 were set as thresholds to obtain the differential expression of lncRNAs (DElncRNAs) and miRNAs (DEmiRNAs) between CCRCC patients and normal patients. Subsequently, cuproptose-related miRNAs (CRMs) were predicted using the miRTarBase database (Huang et al., 2020). DEmiRNAs and CRMs were intersected to obtain cuproptosis-related DEmiRNAs (CRDEMs). Cuproptosis-related lncRNAs (CRLs) were predicted using the starBase database (Li et al., 2014). CRLs and DElncRNAs were intersected to obtain the cuproptosis-related DElncRNAs (CRDELs). Subsequently, we integrated the interactions between CRDEMs, CRDELs, and cuproptosis-related genes to construct a ceRNA regulatory network. Finally, Cytoscape (version 3.8.0) software was used to visualize the ceRNA regulatory network.

All statistical analyses were conducted using R version 4.1.2 software (https://www.r-project.org/). Univariate Cox hazard regression analyses were performed to identify the independent prognostic cuproptosis-related genes. Survival analysis was conducted by the Kaplan-Meier (K-M) method with the log-rank test. We also compared the expression of cuproptosis-related genes at different clinical stages by using the Wilcoxon rank-sum test.

We compared the expression of cuproptosis-related genes between 541 CCRCC samples and 72 normal samples using the Wilcoxon rank-sum test in the TCGA cohort. We found that LIAS and CDKN2A were significantly upregulated in tumor samples, and FDX1, DLD, DLAT, PDHA1, PDHB, MTF1, and GLS were significantly downregulated in tumor samples (Figure 1A).

To obtain reliable survival results for CCRCC, we first excluded samples with a survival time of less than 30 days. In total, 518 samples were obtained (Table 1). Nine differentially expressed cuproptosis-related genes (FDX1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, LIAS, and CDKN2A) were identified. Through the KM curve, we found that 9 selected genes all had an impact on the prognosis of CCRCC, including FDX1 (hazard ratio, HR = 0.49; 95% confidence interval, 95%CI = 0.36–0.67; p < 0.001), LIAS (HR = 0.61; 95%CI = 0.45–0.83; p = 0.002), DLD (HR = 0.53; 95%CI = 0.39–0.72; p < 0.001), DLAT (HR = 0.46; 95%CI = 0.33–0.63; p < 0.001), PDHA1 (HR = 0.63; 95%CI = 0.46–0.85; p = 0.003), PDHB (HR = 0.57; 95%CI = 0.42–0.78; p < 0.001), MTF1 (HR = 0.59; 95%CI = 0.43–0.80; p = 0.001), GLS (HR = 0.65; 95%CI = 0.48–0.88; p = 0.005), and CDKN2A (HR = 1.49; 95%CI = 1.10–2.02; p = 0.011) (Figure 1B).

We compared the expression of cuproptosis-related genes at different clinical stages using the Wilcoxon rank-sum test. FDX1, LIAS, DLD, DLAT, PDHA1, PDHB, MTF1, and GLS were highly expressed in T1 compared with T3, and CDKN2A was expressed at lower levels in T1 than in T3 (Supplementary Figure S1). Compared to N1, FDX1 and LIAS were highly expressed in N0 (Supplementary Figure S2). FDX1, LIAS, DLD, DLAT, PDHA1, PDHB, MTF1, and GLS were highly expressed in M0 compared with M1, and CDKN2A was expressed at lower levels in M0 than in M1 (Supplementary Figure S3). FDX1, LIAS, DLD, DLAT, PDHA1, PDHB, and MTF1 were highly expressed in Stage 1 compared to Stage 3 and 4, and CDKN2A had lower expression in Stage 1 compared to Stage 3 and 4 (Supplementary Figure S4).

For significant prognosis-related cuproptosis-related genes, we performed gene co-expression analysis in TCGA tumor patients and set the absolute value of the correlation coefficient to greater than 0.4 with a p-value less than 0.001 to obtain the co-expression genes. For co-expression genes, we performed GO and KEGG enrichment analyses and sorted them by p values (Figure 2). We found that co-expression genes were significantly enriched in the mitochondria during cell localization. The TCA cycle is thought to be associated with cancer progression, the site of biological processes in the mitochondria. This is consistent with Tsvetkov et al. (2022)’s view that Cu causes cell death by influencing the TCA cycle. The enrichment analysis results showed that co-expression of genes is correlated with autophagy and ubiquitin-mediated proteolysis, which provides a research direction for further exploration of the mechanism of cuproptosis.

We found no significant differences between the training and testing cohorts in terms of age, sex, or TNM staging (Table 1). In the training cohort, univariate Cox regression analysis yielded eight cuproptosis-related genes that were significantly associated with prognosis (Figure 3A). Using lasso regression method, six optimal variables were obtained from the above 8 cuproptosis-prognostic-related gene (Figures 3B,C). By Cox regression analysis, the signature was finally established: risk score = EXP FDX1 ∗ −0.499501220246694 + EXP DLD ∗ −0.59322127824406 + EXP DLAT ∗ −0.659153532219121 + EXP CDKN2A ∗ 0.199116740963518. The KM curve showed that the prognosis of the high-risk group was worse than that of the low-risk group (Figure 3D, log-rank p < 0.001; HR = 2.55, 95%CI = 1.73–3.76). ROC curves were used to assess the accuracy of the established models in predicting overall survival (OS) in patients with CCRCC. As shown in Figure 3H, the AUC values at 1, 3, and 5 years were 0.684, 0.688, and 0.670, respectively, indicating the robustness and accuracy of the model in predicting patient prognosis.

In the testing cohort, the KM curve showed that the prognosis of the high-risk group was worse than that of the low-risk group (Figure 4A, log-rank p = 0.006; HR = 2.31, 95%CI = 1.28–4.17). ROC curves were used to assess the accuracy of the established models in predicting OS in patients with CCRCC. As shown in Figure 4E, the AUC values at 1, 3, and 5 years were 0.665, 0.632, and 0.666, respectively, indicating the robustness and accuracy of the model in predicting patient prognosis.

We validated the differences in the expression of FDX1, DLD, DLAT, and CDKN2A between CCRCC and normal samples using GSE53757. The results showed that FDX1, DLD, and DLAT exhibited low expression in CCRCC, whereas CDKN2A was highly expressed in CCRCC (Figures 4F–I). This is consistent with the results obtained from the TCGA dataset. We further verify the above prediction method in external data cohorts” GSE22541” and “TCGA-LGG”. In GSE22541 validation cohort, we divided the CCRCC patients into high-risk and low-risk groups based on the risk score. Survival comparison showed that low-risk group had significantly better prognosis outcome than high-risk group (Figure 4J). In addition, In TCGA-LGG validation cohort, we also found that based on the high and low-risk groups divided by risk score. Different groups have significantly different prognostic outcomes (Figure 4K). This demonstrates the generalization power of cuproptosis-related signature and has some value for the prediction of other cancers.

Univariate and multivariate Cox regression analyses showed that risk score, stage, and age were prognostic predictors of TGGA-KIRC (Figures 5A,B). Nomograms are widely used for the prognostic assessment of tumors. Various clinical features have prognostic value in clinical practice. Therefore, we established a nomogram containing multiple clinicopathological characteristics and risk scores. The scores for each variable were calculated and combined to predict the prognosis of patients with CCRCC (Figure 5C). DCA (Figure 5E) also proved that the nomogram combined with various clinical features had a better clinical application value.

We used the CIBERSORT algorithm to estimate differences between 22 tumor-infiltrating immune cells between the low- and high-risk groups. Figure 6D shows that plasma cells, T cells CD8, T cells, regulatory T cells (Tregs), and activated NK cells were more enriched in high-risk groups, while naive B cells, T cells, CD4 memory monocytes, macrophages M0, macrophages M1, and macrophages M2 were more enriched in the low-risk group. This indicates that there are differences in immune cell infiltration in different risk groups, suggesting that cuproptosis -related genes are closely related to immune cell infiltration.

The ssGSEA method was applied to KIRC patients in the high- and low-risk groups to assess the differences in immune function between the high- and low-risk groups. Figure 7A shows that Type-I-IFN-Response, HLA, Cytolytic activity, Inflammation-promoting, T-cell-co-inhibition, Checkpoint, T-cell-co-stimulation, APC-co-stimulation, CCR, and parainflammation were upregulated in the high-risk group, suggesting that cuproptosis-related genes are involved in immune regulation.

We further investigated the links between FDX1, DLD, DLAT, and CDKN2A and the TIL gene markers in the TIMER database (Supplementary Material S1). DLD was strongly correlated with STAT3 (rho = 0.419), STAT5B (rho = 0.47), and CD4 (rho = −0.401) (Supplementary Figure S5). There was a strong correlation between DLAT and TGFBR2 (rho = 0.49), STAT3 (rho = 0.501), AHR (rho = 0.449), STAT5B (rho = 0.571), MRC1 (rho = 0.437), CD7 (rho = −0.417), and TGFBR2 (rho = 0.49) (Supplementary Figure S5). STAT3, STAT5B, and CD7 simultaneously show a strong correlation with DLD and DLAT, suggesting that STAT3, STAT5B, and CD7 may be associated with cuproptosis genes in an important link with levels of immune infiltration.

We explored the relationship between cuproptosis-related genes and immunomodulators associated with model building and found that FDX1, DLD, DLAT, CDKN2A, and immunostimulants were significantly associated. Specifically, high expression of FDX1 was significantly associated with TNFRSF8 (Rho = −0.402), and high expression of DLD was significantly correlated with TNFRSF4 (Rho = −0.648), TNFRSF18 (Rho = −0.435), and TNFRSF25 (Rho = −0.6). DLAT was significantly correlated with TNFRSF4 (Rho = −0.522), TNFRSF8 (Rho = −0.436), TNFRSF18 (Rho = −0.557), TNFRSF25 (Rho = −0.642), CDKN2A, and TNFRSF18 (Rho = 0.398) (Supplementary Figure S6). Hence, FDX1, DLD, DLAT, and CDKN2A may play important roles in immune interactions and may be associated with tumor immune evasion.

We found that PDCD1 and CTLA4 (immune checkpoints) were elevated in patients with high-risk scores (Figures 6B,C). These discoveries suggested that patients of high-risk scores may be more sensitive to ICB therapy.

Finally, we predicted some drugs for tumor treatment using the pRRophetic R package and obtained some drugs that may show different sensitivities in patients in the high-risk and low-risk groups. Specifically, the low-risk group was more sensitive to AKT. inhibitor.VIII, AP.24534, AS601245, AUY922, axitinib, and AZ628, and the high-risk group were more sensitive to A.443654, ABT.888, AG.014699, AICAR, and AMG.706 (Supplementary Figure S7).

We analyzed DElncRNAs and DEmiRNAs between 541 KIRC samples and 72 normal samples and obtained a total of 255 DEmiRNAs, of which 129 were upregulated, 126 were downregulated; and 1010 DElncRNAs, of these, 777 were upregulated and 233 were downregulated. Using the miRTarBase database, 148 CRMs were identified. Thirteen CRDEMs were obtained by intersecting the CRMs with the DEmiRNAs. Subsequently, 795 CRLs were predicted using the starBase database, and 31 CRDELs were obtained by intersecting CRLs with DElncRNAs. Based on the existing theory that lncRNA inhibits the degradation of mRNA by miRNA through competitive binding of miRNA, we constructed a ceRNA network containing one mRNA, two miRNAs, and 12 lncRNAs (Figure 7E).