Prognostic significance of key immune cell functional alterations in clear cell renal cell carcinoma

Introduction

Renal cell carcinoma (RCC), ranking as the eighth most prevalent malignancy globally, is projected to account for an estimated 81,610 new cases and 14,390 cancer-related fatalities in the United States during the year 2024 (1). Clear cell renal cell carcinoma (ccRCC), also known as kidney renal clear cell carcinoma (KIRC), is the most common histological type of RCC. The constellation of risk factors associated with the pathogenesis of ccRCC encompasses cigarette smoking, obesity, and hypertension. In parallel, genetic mutations, chromosomal aberrations, transcriptional dysregulation, and epigenetic alterations are recognized as pivotal contributors to the multifaceted process of ccRCC.

ccRCC is classified as an immunogenic tumor (2) and is characterized by abundant leukocyte infiltration, which includes CD8⁺ T cells, CD4⁺ T cells, natural killer (NK) cells, and myeloid cells (3,4). Among the tumor-infiltrated lymphocytes (TILs) in ccRCC, T lymphocyte is the predominant immune cell type within the tumor microenvironment (TME), with an average prevalence of 51% (5). The proliferative activity of CD8⁺ T cells has been reported to be a better predictor of longer patient survival (6). NK cells are cytotoxic lymphocytes that are able to kill tumor cells in the absence of antigenic priming or prior activation. NK cells comprise approximately 25% of lymphocytes in healthy renal tissue. And in ccRCC, an elevated proportion of NK cells among TILs has been associated with improved patient prognosis (4,7,8). Moreover, myeloid cells are a heterogeneous population of bone marrow-derived cells that contribute to tumor progression and metastasis by promoting angiogenesis and suppressing anti-tumor immune responses. In ccRCC, tolerogenic dendritic cells within the TME have been correlated with an unfavorable prognosis, significantly contributing to the immunosuppressed state of the ccRCC TME (9). Cytokines serve as crucial immunomodulatory molecules that orchestrate cellular communication within the immune system and the TME. Yet, to the best of our knowledge, the precise functional alterations of specific immune cell types and their intricate crosstalk with other cellular components within the context of ccRCC remain to be fully elucidated.

With more cases detected incidentally during routine health examinations at an early stage and the advent of novel therapies, the mortality rate of ccRCC was declining over the past decades. And the 5-year overall survival (OS) probability increased from 50% in the 1970s to 78% in the 2010s. However, it is noteworthy that approximately 17–30% ccRCC cases are diagnosed at an advanced stage, which is often associated with more aggressive disease progression and poorer treatment outcomes. Worse still, nearly 30% of local ccRCC cases developed into advanced disease post-operation with a poor prognosis. Up to now, there is a lack of validated molecular biomarkers for prognostication in the clinical practice of ccRCC. Thus, the development of new effective molecular signatures to predict the prognosis and guide further individual precise treatment of ccRCC patients is necessary. Despite an increasing number of studies (10,11) having explored the potential prognostic markers of ccRCC with promising predictive ability, most of them were based on bulk RNA sequencing (RNA-seq), which mainly focused on the “average” expression of all cells in the tumor samples and could not mirror the exact cellular and molecular changes. As a new advanced sequencing technology, single cell RNA sequencing (scRNA-seq) is capable of quantifying gene expression at the single-cell level and revealing the heterogeneity of gene expression in individual cells. This technology holds significant potential to enhance our comprehension of the immune cells infiltrating the TME of ccRCC and to discover effective prognostic markers correlated with immune cell dynamics.

In this study, we aimed to identify the specific immune cell types infiltrated in ccRCC by integrating scRNA-seq and bulk RNA-seq data. Then polarization of the key immune cells, enrichment of corresponding cytokines and cellular communication alterations in ccRCC samples based on scRNA-seq were explored. Finally, this study sought to develop a prognostic model utilizing genes associated with key immune cells to predict the prognosis and potential effective anticancer drugs for ccRCC patients. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-971/rc).

Methods

Bulk RNA-seq data source and preprocessing

