RNA sequences and clinical data from 530 ccRCC patients and 72 normal kidney tissues were acquired from the Tumor Genome Atlas (https://portal.gdc.cancer.gov). Additionally, ninety-one RNA-SEQ datasets containing survival data for ccRCC patients were obtained from the International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org) for external validation. The GSE53757 and GSE15641 datasets were retrieved from the Gene Expression Overview database (https://www.ncbi.nlm.nih.gov/geo). The GSE53757 dataset (GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array) comprised 72 normal samples and 72 ccRCC samples, while the GSE15641 dataset (GPL96 [HG-U133A] Affymetrix Human Genome U133A array) included 23 normal and 32 ccRCC samples. Differential gene analysis was conducted using the GSE53757 dataset, while the performance of the diagnostic model and the expression of risk model genes were validated using the GSE15641 dataset. Subsequently, 282 immune checkpoint-related genes (ICRGs) were identified from the literature (Zhao et al., 2021).

Differentially expressed genes (DEGs) between ccRCC and normal groups in TCGA, as well as ccRCC and normal groups in GSE53757, were identified based on P - value < 0.05 and |log₂FC|>1 criteria utilizing the R language package “limma” (version 3.42.2) (Ritchie et al., 2015). Heatmaps were generated using the pheatmap package (version 1.0.12), while box plots were created using the ggplot2 package (version 3.3.2) (Ito and Murphy, 2013) In this study, the method section employs the Wilcoxon rank - sum test (P - value <0.05) to compare differences among different groups and evaluate the significance of differences in gene expression levels or other variables when creating boxplots.

The expression data of DICRGs were obtained from the ccRCC sample expression data in TCGA. This data was then integrated with clinical information to derive the clinical expression profiles of 526 ccRCC samples, excluding those with missing survival data. Subsequently, in order to evaluate the generalization ability of the model, the 526 samples were divided into a training cohort and a test cohort at a ratio of 6:4 (training cohort = 316, test cohort = 210). In the training cohort, a univariate Cox algorithm was employed to identify prognostically relevant DICRGs with significance level p < 0.2 (Ye et al., 2021). Following this, a multivariate Cox analysis was conducted using the step function. The Cox model formula utilized was:[ h(t/X) = h_0(t) exp (β_1 X_1 + β_2 X_2 + …… + β_p X_p) ] where (h_0(t)) represents the baseline risk function at time (t) when all variables are zero, (X_1, X_2, … , X_p) denote the influencing factor variables, and (β_1, β_2, … , β_p) refer to the regression coefficients. Patients were then scored based on the risk coefficients and expression values derived from the multivariate Cox analysis. The surv_cutpoint function within the survival package was utilized to determine the optimal cutoff value for the continuous independent variable of the survival profile (train:0.85; test:1.42; ICGC:1.18). Based on this cutoff, patients were categorized into high-risk and low-risk groups. Survival analyses were conducted using the Survival package (version 3.2-7) (N and Lee, 2019). Additionally, the SurvivalROC package (version 1.0.3) (Robin et al., 2011) was employed to calculate the area under the curve (AUC) values of the ROC curves, serving as a measure of the predictive model’s accuracy. Finally, the prognostic model underwent validation using both a test cohort and an external validation cohort.

In order to verify the expression level of prognosis-related DICRGs. First, the expression levels of the four prognosis-related DICRGs were validated using the Wilcoxon test method in the GSE15641 dataset. Then, the protein expression levels of DICRGs in both ccRCC and adjacent normal tissues were confirmed through immunohistochemical staining. Immunohistochemistry (IHC) results were acquired from the Human Protein Atlas database (HPA) available at https://www.proteinatlas.org/.

The correlation between clinicopathological factors and risk models was examined in the TCGA training cohort. Various variables such as age (≤65 or >65 years), sex (female or male), stage (stage I, II, III, or IV), pathologic T (T1, T2, T3, or T4), pathologic N (N0 or N1), pathologic M (M0 or M1), and grading (G2, G3, or G4) were used to categorize the training cohort into different groups. The findings of the correlation analysis were presented through box plots.

Independent prognostic analyses of clinicopathologic factors were conducted, and risk models were developed using univariate and multivariate Cox analyses. Subsequently, column line plots were created and visualized utilizing the rms R package (Qiu et al., 2022) based on 149 samples from the TCGA training cohort in ccRCC patients. Overall survival (OS) is often the primary observation index of choice in phase III clinical trials and has important clinical significance. Based on this, the index used in survival analysis in this study was OS. Furthermore, the model’s performance was evaluated using calibration curve analysis and ROC curves.

To better understand the relationship between risk model genes and risk scores, a correlation analysis was conducted using the Pearson, Spearman and Bayes methods, and scatter plots illustrating the correlation as well as histograms depicting the data distribution are created utilizing the R package “ggplot2 (3.2.1)” (Wickham, 2016). Subsequently, the expression levels of prognosis-related DICRGs were compared between the high and low-risk groups by Wilcox. test (P - value <0.05).

Using the TCGA and GSE53757 datasets, the expression levels of prognosis-related DICRGs in ccRCC and normal groups were determined using the Wilcox. test. Following the identification of prognostically relevant DICRGs, logistic regression (LR) and support vector machine (SVM) machine learning algorithms were employed to develop diagnostic models in both the TCGA and GSE15641 datasets. These models were subsequently validated in the GSE15641 dataset. Finally, the diagnostic models’ efficacy was assessed using ROC curves.

Within the training cohort, comprising 184 high-risk samples and 132 low-risk samples, risk-related DEGs were identified using the limma package (Ritchie et al., 2015) with criteria of P - value <0.05 and |log₂FC|>1. Volcano plots and heat maps were utilized to visualize the results. Functional enrichment analysis, including GO annotation and KEGG pathway analysis, was conducted using the R package “clusterProfiler” (Wu et al., 2021). The enrichment outcomes were presented visually using the ggplot2 package (Wickham, 2016).

The proportion of 22 immune cell types per ccRCC sample in both high-risk and low-risk groups was determined by estimating the relative abundance of RNA transcripts using the cell type identification algorithm Cibersorte (version 1.03) (Chen et al., 2018). Corresponding statistical values were computed, and samples with P - value <0.05 were selected for subsequent analysis. A heatmap illustrating the scores of the 22 immune cell types was generated based on their respective scores. Violin plots were created using the ggplot2 package. The immune phenotype score (IPS) of ccRCC patients was extracted from the TCGA database, and differences in IPS between the high-risk and low-risk groups were assessed. Additionally, correlation analysis was performed using Spearman’s correlation coefficient, and heatmaps were generated using the ggplot2 package (version 3.3.3).

Risk scores were calculated for the NCI60 cell line in the CellMiner database (60 cell lines), and patients with ccRCC were categorized into high-risk and low-risk groups using the median as the cut-off point. Correlation analysis of federal drug risk scores U.S. Food and Drug Administration (FDA)-approved drugs used in 60 cell lines, IC₅₀ was performed using Spearman, |cor|>0.4, P - value <0.05.

Differential analysis was conducted on extensive RNA sequencing data from TCGA (comprising 530 ccRCC and 72 control samples) and the GSE53757 dataset (with a Tumor: Normal ratio of 72:72). A total of 1,375 DEGs were identified between the ccRCC and Normal groups, comprising 624 upregulated genes and 751 downregulated genes in the ccRCC group (Figure 1A). Differential gene selection criteria included |log₂fold change|> 1 and P - value <0.05 Subsequently, 14 DICRGs were identified by intersecting DEGs and ICRGs. Among these 14 DICRGs, the expression of 13 genes was upregulated, while 1 gene was downregulated (Figure 1B). Visualization of the expression patterns of the 14 DICRGs in the TCGA and GSE53757 datasets was achieved through heatmaps and box plots (Figures 1C–F).

Four DICRGs (EGFR, TRIB3, ZAP70, CD4) prognostically relevant were identified through univariate and multivariate Cox analyses (P - value <0.2) (Table 1; Table 2). Among these, EGFR and CD4 exhibited a protective role in renal cell carcinoma, whereas TRIB3 and ZAP70 were identified as risk factors (P - value <0.05) (Figure 2A). Patients were stratified into high-risk and low-risk groups based on the optimal cutoff values (training cohort = 0.85; test cohort = 1.42; ICGC = 1.18) (Supplementary Tables S3-S5), and risk curves along with heatmaps were generated (Figure 2E). Survival analysis indicated a higher survival rate among the low-risk group in the training cohort (Figure 2H). The ROC curve analysis results for the training cohort demonstrated superior predictive performance of the model, with AUC values exceeding 0.643 at 1, 3, and 5 years (Figure 2B). Subsequently, we evaluated the predictive performance of the model using the test cohort and the external validation cohort. The findings were consistent with those of the training cohort, as evidenced by the risk profile, heatmap, and survival curves displayed for the test cohort (Figures 2F, I). The AUCs for 1-, 3-, and 5-year survival in the test cohort were 0.665, 0.642, and 0.651, respectively (Figure 2C). Similarly, in the external validation cohort, the results were comparable (Figures 2G, J). The AUC values for the 1-, 3-, and 5-year ROC curves in the external validation cohort were all above 0.599 (Figure 2D).

In the GSE15641 dataset, the expression levels of the four prognosis-related DICRGs exhibited significant differences between the ccRCC and normal groups, aligning with the expression trends observed in the TCGA dataset (Figure 3A). Within the HPA database, researchers conducted detailed examinations of each protein’s expression across 64 cell lines, 48 human normal tissues, and 20 tumor tissues, utilizing highly specific antibodies alongside immunodetection techniques such as immunoblotting, immunofluorescence, and immunohistochemistry. We conducted a search within the database to retrieve the immunohistochemical results for the four prognostic model genes in ccRCC tissues, presenting the expression patterns in both normal and ccRCC tissues(Figures 3B–E). The immunohistochemical findings revealed elevated expression levels of EGFR, TRIB3, ZAP70, and CD4 in ccRCC tissues, consistent with the observations from the TCGA database.

To further validate the expression changes of the differentially expressed immune checkpoint-related genes in clinical samples, and to enhance the reliability and accuracy of the research findings, we performed qPCR analysis of the expression differences of EGFR, TRIB3, CD4, and ZAP70 in 10 pairs of ccRCC tumor tissues and their adjacent normal tissues. The results showed that, compared to normal kidney tissues, EGFR, TRIB3, and CD4 were significantly upregulated in ccRCC tumor tissues (p < 0.05), while ZAP70 showed no significant difference (p = 0.9744), suggesting that it may not directly participate in the pathological process of ccRCC. Although these findings provide new insights into the molecular mechanisms of ccRCC, the small sample size (n = 10) and the limitations of qPCR technology require further functional experimental validation. The high expression of genes like EGFR may have a pro-cancer effect and serve as potential biomarkers for ccRCC molecular subtyping and targeted therapy, but larger-scale clinical validation is still needed. (Table 3; Figures 4A–D).

Significant (P - value <0.05) disparities in grading, pathologic T, pathologic M, and stage were detected by comparing the proportions of high-risk and low-risk patients across different subgroups (Supplementary Table S6). The Wilcox test revealed that the risk scores among different subgroups were notably significant (P - value <0.01) for stages G2-G4, T1-T3, and I-IV (Supplementary Figures S1A, B). Moreover, both univariate and multivariate Cox regression analyses demonstrated that the risk score independently influenced the prognosis of ccRCC patients (Supplementary Figure S1C). Column line plots, based on the four prognostically relevant DICRGs, were constructed to forecast patients’ overall survival at 1, 3, and 5 years (Figure 5A). The calibration curves’ slopes at 1, 3, and 5 years approximated 1, indicating a high compatibility between the predictions and actual outcomes (Figure 5B). Additionally, the area under the curves (AUCs) of 0.676, 0.666, and 0.669 were achieved at 1, 3, and 5 years, respectively, suggesting that the risk model exhibited significant prognostic value for ccRCC patients (Figure 5C).

A significant disparity in the expression levels of prognosis-related DICRGs was observed between the high-risk and low-risk groups. Specifically, the expression of TRIB3 and ZAP70 was markedly higher in the high-risk group compared to the low-risk group, whereas EGFR and CD4 exhibited significantly higher expression levels in the low-risk group than in the high-risk group (Figure 6A). Correlation analysis further revealed a positive correlation between TRIB3, ZAP70, and risk scores, while EGFR and CD4 showed a negative correlation (Figures 6B–E; Supplementary Figure S2). These findings indicate that TRIB3 and ZAP70 may have a protective role in renal cell carcinoma, whereas EGFR and CD4 could be considered as risk factors.

Following ANOVA analysis, all four prognostically relevant DICRGs within the TCGA and GSE53757 datasets exhibited significant upregulation in the ccRCC group (Figures 7A, B). The ROC curve results indicated that the AUC values of the four prognostically relevant DICRGs within the TCGA and GSE15641 datasets exceeded 0.85, demonstrating their robust diagnostic capability (Figures 7C, D). Subsequently, two machine learning algorithms, LR and SVM, were employed to develop the diagnostic models. The AUC values of these diagnostic models surpassed 0.9, indicating their effectiveness in accurately diagnosing ccRCC (Figure 7E). The validation outcomes within the GSE15641 dataset corroborated the performance of the diagnostic models, with the ROC curves demonstrating AUC values above 0.9 for both diagnostic models (Figure 7F).

A total of 20 risk-associated genes were identified in both the high-risk and low-risk groups, comprising 2 upregulated genes and 18 downregulated genes (Figures 8A, B). The differentially expressed genes (DEGs) were predominantly enriched in terms related to organic anion transport, anion transmembrane transport, and vascular processes within the circulatory system, as illustrated in Figure 8C and detailed in Supplementary Table S7. Moreover, the results of KEGG enrichment analysis revealed significant associations with the metabolism of ascorbic acid and glyoxylate, as well as the inter-conversion of pentose and glucuronide pathways, suggesting potential functional roles of the DEGs (Figure 8D; Supplementary Table S8).

The heatmap in Figure 9A displays the proportion of immune cells in each sample. Ten immune cell types exhibited significant differences (P - value <0.05) between the high-risk and low-risk groups, including CD8 T cells, resting memory CD4 T cells, T follicular helper cells, regulatory T cells (Tregs), NK cells, monocytes, M0 macrophages, M2 macrophages, resting mast cells, and neutrophils (Figure 9B). IPS was notably higher in the high-risk group compared to the low-risk group, as demonstrated in Figure 9C. Correlation analysis revealed that ZAP70 exhibited a positive correlation with CD8 T cells and a negative correlation with M2 macrophages, while CD4 showed a positive correlation with regulatory T cells (Tregs) and a negative correlation with resting memory CD4 T cells. Additionally, risk scores displayed a positive correlation with M0 macrophages and T follicular helper cells, and a negative correlation with M2 macrophages (Figure 9D).

The correlation analysis of cell-related risk scores and IC50 revealed significant associations with 24 drugs (arsenic trioxide, asparaginase, batrachotoxin, bendamustine, carmustine, loratadine, chlorambucil, cyclophosphamide, dasatinib, dexamethasone Decadron, dimethylaminobenzylamine, felitinib, fludarabine, fluphenazine, hydroxyurea, ifosfamide, imefamide, imethadone, inofosfene, nelarabine, oxaliplatin, pipecolonium bromide, PX-316, uracil mustard, and XK-469). These drugs exhibited significant correlations with the prognostic model ((|cor |> 0.4 and P - value <0.05), as illustrated in Figure 10.