The whole genomic expression data, clinical data and single nucleotide variation data of 542 KIRC and 72 normal samples were acquired from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). With the samples of 31 patients excluded due to a short follow-up time, a total of 511 tumor samples with clinical survival data were obtained. Gene expression profiles, accessioned as Gene Expression Omnibus Series (GSE)105261, comprising a cohort of nine ccRCC tumor samples and nine corresponding normal renal tissue samples, were extracted from the Gene Expression Omnibus (GEO) database. Additionally, Microarray expression profiles of 101 ccRCC tumor samples was downloaded from the E-MTAB-1980 cohort through ArrayExpress official website (https://www.ebi.ac.uk/biostudies/arrayexpress), with corresponding clinical data obtained from the literature published by Sato and colleagues (12). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The workflow of the whole study is shown in Figure 1.

Figure 1 The workflow of the present study. ccRCC, clear cell renal cell carcinoma; GSE, Gene Expression Omnibus Series; QC, quality control; ROC, receiver operating characteristic; TCGA-KIRC, The Cancer Genome Atlas-Kidney Renal Clear Cell Carcinoma; TMB, tumor mutation burden.

Single-cell data source and preprocessing

The scRNA-seq dataset GSE210038, which encompasses seven tumor tissue samples of ccRCC and two control samples of para-cancerous tissue, was retrieved from the GEO database. The single-cell raw data in GSE210038 were processed using the “Seurat” R package (version 4.2.0) in R software. Firstly, quality control criteria were as follows: (I) cells with genes between 200 and 4,000 were retained; (II) cells with unique molecular identifier (UMI) reads below 20, 000 were retained; (III) cells with less than 20% of mitochondrial genes were preserved, and (IV) cells with gene number per read counts more than 0.8 were preserved. After the data was normalized, the high variability genes of single cells were identified by balancing the relationship between average expression and dispersion. Then, principal component analysis (PCA) was performed and a significant PC was used as an input for the graph-based clustering. The batch effect between different samples was effectively mitigated in this study utilizing the “harmony” computational method. For clustering, we employed the “FindClusters” function, which effectively delineated 23 distinct clusters on the basis of 30 PCs, utilizing a resolution parameter set at 0.6. This unsupervised clustering approach was grounded on the “shared nearest neighbor (SNN)” modularity-optimized clustering algorithm. We then applied the “RunUMAP” function to perform Unified Manifold Approximation and Projection (UMAP). Cell aggregation was proved using UMAP-1 and UMAP-2. To identify differentially expressed genes (DEGs) in each cell cluster, we applied the “FindAllMarkers” function to the normalized gene expression data with parameters selected from the default parameters set by “Seurat”. Subsequently, the cell clusters were identified by cell type-specific biomarkers, and the proportion of cell types was calculated and assessed.

Deconvolution of the scRNA-seq data

Based on the gene expression data of various types of immune cells in the scRNA-seq of ccRCC, the relative enrichment score of each cell type was quantified according to the gene expression profile of each sample in the dataset integrating TCGA-KIRC and GSE105261. Batch effects caused by non-biotechnological bias were corrected using the “ComBat” function of R package “sva”. PCA was applied to check the degree of correction. In addition, we performed correlation analysis among different immune cell types.

Selection of key immune cells

Support vector machine-recursive feature elimination (SVM-RFE) was employed to train subsets of features from different categories to narrow down the feature set and find the most predictive features. In order to compute and retain a linear regression model that preserves the most informative variables, we employed the “glmnet” package (version 4.1.4) for least absolute shrinkage and selection operator (LASSO) regression. The binomial distribution variables were then used for LASSO classification and the model with the minimum cross-validated error was selected for its optimal balance between predictive performance and model complexity, which was ascertained through a rigorous 10-fold cross-validation process. The “RandomForest” function was deployed to perform random forest (RF) analysis. The minimum error was chosen as the mtry node value, and the image value that tends to be stable was chosen as the ntree. According to the mean reduction precision (MDA) and mean reduction Gini rate (MDG) of the feature weights, the top three immune cell types were selected as the key immune cell types for RF screening. Then, by integrating the analytical outcomes derived from LASSO regression, SVM-RFE and RF analysis, the most significant characteristic immune cell types were selected as the key immune cells in this study.

Immune response enrichment analysis

Based on the DEGs of key immune cells between ccRCC group and healthy controls, we employed Immune Response Enrichment Analysis (IREA) (13) to delineate the predominant cytokines to which key immune cells respond with heightened sensitivity. This analytical approach also enabled us to ascertain the polarization states of the key immune cells within the ccRCC microenvironment.

Cellular communication analysis and ligand-receptor expression

CellChat object is created based on UMI of each group (ccRCC group and healthy control group) by R package “CellChat” (https://www.github.com/sqjin/CellChat). The “CellChatDB.human” ligand-receptor interaction database was used as the reference database for intercellular communication analysis with default parameters. Comparisons of the total number and strength of interactions were obtained by combining the CellChat objects of each group with the “mergeCellChat” function. The “netVisual_diffInteraction” function was employed to visualize the difference in the number or strength of interactions between different cell types in each group.

Construction and validation of prognostic model

The ‘Seurat’ R package was utilized to rigorously identify the DEGs that distinguish the key immune cells between the ccRCC and healthy control groups [|log₂fold change (FC)| >0.25, adjusted P<0.05]. The R package “limma” (version 3.50.0) was applied to identify DEGs between ccRCC samples and healthy controls [|log₂(FC)| >0.25, adjusted P<0.05]. Overlapped genes of the two groups of DEGs were obtained. To evaluate the prognostic value of the overlapped genes, tumor samples from TCGA-KIRC and E-MTAB-1980 were utilized as training set (n=511) and independent validation set (n=101), respectively. Univariate Cox regression analysis was conducted to assess the correlation between the expression of each gene and OS within the training set. Genes that demonstrated a statistically significant association with OS, with P<0.05, were deemed to be significantly correlated with survival outcomes and were subsequently selected for further in-depth analysis. The LASSO Cox regression model in R package “glmnet” was then employed to narrow down the candidate genes and constructed a prognostic model. The penalty parameter (λ) was determined by the minimum criterion. The risk score (RS) is calculated using the following formula:

Riskscore=∑i=1nCoef(genei)∗Expression(genei)[1]

In the formula, “Coef(genei)” represents the risk coefficient of the selected genes in the prognostic model and “Expression(genei)” refers to their expression value. “surv_cutpoint” function was utilized to select the best truncation value to divide the training set samples into low-risk group and high-risk group. The Kaplan-Meier method was employed to plot survival curves to assess prognostic prediction value of the model, and log-rank test was performed to ascertain the significance of differences in survival curves between groups. The efficacy of the prognostic model was substantiated utilizing the receiver operating characteristic (ROC) curve analysis. The validation cohort was stratified into low- and high-risk groups based on the prognostic model to further ascertain its predictive accuracy.

Construction and validation of the nomogram

Clinical information of patients was extracted from the TCGA cohort. Univariate and multivariate Cox regression analyses were conducted to integrate clinical covariates with the RS derived from the prognostic model. Subsequently, we integrated RS with pertinent clinical parameters to construct a predictive nomogram, utilizing the “rms” R package (version 6.3.0). The predictive accuracy of the nomogram was rigorously assessed through calibration curves.

Analysis of somatic mutation and drug sensitivity

A landscape of genomic variation was depicted based on mutation data from TCGA-KIRC datasets. The R package “maftools” (version 2.10.05) was employed to show somatic variation, including single nucleotide polymorphisms, insertions and deletions, tumor mutation burden (TMB), and mutation frequency between different groups. In general, frequently mutated genes with the top 20 mutation frequencies are considered to be the main driver genes of malignancies.

Based on the maximum half inhibitory concentration (IC₅₀) data and corresponding gene expression data downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/), the R package “oncoPredict” (version 0.2) was employed to predict potential therapeutic drug sensitivity in patients stratified as high-risk and low-risk for ccRCC.

In order to explore the potential application of existing ccRCC drugs, the prognostic genes were submitted to The Drug Gene Interaction Database (DGIdb: https://www.dgidb.org) as potential targets, to search for the existing agonist or antagonists.

Validation of expression patterns of prognostic genes via the human protein atlas

The protein expression of prognostic genes in ccRCC and normal kidney tissues was determined using immunohistochemistry from The Human Protein Atlas (HPA) (https://www.proteinatlas.org/), which is a valuable database providing extensive transcriptome and proteomic data for specific human tissues and cells.

Statistical analysis

R software (version 4.1.2) was utilized for statistical analysis in this study. Spearman correlation test was employed to ascertain the correlation between two parameters. The Wilcoxon test was utilized to compare the differences between the two groups. Bilateral P<0.05 were considered statistically significant.

Results

Analysis of scRNA-seq data

We utilized the scRNA-seq dataset GSE210038 to perform a comprehensive analysis of cellular clusters in ccRCC. After quality control, visualization of UMAP showed 55,973 cells from the ccRCC samples and normal samples were clustered into 23 subgroups (Figure 2A). Annotating cell types based on the gene expression characteristics and marker genes (14) (available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-1.xlsx), these identified clusters were defined as 13 major cell types, including CD4⁺ T cells, CD8⁺ T cells, B cells, NK cells, and mast cells, etc. (Figure 2B). The dot plot was employed to show the expression levels of marker genes in each cell type (Figure 2C). Finally, the distributional profiles of disparate cell types between the ccRCC group and the paracancerous control group was shown (Figure 2D). We observed that following tumorigenesis, there was a significant increase in the population of CD4⁺ T cells and CD8⁺ T cells, while the proportion of renal tubules and endothelium, among other cellular components, notably decreased in ccRCC.

Figure 2 Identification of cell subgroups by scRNA-seq analysis. (A) UMAP plot showing the distribution of cell subsets in the ccRCC group and the paracancer control group. (B) UMAP plot showing annotated results of the cell subgroups in the ccRCC group and the paracancer control group. (C) Dot plot showing the expression of the marker genes of each cell subgroup. (D) Cumulative histograms showing the distribution of cell types in the ccRCC group and the paracancer control group. CAF, cancer-associated fibroblast; ccRCC, clear cell renal cell carcinoma; NK, natural killer; PST, proximal straight tubule; scRNA-seq, single-cell RNA sequencing; UMAP, Unified Manifold Approximation and Projection.

Analysis of immune cells infiltrated in ccRCC

The immunophenotypic landscape characterized by the gene expression profiles of diverse immune cell types, as delineated in the scRNA-seq datasets, was subsequently deconvoluted and projected onto the bulk RNA-seq dataset. Consequently, the levels of immune cell infiltration within the ccRCC cohort, encompassing B lymphocytes, CD4⁺ T lymphocytes, CD8⁺ T lymphocytes, macrophages, mast cells, monocytes, and NK cells, were markedly elevated compared to those observed in the control group (P<0.001, Figure 3A). A positive correlation was observed among the various immune cell populations infiltrating the TME, e.g., CD4⁺ T cells exhibited a significant positive correlation with other immune cell types, including CD8⁺ T cells, B lymphocytes, NK cells, monocytes and macrophages (Figure 3B).

Figure 3 Infiltration levels of immune cells between ccRCC group and control group. (A) Box plot showing the estimated proportion of infiltrated immune cell types between ccRCC group and control group. (B) Correlation between different infiltrated immune cell types. **, P<0.01; ****, P<0.0001. ccRCC, clear cell renal cell carcinoma; NK, natural killer.

Screening of key immune cells with machine learning algorithms

Key immune cell types associated with ccRCC were further screened with SVM-RFE, LASSO regression and RF analyses. Utilizing LASSO regression analysis, we successfully identified seven pivotal immune cell types (Figure 4A,4B). By employing the RF algorithm, we identified the intersection of the first three immune cell types based on their MDA and MDG indices, which were subsequently selected as the key immune cell types (Figure 4C-4E). Then, seven key immune cell types were selected by SVM-RFE (Figure 4F,4G). Finally, we performed an intersection analysis of the key immune cell types identified by each computational method. This integrative approach facilitated the identification of the two most critical immune cell types, CD8⁺ T cells and NK cells, which we called the key immune cells subsequently (Figure 4H). The diagnostic efficacy of the key immune cells was rigorously assessed utilizing ROC curve analysis. The analysis revealed that the area under curve (AUC) values for the key immune cells all exceeded 0.8, a threshold that denotes excellent diagnostic performance, indicating CD8⁺ T cells and NK cells hold promise as reliable diagnostic markers to distinguishing between diseased and non-diseased states with high accuracy (Figure 4I,4J).

Figure 4 Screening of key immune cell types associated with ccRCC utilizing machine learning algorithms. (A) LASSO coefficient profiles of seven immune cell types infiltrated in ccRCC, immune cell types are represented by different colors. (B) LASSO regression with 10-foldCV, and selection of the optimal parameter (Lambda) in the LASSO model. (C) Line chart showing effect of the number of decision trees on the error rate in RF analysis. (D,E) Dot plot showing mean decrease accuracy and mean decrease Gini of potential targets, the top three immune cell types were selected as the key immune cell types for RF screening. (F,G) SVM models with the highest accuracy and lowest error rate were established on seven immune cell types. (H) Venn diagram showing the intersection of three machine learning algorithms. (I,J) ROC curve of CD8⁺ T cells and cells in predicting ccRCC. AUC, area under the curve; ccRCC, clear cell renal cell carcinoma; CV, cross validation; LASSO, least absolute shrinkage and selection operator; NK, natural killer; OOB, out-of-bag error; RF, random forest; ROC, receiver operating characteristic; SVM-RFE, support vector machine-recursive feature elimination.

IREA of the key immune cells

To elucidate the immune polarization phenotypes of the key immune cells in ccRCC, DEGs of the key immune cells between ccRCC group and healthy control group were explored. The heat map showed in detail the expression of ten characteristic genes of the key immune cells with significant differences between the two groups (Figure 5A, available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-2.xlsx). Subsequently, IREA was performed based on the significantly upregulated genes in ccRCC. The IREA radar plot delineated that CD8⁺ T cells predominantly exhibited a T8-c polarization state (Figure 5B, available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-3.xlsx), while NK cells displayed an NK-f polarization phenotype (Figure 5C, available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-4.xlsx) in the ccRCC group. The IREA cytokine enrichment plot showed that CD8⁺ T cells were found to mainly enrich cytokines such as interferon (IFN)-α1 and interleukin (IL)-36α (Figure 5D, available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-5.xlsx). Concurrently, the IREA cytokine enrichment plot demonstrated that NK cells were primarily associated with the enrichment of IL-18 and IL-2 within the ccRCC cohort (Figure 5E, available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-6.xlsx).

Figure 5 IREA of key immune cell types. (A) Heatmap showing ten significant DEGs of key immune cell types between ccRCC group and healthy control group. (B) IREA radar plot showing the polarization state of CD8⁺ T cells in ccRCC. (C) IREA radar plot showing the polarization state of NK cells in ccRCC. (D) IREA cytokine enrichment plot showing cytokine enrichment of CD8⁺ T cells in ccRCC. (E) IREA cytokine enrichment plot showing cytokine enrichment of NK cells in ccRCC. ccRCC, clear cell renal cell carcinoma; DEGs, differentially expressed genes; IREA, Immune Response Enrichment Analysis; NK, natural killer.

Analysis of cellular communication

To further investigate the interaction networks of cells in ccRCC and control groups, as well as to investigate how CD8⁺ T cells and NK cells modulate intercellular communication in influencing tumorigenesis, we applied the R package “Cellchat” to reveal changes of cell crosstalk in the two groups, particularly focusing on the signaling events involving CD8⁺ T cells or NK cells as key communicators. Compared with the control group, an escalated total number and strength of interactions were observed in the ccRCC group (Figure 6A-6D). Subsequently, a detailed comparative analysis of the signaling patterns between the ccRCC group and the control group was undertaken. The overall signaling patterns derived from the ccRCC and control groups are distinctly illustrated in Figure 6E. For instance, the strength of the ANNEXIN signaling cascade, originating from NK cells, and the protease activated receptors (PARs) signaling cascade, emanating from CD8⁺ T cells, exhibit significant alterations in the ccRCC group. The distinct patterns of incoming signals within the ccRCC group and the control group are explicitly depicted in Figure S1. For instance, the strength of the pleiotrophin signal, which exerts its effects on CD8⁺ T cells, is observed to be augmented in the ccRCC group. The patterns of outgoing signaling emanating from the ccRCC group and the control group are clearly delineated in Figure S2. For instance, the strength of the ANNEXIN signal, originating from NK cells, is observed to be attenuated in the ccRCC group. Conversely, the signal strength of PARs, derived from CD8⁺ T cells, is augmented in the ccRCC group. Subsequently, our investigative focus narrowed onto the alterations in cellular communications within the key immune cells. An in-depth analysis was conducted to identify the receptor-ligand pairs that potentially mediate the intercellular communications between the key immune cells and other cellular components within the TME. The PARs signaling pattern has been identified as a critical immunomodulatory axis in ccRCC. This is evidenced by an augmented interaction strength between Granzyme A (GZMA) ligands, which are secreted by NK cells and CD8⁺ T cells, and their corresponding receptors on cancer-associated fibroblasts (CAFs) (Figure 7A). Annexin A1 (ANXA1) ligands, which originate from NK cells and CD8⁺ T cells, interact with their respective receptors on macrophage cells. In the ccRCC group, there is a significant decrease in the frequency of this receptor-ligand engagement (Figure 7B). Finally, the PARs and ANNEXIN signaling pattern in the ccRCC group were further investigated. Regarding the PARs signaling pattern, the results showed that CD8⁺ T cells in the ccRCC group engaged in communication with a spectrum of cellular elements through the GZMA-coagulation factor II thrombin receptor (F2R) and GZMA-F2R like thrombin or trypsin receptor 3 (F2RL3) ligand-receptor pairs (Figure 8A). The expression levels of receptors and ligands constituting the PARs signaling pattern, as articulated by disparate cell types within the ccRCC group, are graphically represented in Figure 8B. The elevated expression of F2R in CAFs within the ccRCC group elucidates the augmented PARs signaling between key immune cells and CAFs. With regard to the ANNEXIN signaling pattern, the results showed that CD8⁺ T cells in the ccRCC group engaged in multifaceted cellular communications via the ANXA1-formyl peptide receptor 1 (FPR1) ligand-receptor axis (Figure 8C). The levels of ANNEXIN receptor and ligand expression by different cell types in the ccRCC group are shown in Figure 8D. The diminished expression of ANXA1 accounts for the attenuated ANNEXIN signal observed between key immune cells and macrophages within the ccRCC group.

Figure 6 Number and intensity of cell communication in ccRCC group and control group. (A) Bar chart showing the number of interactions between cell types in the ccRCC group and the control group. (B) Bar chart showing the intensity of the interactions between cell types in the ccRCC group and control group. (C) Circle plot showing number of significant ligand-receptor pairs between any pair of two cell populations. Circle sizes are proportional to the number of cells in each cell group and edge width represents the number of interactions. Edge colors are consistent with the signaling source. (D) Circle plot showing intensity of significant ligand-receptor pairs between any pair of two cell populations. Circle sizes are proportional to the number of cells in each cell group and edge width represents the communication probability. Edge colors are consistent with the signaling source. (E) Heat maps delineating the signaling patterns that contribute the most to the overall signaling in the ccRCC group and the control group. CAF, cancer-associated fibroblast; ccRCC, clear cell renal cell carcinoma.

Figure 7 Intercellular ligand-receptor pairs between key immune cell types and other cell types. (A) The ligand-receptor pairs contribute to the increased signaling probability sending from CD8⁺ T cells and NK cells to other cell types in ccRCC. (B) The ligand-receptor pairs contribute to the decreased signaling probability sending from CD8⁺ T cells and NK cells to other cell types in ccRCC. ccRCC, clear cell renal cell carcinoma; NK, natural killer.

Figure 8 Key signaling patterns between key immune cell types and other cell types. (A) Ligand-receptor pairs of PARs signaling pattern in ccRCC. (B) The expression distribution of ligands and receptors of PARs signaling pattern in ccRCC group and control group. (C) Ligand-receptor pairs of ANNEXIN signaling pattern in ccRCC. (D) The expression distribution of ligands and receptors of ANNEXIN signaling pattern in ccRCC group and control group. CAF, cancer-associated fibroblast; ccRCC, clear cell renal cell carcinoma; PARs, protease activated receptors; PST, proximal straight tubule.

Construction and validation of prognostic model

A total of 2,246 DEGs were identified by comparing the ccRCC group and the control group from the bulk RNA-seq datasets (Figure S3A and available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-7.xlsx). Seven hundred and forty-two DEGs of key immune cells between ccRCC group and the control group were identified in the scRNA-seq datasets (available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-8.xlsx). The 175 overlapped genes of the two groups of DEGs were selected (Figure S3B, available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-9.xlsx). To explore the predictive value of the overlapped genes with respect to clinical outcomes in ccRCC, univariate Cox proportional hazards regression analysis of these genes was conducted using the TCGA cohort as a training set. As a result, 43 genes that were significantly associated with OS were obtained (available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-10.xlsx). To eliminate the effect of overfitting, LASSO analysis was employed. As a result, a subset of ten genes was identified and subsequently designated as the prognostic genes (available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-11.xlsx). The RS was derived from a formula incorporating the expression levels of the prognostic genes along with their respective coefficients: RS = PTTG1 × 0.041720319 + CLEC2D × 0.064191009 + PLIN2 × (−0.063698179) + LRBA × (−0.037693656) + ABCB1 × (−0.001437449) + MIR155HG × 0.005229499 + ADAM8 × 0.246465931 + P2RY8 × (−0.214111662) + SORL1 × (−0.1341849) + CD82 × 0.097135005. The “surv_cutpoint” function was employed to determine the best truncation value to distinguish high- and low-risk groups. The distribution of RS, survival time and expression levels of the prognostic genes in different groups of the training set were shown in Figure 9A. The high-risk group exhibited a significantly inferior prognosis. With the increasing of RS, the expression of PTTG1, CLEC2D, MIR155HG, ADAM8, and CD82 were observed to be upregulated, while the expression of PLIN2, LRBA, ABCB1, P2RY8, and SORL1 were noted to be downregulated. In the TCGA cohort, the median OS was significantly greater in the low-risk group as compared to the high-risk group (Figure 9B). The AUC values for predicting 1-, 3-, and 5-year survival in the TCGA cohort were 0.751, 0.734, and 0.747, respectively (Figure 9C). To assess the stability of the prognostic model developed from the training cohort with the prognostic genes, the E-MTAB-1980 cohort was introduced as an independent validation set. The RS for each patient in the validation set was calculated utilizing the identical formula that was applied in the TCGA cohort. For the validation cohort, patients in the high-risk group had poorer prognosis than those in the low-risk group (Figure 9D). Furthermore, based on the ROC analysis, the RS demonstrated robust predictive accuracy for OS in the validation cohort, with respective AUC values for 1-, 3-, and 5-year OS being 0.742, 0.740, and 0.778 (Figure 9E).

Figure 9 Construction and validation of key immune-cell-associated prognostic model. (A) RS distribution, survival status and heat map of the expression of 10 prognostic genes in the training cohort. (B) Survival curves of patients in high- and low-risk groups in the training cohort. (C) Time-dependent ROC curves of 1-, 3-, and 5-year of ccRCC patients in the training cohort. (D) Survival curves for patients in the high- and low-risk groups of the validation cohort. (E) Time-dependent ROC curves of 1-, 3-, and 5-year of ccRCC patients in the validation cohort. AUC, area under the curve; ccRCC, clear cell renal cell carcinoma; ROC, receiver operating characteristic; RS, risk score.

Construction of nomogram

To verify whether the RS could serve as an independent prognostic factor, univariate and multivariate Cox regression analyses were conducted, adjusting for clinical covariates including age, sex, and tumor stage. The results showed that RS and metastatic stage (M stage) were independent prognostic risk factors (Figure 10A,10B). Multivariate Cox regression analysis was employed to develop a nomogram, demonstrating that RS and M stage were significant predictors of clinical outcomes (Figure 10C). Calibration curves revealed that the nomogram model exhibited satisfactory stability and accuracy in predicting OS at the 1-, 3-, and 5-year milestones (Figure 10D).

Figure 10 Construction and validation of the nomogram. (A) Forest plot showing the result of univariate Cox regression analysis. (B) Forest plot showing the result of multivariate Cox regression analysis uncovered that the risk score is an independent prognostic risk factor for ccRCC patients. (C) Nomogram based on RS and M stage. (D) Calibration curve illustrating the consistency between predicted and observed 1-, 3-, and 5-year survival rates in ccRCC patients depending on the prognostic nomogram. ccRCC, clear cell renal cell carcinoma; CI, confidence interval; HR, hazard ratio; M stage, metastatic stage; N, node; RS, risk score; T, tumor.

Analysis of TMB and drug sensitivity

A spectrum of specific genetic mutations was scrutinized in ccRCC, and the top 20 genes exhibiting the highest mutation prevalence were graphically represented. von Hippel-Lindau (VHL) had the highest frequency of mutation in both groups, followed by polybromo 1 (PBRM1) (Figure 11A,11B). TMB is a key biomarker for predicting the success of immunotherapy. Therefore, an analysis of somatic mutations associated with ccRCC was conducted, revealing no significant difference in TMB between the high and low risk groups (Figure 11C).

Figure 11 Difference of TMB and drug sensitivity between high- and low-risk group. (A,B) Heatmap showing the top 20 genes with the highest mutation frequency in the high- and low-risk group. (C) TMB differences between the high and low risk groups. (D-L) Differences of sensitivity to therapy drugs between the high- and low-group. TMB, tumor mutation burden.

We analyzed whether the RS could accurately predict drug sensitivity in patients afflicted with ccRCC. The clinical efficiency of Afatinib_1032, Dihydrorotenone_1827, BMS.754807_2171, BI.2536_1086, Osimertinib_1919, Buparlisib_1873, AZD7762_1022, Sabutoclax_1849, and Topotecan_1808 (Figure 11D-11L) in the treatment of ccRCC were studied (available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-12.xlsx). Patients with low RS were likely to be more sensitive to drugs such as Afatinib_1032, Dihydrorotenone_1827, BMS.754807_2171, BI.2536_1086, Osimertinib_1919, and Buparlisib_1873, While, patients with high RS showed higher sensitive to drugs such as AZD7762_1022, Sabutoclax_1849, and Topotecan_1808.

A comprehensive screening identified a total of 244 pharmacological agents with potential interactions with the 10 prognostic genes implicated in ccRCC within the DGIdb database (available online: https://cdn.amegroups.cn/static/public/tcr-2025-971-13.xlsx). Among these, a select subset of five drugs exhibited an interaction score surpassing the threshold of 10, indicating a more significant pharmacogenomic association (Table 1). Among the five drugs, YL-365 and LYSOPHOSPHATIDYLSERINE were the drugs targeting gene A disintegrin and metalloproteinase domain 8 (ADAM8), and RTI-122, RTI-13951-33 and COMPOUND 2 were the drugs targeting the gene perilipin 2 (PLIN2) (15).

Validation of the expression levels of prognostic genes in clinical samples

Immunohistochemistry results of the protein expression of the prognostic genes from HPA database are displayed in Figure 12. Pituitary tumor-transforming gene 1 (PTTG1) demonstrates low expression in normal renal tissues, with a range from low to moderate expression in ccRCC tumor tissues. C-type lectin domain family 2 member D (CLEC2D) is moderately expressed in normal renal tissues, yet its expression is predominantly low or undetectable in tumor tissues. PLIN2 is characterized by low expression in normal renal tissues, in contrast to its predominantly high expression in tumor tissues. LPS-responsive beige-like anchor protein (LRBA) shows high expression in normal renal tissues, with a variable expression pattern in tumor tissues, spanning from undetectable to low and high expression levels. ATP-binding cassette subfamily B member 1 (ABCB1) exhibits medium and low expression levels in normal renal and tumor tissues, respectively. ADAM8 is moderately expressed in normal renal tissues, but its expression in tumor tissues is heterogeneous, ranging from undetectable to low and medium levels. P2Y receptor family member 8 (P2RY8) is not detected in both normal and tumor tissues. Sortilin-related receptor 1 (SORL1) is highly expressed in normal renal tissues, yet its expression in tumor tissues is diverse, varying from undetectable to low and medium levels. CD82 is undetectable in normal renal tissues but exhibits medium expression in tumor tissues. Data for MIR155 host gene (MIR155HG) were unavailable and thus not presented in the analysis.

Figure 12 Immunohistochemistry of prognostic genes in ccRCC and normal samples from the HPA database. IHC staining of prognostic genes in ccRCC and normal samples from the HPA database. Image credit goes to the HPA (https://www.proteinatlas.org/). The links to the staining of proteins encoded by prognostic genes in ccRCC tissues and corresponding control tissues are as follows: PTTG1_Tumor (https://www.proteinatlas.org/ENSG00000164611-PTTG1/cancer/renal+cancer#img), PTTG1_Normal (https://www.proteinatlas.org/ENSG00000164611-PTTG1/tissue/kidney#img), CLEC2D_Tumor (https://www.proteinatlas.org/ENSG00000069493-CLEC2D/cancer/renal+cancer#img), CLEC2D_Normal (https://www.proteinatlas.org/ENSG00000069493-CLEC2D/tissue/kidney#img), PLIN2_Tumor (https://www.proteinatlas.org/ENSG00000147872-PLIN2/cancer/renal+cancer#img), PLIN2_Normal (https://www.proteinatlas.org/ENSG00000147872-PLIN2/tissue/kidney#img), LRBA_Tumor (https://www.proteinatlas.org/ENSG00000198589-LRBA/cancer/renal+cancer#img), LRBA_normal (https://www.proteinatlas.org/ENSG00000198589-LRBA/tissue/kidney#img), ABCB1_Tumor (https://www.proteinatlas.org/ENSG00000085563-ABCB1/cancer/renal+cancer#img), ABCB1_Normal (https://www.proteinatlas.org/ENSG00000085563-ABCB1/tissue/kidney#img), ADAM8_Tumor (https://www.proteinatlas.org/ENSG00000151651-ADAM8/cancer/renal+cancer#img), ADAM8_Normal (https://www.proteinatlas.org/ENSG00000151651-ADAM8/tissue/kidney#img), P2RY8_Tumor (https://www.proteinatlas.org/ENSG00000182162-P2RY8/cancer/renal+cancer#img), P2RY8_Normal (https://www.proteinatlas.org/ENSG00000182162-P2RY8/tissue/kidney#img), SORL1_Tumor (https://www.proteinatlas.org/ENSG00000137642-SORL1/cancer/renal+cancer#img), SORL1_Normal (https://www.proteinatlas.org/ENSG00000137642-SORL1/tissue/kidney#img), CD82_Tumor (https://www.proteinatlas.org/ENSG00000085117-CD82/cancer/renal+cancer#img), CD82_Normal (https://www.proteinatlas.org/ENSG00000085117-CD82/tissue/kidney#img). Scale bar, 200 µm. ABCB1, ATP-binding cassette subfamily B member 1; ADAM8, A disintegrin and metalloproteinase domain 8; ccRCC, clear cell renal carcinoma; CLEC2D, C-type lectin domain family 2 member D; HPA, Human Protein Atlas; IHC, immunohistochemical; LRBA, LPS-responsive beigelike anchor protein; PLIN2, perilipin 2; PTTG1, Pituitary tumor-transforming gene 1; P2RY8, P2Y receptor family member 8; SORL1, Sortilin-related receptor 1.

Discussion

In this study, we identified seven immune cells infiltrated in the TME of ccRCC, namely B lymphocytes, CD4⁺ T cells, CD8⁺ T cells, macrophages, mast cells, monocytes, and NK cells, from the scRNA-seq datasets. When deconvolving the scRNA-seq expression profiles to the bulk RNA-seq dataset, we observed that the infiltration levels of the seven immune cell types elevated, while the proportion of renal tubules and endothelium, among other cellular components, notably decreased in the ccRCC group relative to the healthy control group. This observation reflects the heightened inflammatory activity and concomitant loss of renal tubular structure in ccRCC and aligns with the immunogenic characteristic of ccRCC (16). Notably, CD8⁺ T cells and NK cells were verified as the key immune cells in ccRCC with three machine-learning algorithms. CD8⁺ T cells and NK cells can be induced by various cytokines to differentiate into distinct polarized states. To provide a global view of the cellular responses of each immune cell type to each cytokine, Cui and his colleagues created immune dictionary and developed companion software, IREA, for assessing cytokine activities and immune cell polarization from gene expression data. And IREA classified CD8⁺ T cells into T8-a, T8-b, T8-c, T8-e, and T8-f subsets, and NK cells into NK-a, NK-c, NK-e, and NK-f polarized states, using DEGs from scRNA-seq. In our study, IREA automatically inferred that CD8⁺ T cells and NK cells had polarized into the T8-c and NK-f states, respectively. From a traditional perspective, effector CD8⁺ T cells are typically perceived as a homogenous population of cytotoxic lymphocytes that secret IFN-γ and protease granzyme B. However, recent findings by Michael St. Paul and colleagues (17) delineated multiple subsets of CD8⁺ T cells within TME, each possessing distinct effector functions and varying levels of cytotoxic potential. These subsets were referred to as the Tc subsets. Similarly, in our study, the CD8⁺ T cells with a T8-c polarization state were induced mainly by IFN-α1 and IL-36α. Functionally, T8-c polarized CD8⁺ T cells were characterized by elevated levels of transcription factor 7 (TF7) and interferon gamma receptor 1 (IFNGR1) expression, both of which are pivotal in the differentiation of T-lymphocytes. The induction of NK cells into the NK-f polarization phenotype is mediated by IL-18, a cytokine for which the receptor is highly expressed on the surface of NK cells. IL-18 can stimulate NK cells to secrete IFN-γ and exert broad immunomodulatory effects (18). And NK cells exhibiting this polarization phenotype are characterized by elevated expression of MYC proto-oncogene (MYC), which is essential for the development and proliferation of NK cells (19). Beyond evaluating immune cell polarization, IREA also revealed the pivotal role of cytokines in immune modulation. IL-15, IFN-α1, IL-2, and IFN-β were significantly enriched in NK cells and CD8⁺ T cells within ccRCC, indicating an activated immune state in ccRCC. These cytokines play crucial roles in the pathogenesis of ccRCC by enhancing antitumor immune functions, activating T lymphocytes and NK cells, inhibiting tumor cell proliferation, and promoting tumor cell death.

The intercellular communications among immune cells had also undergone changes accompanying the alterations of immune cells and cytokines in ccRCC samples. The cell communication analysis revealed that the PARs signaling pattern derived from CD8⁺ T cells and NK cells to CAFs via GZMA-F2R ligand-receptor pair elevated in comparison to the PARs signaling from CD8⁺ T cells and NK cells to normal kidney-associated fibroblasts. Previous study demonstrated that GZMA secreted by cytotoxic T lymphocytes and NK-cells was demonstrated to interact with F2R expressed by the tumor cells both in vivo and in vitro (20). This communication process is also intimately associated with cytokines such as IL-15, IFN-α1, IL-2, and IFN-β. IL-15 significantly amplifies the proliferation and activation of NK cells and CD8⁺ T cells, thereby augmenting the cytotoxic functions of these immune cells and leading to an increase in the secretion of GZMA. This enhancement may, in turn, boost the F2R signal transduction, facilitating the migration of immune cells into the TME and intensifying the antitumor immune response. Similarly, IFN-α1 and IFN-β are capable of activating immune cells, including dendritic cells and NK cells, and may promote the secretion of GZMA, thereby enhancing the cytotoxic effects of immune cells. These cytokines may also modulate the function of CAFs through F2R, altering their immunomodulatory effects and subsequently impacting the inflammatory state of the TME, which can either support or suppress tumor growth. IL-2, a key cytokine for T cell proliferation, further intensifies the immune response by enhancing the secretion of GZMA from activated CD8+ T cells, potentially influencing the interaction with CAFs through F2R. This intricate dynamic communication between immune cells and tumor stroma underscores the pivotal role of PARs signaling pattern in regulating tumor immunity. Accordingly, the therapeutic implications are substantial, as enhancing the interaction between immune cells and the stroma may help overcome mechanisms of immune evasion. Furthermore, the combined use of IL-15, IFN-α1, IL-2, and IFN-β with GZMA or F2R targeted therapies may provide new opportunities to improve immune checkpoint blockades (ICBs) therapeutic efficacy, opening new avenues for more effective cancer immunotherapy. Conversely, the ANNEXIN signal derived from CD8+ T cells and NK cells to macrophages via ANXA1-FPR1 ligand-receptor pair was observed to be diminished. ANXA1 is a calcium-dependent phospholipid binding protein with a molecular weight of 37 kDa. It consists of two different regions, the tail of the N-terminus and the core domain of the C-terminus (21). The core domain has Ca²⁺ and cell membrane binding sites, which are common to the annexin family, while each member of the N-terminus is different and mainly regulates the interaction of ANXA1 with ligands. ANXA1 works mainly through N-terminal binding to other molecules, and its unique structure is the basis for many biological activities within cells (22). Under physiological conditions, ANXA1 facilitates immunosuppression and inhibits inflammatory processes, potentially due to its capacity to augment the differentiation of anti-inflammatory (M2) macrophages and promote IL-10 expression (23). And ANXA1 can exert its effects on the FPR1 to induce the activation of macrophages toward an M2-type polarization. Subsequently, in our study, the diminished ANNEXIN signaling indicated a promoted immune response in ccRCC.

Furthermore, we developed a prognostic signature comprising 10 prognostic genes (i.e., PTTG1, CLEC2D, PLIN2, LRBA, ABCB1, MIR155HG, ADAM8, P2RY8, SORL1, and CD82) to determine the RS in the TCGA cohort. Upon stratifying the ccRCC patients into high- and low-risk groups based on the RS, we observed a robust predictive efficacy of the prognostic signature in the independent validation cohort E-MTAB-1980. PTTG1 is an oncogene localized to the chromosome 5q region, which is frequently amplified in ccRCC. PTTG1 has the potential to directly induce carcinogenesis by cell transformation, activating proto-oncogenes, and growth factors. Moreover, overexpression of PTTG1 in ccRCC tissue is associated with high Fuhrman grade, high tumor stage, signifying more aggressive and invasive tumor behaviors (24). CLEC2D, the ligand for CD161, is a surface molecule expressed by myeloid cells and malignant cells. According to the single-cell analysis conducted by Mathewson and colleagues (25), killer cell lectin like receptor B1 (KLRB1), which encodes CD161, is expressed by CD8⁺ T cells with high cytotoxicity in diffuse gliomas. In vitro and in vivo experiments demonstrated that the CD161 receptor can inhibit key functions of T lymphocytes, including cytotoxicity and cytokine secretion. Thus, the CLEC2D-CD161 pathway exhibits certain parallels with the programmed cell death protein 1 (PD-1) pathway in terms of its immunomodulatory effects. Despite the unknown role of IL-18, the genes encoding IL-18 receptor were expressed in KLRB1-expressing T cells in diffuse gliomas. Notably, the IREA analysis in our study also enriched IL-18 and T8-c polarized CD8⁺ T cells. Then the immune suppression mechanism of CLEC2D-CD161 pathway deserves further investigation in ccRCC. PLIN2 is involved in lipid droplet formation. Prior research demonstrated lipid metabolism may play an important role in supporting inflammatory conditions in TME, analogous to the function of classically activated, pro-inflammatory (M1) macrophages (26). When treated with cancer vaccines, macrophages positive for PLIN2 upregulate genes associated with interferon signaling, including interferon regulatory factor 7 (IRF7), C-X-C motif chemokine ligand 2 (CXCL2), interferon induced transmembrane protein 1 (IFITM1) and ISG15 Ubiquitin Like Modifier (ISG15), thus promoting anti-tumor efficacy in TME (27). LRBA is closely related to lysosomal trafficking regulator (LYST) and is known to be an intracellular adaptor protein. LYST is a molecule that regulates the fusion of intracellular vesicles such as lysosomes. LYST deficiency causes Chediak-Higashi Syndrome, a disease that is related to defects in vesicular trafficking in lymphocytes. In cytotoxic lymphocytes, secretory lysosomes are responsible for the transport of the toxic cargo to the cell membrane, where they fuse, delivering their content directly to their target tumor cells (28). ABCB1 encodes a membrane-associated protein which is a member of the multidrug resistance (MDR)/transporter associated with antigen processing (TAP) subfamily belonging to ATP-binding cassette transporters. Members of the MDR/TAP subfamily are involved in multidrug resistance for decreasing drug accumulation in multidrug-resistant cells and mediating the development of resistance to anticancer drugs. Previous study demonstrated that IFN-γ induced downregulation of ABCB1 through the mediation of interferon regulatory factor-1, thereby reversing the multidrug resistant phenotype in gastric cancer cells (29). MIR155HG is a precursor gene encoding miR-155 and is highly expressed in a variety of immune cells, such as CD8⁺ T cells and NK cells. The long RNA transcribed from this gene is expressed at high level in lymphoma and may function as an oncogene. In hepatocellular carcinoma cells, N6-methyladenosine (m6A) modification of MIR155HG upregulates programmed cell death ligand 1 (PD-L1) expression facilitating immune escape (30). ADAM8 encodes a member of the disintegrin and metalloprotease domain family. Members of this family are membrane-anchored proteins structurally related to a variety of biological processes involving cell-cell and cell-matrix interactions. Many studies focus on the prognostic prediction value of ADAM8. For instance, Qu and colleagues (31) found ADAM8 promoted the proliferation, migration, invasion, and tumorigenesis of RCC cells, while, Roemer and colleagues identified ADAM8 as a best predictor of distant metastasis in ccRCC (32). Thus, our findings are consistent with the evidence in the literature. P2RY8 encodes a G protein–coupled receptor that is involved in restraining germinal center B lymphocytes migration and growth. Loss-of-function of P2RY8 in human T cells increased their bone marrow homing (33). According to the latest study by Wang and colleagues (34), dietary fasting drove T cell homing to bone marrow and suppressed T cell activation, proliferation, differentiation and cytokine production in autoimmune mouse models and substantially alleviated disease symptoms. Given the observed negative correlation of P2RY8 with the RS, we hypothesize that P2RY8 may contribute to promoting antitumor immune responses in ccRCC which warrants further investigation. SORL1 facilitates sorting and transport between the Golgi apparatus and the cell surface through an endosomal separator. Previous studies have documented that SORL1 expression is significantly increased in breast cancer and bladder cancer and is predictive of poor prognosis (35). Interestingly, downregulation of SORL1 expression can reverse chemoresistance in ovarian cancer cells by inhibiting the ABCB1 and early endosomal antigen 1 pathway (36). CD82 is a metastasis suppressor gene product. Expression of this gene has been demonstrated to be downregulated during tumor progression in human cancers and can be activated by P53 through a consensus binding sequence present in its promoter region (37). Additionally, CD82 augments NK cell-mediated cytotoxicity. Then CD82 can work synergistically with IL-18 to bolster antitumor immune responses.

In terms of prediction of potential anticancer drugs, the prognostic model played a crucial role. Given the intrinsic resistance to conventional chemotherapy and radiotherapy, significant advancements have been achieved in the treatment of metastatic ccRCC, with the approval of targeted therapy agents and ICBs for clinical use. Particularly, in the new era of immune therapy, the establishment of predictors for ICBs is thus of utmost importance to identify patient populations more likely to respond to such therapies, thus maximizing therapeutic benefits. TMB has been identified as one of the most important predictors for responders to ICBs. However, analysis in our study revealed no difference in TMB between high- and low-risk groups. This finding may align with the consensus that ccRCC is a neoplasm with low TMB, despite its immunogenic characteristics. To some extent, the outcomes of drug sensitivity analysis predicated on the RS from our study aligned with the highly esteemed International mRCC Database Consortium (IMDC) risk score employed in clinical settings, as the low-risk group exhibited heightened sensitivity to tyrosine kinase inhibitors, while the high-risk group showed increased sensitivity to a B-cell lymphoma-2 inhibitor. Furthermore, we identified YL-365 as a pharmaceutical agent targeting ADAM8, one of the ten prognostic genes, with the most significant interaction score as appraised by the DGIdb. Previous investigations have established that the overexpression of G protein-coupled receptor 34 (GPR34) influences the progression and prognosis of human gastric adenocarcinoma and colorectal cancer through the phosphatidylinositol 3-kinase (PI3K)/protein kinase B (AKT) signaling pathway (38). As a potent and selective antagonist of GPR34, YL-365 may be a candidate agent in the treatment of ccRCC.

The present study is not without its limitations. Firstly, all data analyzed in this study were derived from public databases. The polarization states of key immune cells and the prediction value of the prognostic model were not validated in our own clinical practice prospectively. Then, there is a lack of mechanistic research of the ten genes in the prognostic model. Nonetheless, the findings we report possess unprecedentedly substantial implications for elucidating the functional alterations of key immune cells and their associations with prognostic prediction in ccRCC. Future validation of our results in large-scale, prospective clinical and experimental studies are warranted to substantiate their clinical relevance and applicability.

Conclusions

In conclusion, we have identified the key immune cells infiltrated in ccRCC by integrating scRNA-seq and bulk RNA-seq datasets. Then, the functional alterations of the key immune cells indicated a promoted immune response in ccRCC. In addition, we constructed a prognostic signature based on DEGs associated with key immune cells, which has been shown to provide accurate survival and antitumor drug sensitivity predictions for ccRCC patients. Future research should encompass large-scale prospective clinical studies to confirm the functional states of key immune cells and to verify the interactions and alterations between these cells and other cellular components within the TME. Additionally, the predictive capacity of the prognostic model should be validated in clinical practice to ensure its utility and reliability in patient management.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-971/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-971/prf

Funding: This work was supported by the National Natural Science Foundation of China (No. 82172759 to S.Z.) and Tianjin health research project (No. TJWJ2024ZD002 to C.Q.).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-971/coif). S.Z. receives funding from the National Natural Science Foundation of China (No. 82172759). C.Q. receives funding from Tianjin health research project (No. TJWJ2024ZD002). The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Siegel RL, Giaquinto AN, Jemal A. Cancer statistics, 2024. CA Cancer J Clin 2024;74:12-49. [Crossref] [PubMed]
2. Blankenstein T, Coulie PG, Gilboa E, et al. The determinants of tumour immunogenicity. Nat Rev Cancer 2012;12:307-13. [Crossref] [PubMed]
3. Tong L, Kremer V, Neo SY, Seitz C, et al. Cellular and secretome profiling uncover immunological biomarkers in the prognosis of renal cell carcinoma patients. Oncoimmunology 2025;14:2481109. [Crossref] [PubMed]
4. Núñez SY, Trotta A, Regge MV, et al. Tumor-associated macrophages impair NK cell IFN-γ production and contribute to tumor progression in clear cell renal cell carcinoma. Eur J Immunol 2024;54:e2350878. [Crossref] [PubMed]
5. Chevrier S, Levine JH, Zanotelli VRT, et al. An Immune Atlas of Clear Cell Renal Cell Carcinoma. Cell 2017;169:736-749.e18. [Crossref] [PubMed]
6. Mei S, Alchahin AM, Tsea I, et al. Single-cell analysis of immune and stroma cell remodeling in clear cell renal cell carcinoma primary tumors and bone metastatic lesions. Genome Med 2024;16:1. [Crossref] [PubMed]
7. Eckl J, Buchner A, Prinz PU, et al. Transcript signature predicts tissue NK cell content and defines renal cell carcinoma subgroups independent of TNM staging. J Mol Med (Berl) 2012;90:55-66. [Crossref] [PubMed]
8. Li X, Zhang Y, Ye Y, et al. NK cells in renal cell carcinoma and its implications for CAR-NK therapy. Front Cell Dev Biol 2025;13:1532491. [Crossref] [PubMed]
9. Ye Z, Zhang Y, Xu J. Integrating Bulk and Single-Cell RNA-Seq Data to Identify Prognostic Features Related to Activated Dendritic Cells in Clear-Cell Renal-Cell Carcinoma. Int J Mol Sci 2024;
10. Wu J, Jin S, Gu W, et al. Construction and Validation of a 9-Gene Signature for Predicting Prognosis in Stage III Clear Cell Renal Cell Carcinoma. Front Oncol 2019;9:152. [Crossref] [PubMed]
11. Rini B, Goddard A, Knezevic D, et al. A 16-gene assay to predict recurrence after surgery in localised renal cell carcinoma: development and validation studies. Lancet Oncol 2015;16:676-85. [Crossref] [PubMed]
12. Sato Y, Yoshizato T, Shiraishi Y, et al. Integrated molecular analysis of clear-cell renal cell carcinoma. Nat Genet 2013;45:860-7. [Crossref] [PubMed]
13. Cui A, Huang T, Li S, et al. Dictionary of immune responses to cytokines at single-cell resolution. Nature 2024;625:377-84. [Crossref] [PubMed]
14. Davidson G, Helleux A, Vano YA, et al. Mesenchymal-like Tumor Cells and Myofibroblastic Cancer-Associated Fibroblasts Are Associated with Progression and Immunotherapy Response of Clear Cell Renal Cell Carcinoma. Cancer Res 2023;83:2952-69. [Crossref] [PubMed]
15. Jin C, Decker AM, Huang XP, et al. Synthesis, pharmacological characterization, and structure-activity relationship studies of small molecular agonists for the orphan GPR88 receptor. ACS Chem Neurosci 2014;5:576-87. [Crossref] [PubMed]
16. Díaz-Montero CM, Rini BI, Finke JH. The immunology of renal cell carcinoma. Nat Rev Nephrol 2020;16:721-35. [Crossref] [PubMed]
17. St Paul M, Ohashi PS. The Roles of CD8(+) T Cell Subsets in Antitumor Immunity. Trends Cell Biol 2020;30:695-704. [Crossref] [PubMed]
18. Shen J, Zhang Y, Tang W, et al. Short IL-18 generated by caspase-3 cleavage mobilizes NK cells to suppress tumor growth. Nat Immunol 2025;26:416-28. [Crossref] [PubMed]
19. Khameneh HJ, Fonta N, Zenobi A, et al. Myc controls NK cell development, IL-15-driven expansion, and translational machinery. Life Sci Alliance 2023;
20. Gao Y, Xu Q, Li X, et al. Heterogeneity induced GZMA-F2R communication inefficient impairs antitumor immunotherapy of PD-1 mAb through JAK2/STAT1 signal suppression in hepatocellular carcinoma. Cell Death Dis 2022;13:213. [Crossref] [PubMed]
21. Dudas EF, Tully MD, Foldes T, et al. The structural properties of full-length annexin A11. Front Mol Biosci 2024;11:1347741. [Crossref] [PubMed]
22. Gerke V, Gavins FNE, Geisow M, et al. Annexins-a family of proteins with distinctive tastes for cell signaling and membrane dynamics. Nat Commun 2024;15:1574. [Crossref] [PubMed]
23. Li Y, Cai L, Wang H, et al. Pleiotropic regulation of macrophage polarization and tumorigenesis by formyl peptide receptor-2. Oncogene 2011;30:3887-99. [Crossref] [PubMed]
24. Zhang X, Ji H, Huang Y, et al. Elevated PTTG1 predicts poor prognosis in kidney renal clear cell carcinoma and correlates with immunity. Heliyon 2023;9:e13201. [Crossref] [PubMed]
25. Mathewson ND, Ashenberg O, Tirosh I, et al. Inhibitory CD161 receptor identified in glioma-infiltrating T cells by single-cell analysis. Cell 2021;184:1281-1298.e26. [Crossref] [PubMed]
26. Rosas-Ballina M, Guan XL, Schmidt A, et al. Classical Activation of Macrophages Leads to Lipid Droplet Formation Without de novo Fatty Acid Synthesis. Front Immunol 2020;11:131. [Crossref] [PubMed]
27. Baharom F, Ramirez-Valdez RA, Khalilnezhad A, et al. Systemic vaccination induces CD8(+) T cells and remodels the tumor microenvironment. Cell 2022;185:4317-4332.e15. [Crossref] [PubMed]
28. Trujillo-Vargas CM. Editorial: Intracellular trafficking in lymphocytes: the role of inner crowds in determining cell fate. Front Immunol 2023;14:1212138. [Crossref] [PubMed]
29. Yuan J, Yin Z, Tan L, et al. Interferon regulatory factor-1 reverses chemoresistance by downregulating the expression of P-glycoprotein in gastric cancer. Cancer Lett 2019;457:28-39. [Crossref] [PubMed]
30. Peng L, Pan B, Zhang X, et al. Lipopolysaccharide facilitates immune escape of hepatocellular carcinoma cells via m6A modification of lncRNA MIR155HG to upregulate PD-L1 expression. Cell Biol Toxicol 2022;38:1159-73. [Crossref] [PubMed]
31. Qu H, Mao M, Wang K, et al. Knockdown of ADAM8 inhibits the proliferation, migration, invasion, and tumorigenesis of renal clear cell carcinoma cells to enhance the immunotherapy efficacy. Transl Res 2024;266:32-48. [Crossref] [PubMed]
32. Chu MC, Zhang ZL, Wang ZQ, et al. Potential prognostic and immunologic significances of ADAM8 in clear cell renal cell carcinoma. Medicine (Baltimore) 2025;104:e41375. [Crossref] [PubMed]
33. Gallman AE, Wolfreys FD, Nguyen DN, et al. Abcc1 and Ggt5 support lymphocyte guidance through export and catabolism of S-geranylgeranyl-l-glutathione. Sci Immunol 2021;
34. Wang L, Cheng M, Wang Y, et al. Fasting-activated ventrolateral medulla neurons regulate T cell homing and suppress autoimmune disease in mice. Nat Neurosci 2024;27:462-70. [Crossref] [PubMed]
35. Pietilä M, Sahgal P, Peuhu E, et al. SORLA regulates endosomal trafficking and oncogenic fitness of HER2. Nat Commun 2019;10:2340. [Crossref] [PubMed]
36. Zheng Z, Li X, Yang B, et al. SORL1 stabilizes ABCB1 to promote cisplatin resistance in ovarian cancer. Funct Integr Genomics 2023;23:147. [Crossref] [PubMed]
37. Lee JW, Hur J, Kwon YW, et al. KAI1(CD82) is a key molecule to control angiogenesis and switch angiogenic milieu to quiescent state. J Hematol Oncol 2021;14:148. [Crossref] [PubMed]
38. Yu W, Ma S, Wang L, et al. Upregulation of GPR34 expression affects the progression and prognosis of human gastric adenocarcinoma by PI3K/PDK1/AKT pathway. Histol Histopathol 2013;28:1629-38. [Crossref] [PubMed]