RNA-sequencing data, clinical information, and survival data for 524 ccRCC patients and 72 control samples were obtained from The Cancer Genome Atlas (TCGA-KIRC) database. The GEO dataset GSE29609 provided RNA sequencing data for 39 ccRCC patients. We also used the GSE40435 dataset (101 ccRCC and 101 control samples) to validate key findings. To form a unified discovery cohort, TCGA-KIRC and GSE29609 were merged after normalizing expression data to transcripts per million (TPM) values and applying batch effect correction using the ComBat algorithm from the ‘sva’ R package (version 3.42.0). This resulted in a combined cohort of 563 ccRCC samples for primary analyses, including consensus clustering. All processing was conducted in R (version 4.2.2) using the limma (version 3.50.0) and sva (version 3.42.0) packages. Input files were read as tab-delimited matrices with gene symbols in the first column. Common genes across datasets were identified via intersection (resulting in 15,742 genes). For each dataset: (1) Duplicate genes were averaged using avereps from limma. (2) For TCGA-KIRC, only tumor samples were retained by filtering barcodes (fourth digit = 0), duplicate samples averaged if present, and data log2-transformed if quantile checks indicated raw-scale values (qx[5] > 100 or (qx[6] - qx[1]) > 50 & qx[2] > 0), using log2(x + 1). (3) For GSE29609, quantile normalization was applied via normalizeBetweenArrays from limma. Matrices were then merged on intersecting genes, with dataset origin labeled as batch factors. Batch effects were corrected using the ComBat algorithm from sva, with parametric priors (par.prior = TRUE) and dataset as the batch covariate. Quality control metrics included quantile distributions inspected to confirm normalization needs (e.g., median qx values ∼6–10 post-log for TCGA, ∼8–12 for GEO), genes with zero variance or missing in >50% samples filtered via intersection and avereps (∼20% removal), and post-ComBat principal component analysis (PCA) via prcomp confirming batch removal (PC1 now 24.0% variance, PC2 6.6%, together explaining 30.7% variance, with TCGA and GSE29609 samples largely intermixed and clustering by tumor biology; adjusted Rand index pre- vs. post: ∼0.15). This aligns with the workflow depicted in Figure 1, where PCA supports the identification of PANoptosis clusters in ccRCC by validating data integration prior to consensus clustering and subsequent analyses.

The Wilcoxon rank-sum test was applied to identify differentially expressed PRGs between ccRCC and normal tissues. Prognosis-related PRGs were identified through univariate Cox regression and Kaplan-Meier survival analyses.

Based on 29 PRGs from prior literature (Samir et al., 2020; Wang and Kanneganti, 2021; Wang X. et al., 2022), consensus clustering analysis was conducted using with the “ConsensusClusterPlus” R package to classify ccRCC patients into PANoptosis subtypes. Clustering analysis was performed over 1,000 iterations to ensure robustness, and subtypes were validated using principal component analysis (PCA).

GSVA was conducted with the “GSVA” R package, using the HALLMARK gene set “h.all.v2022.1.Hs.symbols.gmt” downloaded from the MSigDB database (https://www.gsea-msigdb.org/gsea/index.jsp) to explore biological differences among the subtypes (Liberzon et al., 2015).

In addition, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the “org.Hs.eg.db” and “enrichplot” R packages to identify relevant biological functions and pathways. An adjusted p value <0.05 indicates statistical significance.

Differentially expressed genes (DEGs) were identified between PANoptosis subtypes in pairwise comparisons. DEGs with a log2(fold change) > 0.585 and adjusted p value <0.05 were employed to identify gene signatures. Prognosis-related significant genes were identified through univariate Cox regression analysis. To identify additional PRGs for signature development, we divided the patients into two distinct gene clusters based on the expression of DEGs. To minimize the least absolute shrinkage between signatures, we performed LASSO Cox regression analysis using the univariate significant genes. Subsequently, multivariate Cox regression analysis was conducted to select the best gene and estimate the correlation coefficient. The PRG risk score can be calculated using the following formula: PANoptosis score = ∑ i = 1 n exp ⁡ Xi × coef Xi . In this formula, exp (Xi) represents the expression level of the corresponding genes, and coef (Xi) represents the coefficient of the three genes. Patients were stratified into high- and low-risk groups according to the median risk score, and Kaplan-Meier analysis was performed to evaluate survival differences. Predictive accuracy was validated using time-dependent ROC curves, and a nomogram was developed for 1-, 2-, and 3-year survival prediction, incorporating clinical characteristics and risk groups. Calibration plots were utilized to validate the consistency between the predicted and observed survival rates of patients with ccRCC.

Single-sample gene set enrichment analysis (ssGSEA) and the CIBERSORT algorithm were employed to evaluate immune cell infiltration, while the ESTIMATE algorithm quantified stromal and immune scores within the TME which refers to the complex, dynamic environment surrounding a tumor.

Gene mutation data were analyzed using the “maftools” R package to explore mutation profiles in different groups of ccRCC. Additionally, the relationship between risk scores and the cancer stem cell index was also investigated.

Drug sensitivity predictions were made using the “oncoPredict” R package, based on IC50 values from the Genomics of Drug Sensitivity in Cancer v2(GDSC) database for 198 drugs (https://www.cancerrxgene.org/). P-values <0.05 were considered statistically significant.

A total of 22 pairs of ccRCCs and their adjacent normal samples were collected from Liuzhou Workers’ Hospital, with approval from the Ethics Committee (KY2023412). Total RNA was isolated from the human samples, normal renal cells (Hk-2), and ccRCC cells (Caki-1, Caki-2, A498, 786-o) using TRIzol reagent (Thermo Fisher Scientific, United States) following the manufacturer’s instructions. Reverse transcription was carried out using the PrimeScriptTM RT Reagent Kit (Takara, Japan). Quantitative reverse transcription polymerase chain reaction (qRT‒PCR) was performed on an FX Connect system (Bio-Rad, United States) using SYBR® Green Supermix (Bio-Rad, United States) to measure the expression levels of hub genes. The primer sequences used for qRT‒PCR are provided in Supplementary Table S1. Unpaired t tests were used to assess the differential expression levels between ccRCC and adjacent noncancer tissues. β-Actin was used as the reference gene.

All statistical tests were conducted using R version 4.2.2. For two-group comparisons, the Student’s t-test or Wilcoxon test was applied based on data distribution, while the Kruskal–Wallis test was used for multiple group comparisons. Statistical significance was defined as p < 0.05. In the figures, the asterisks reflect the statistical P value (ns, P ≥ 0.05, *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001).

A total of 29 PRGs were included in this analysis. We compared the differential expression of PRGs in ccRCC tumor samples and healthy samples collected from the TCGA-KIRC database. The results indicated that CASP8, FADD, CASP6, NLRP3, PSTPIP2, TNFAIP3, GSDMD, MLKL, IRF1, AIM2, ZBP1, CASP1, RIPK1, RIPK3, TRADD, MEFV, PYCARD, and NLRC4 were upregulated in ccRCC samples compared to normal samples. However, only TAB2 and TAB3 were downregulated in ccRCC (Figure 2A).

Somatic mutation analysis of the PRGs revealed a low mutation frequency in the ccRCC samples. Only 17 (4.23%) of 402 ccRCC samples exhibited somatic mutations (Figure 2B). Additionally, somatic copy number variation (CNV) analysis was performed on these PRGs, revealing frequent CNV alterations. Among them, GSDMD, FADD, NLRP3, ZBP1, TNFAIP3, AIM2, PARP1, TAB2, PYCARD, CASP8, RIPK3, MLKL, CASP6, and TRADD showed increased CNV, whereas CASP7, RIPK1, MEFV, IRF1, NLRC4, PSTPIP2, TAB3, and CASP1 showed relatively decreased CNV (Figure 2C). The specific locations of the CNV alterations in PRGs on the chromosomes of PRGs are shown in Figure 2D.

A total of 524 CCRCC patients from TCGA and 39 patients from the GEO database (GSE29609) were enrolled to explore the expression pattern of PRGs in ccRCC. Kaplan‒Meier curves revealed that 18 PRGs were associated with the overall survival rate (OS) of ccRCC. High expression of AIM2, GSDMD, MEFV, PSTPIP2, PYCARD, RIPK3, and ZBP1 was associated with poor prognosis, whereas high expression of CASP1, CASP6, CASP7, FADD, NLRC4, NLRP3, PARP1, RIPK1, TAB2, TAB3, and TNFAIP3 was linked to better survival rate than low expression (Figures 3A–G; Supplementary Figure S1). A network was constructed to investigate the intersection between PRGs and their prognostic significance (Figure 3H).

The consensus clustering algorithm was applied to classify 572 ccRCC patients into three molecular subtypes based on PRG expression. Optimal clustering was achieved with K = 3, dividing the patients into three distinct subgroups: PRGcluster A (294 patients), PRGcluster B (185 patients), and PRGcluster C (93 patients) (Figures 4A,B). PCA confirmed the satisfactory separation between these clusters (Figure 4C). We further analyzed the correlation between PRGs expression and TNM stage, clinical stage, and age. The heatmap showed that most PRGs showed relatively high expression in the PRGcluster B group but comparatively lower in the PRGcluster C group (Figure 4D). Survival analyses were also performed between the three PRG clusters. The survival analysis revealed a significant difference in OS between these three PRGclusters, with PRGclusters A showed the best prognosis for ccRCC (Figure 4E).

PRGcluster A was significantly enriched in oxidative phosphorylation, adipogenesis, Hedgehog signaling, and TGF-β signaling, while PRGcluster B was significantly enriched in the PI3K-AKT-Mtor pathway, interferon-β, IL6-JAK-STA3 signaling, complement, and apoptosis. PRGcluster C was significantly enriched in KRAS pathway signaling and myogenesis (Figures 4F–H). Subsequently, ssGSEA was conducted to assess the differences in immune cell infiltration among the three PRGclusters. PRGcluster B was enriched in activated B cells, activated CD4 T cells, activated CD8 T cells, activated dendritic cells, CD56 bright nature killer cells, gamma delta T cells, B immature cells, dendritic immature cells, MDSCs, macrophages, mast cells, natural killer T cells, plasmacytoid dendritic cells, regulatory T cells, T follicular helper cells, type 1 T helper cells, and type 2 T helper cells. PRGcluster A was enriched in eosinophils and neutrophils. While PRGcluster C was enriched in CD56 dim nature killer cells (Figure 4I).

DEGs among the three PRGclusters were compared in triplicate to explore the underlying biological functions of these PANoptosis patterns. Across these comparisons, we identified 289 intersecting genes (Figure 5A). Subsequently, we conducted GO and KEGG enrichment analyses on these 289 intersected genes. The GO results indicated that these genes were involved in regulating cell-cell adhesion, leukocyte migration, and leukocyte cell-cell adhesion in the biological process (BP) category. The molecular function (MF) category encompassed amide binding, peptide binding, and immune receptor activity. In the cellular component (CC) category, the genes were significantly associated with the external side of the plasma membrane, endocytic vesicle, and membrane microdomain (Figure 5B). The KEGG results revealed that these genes participated in the PI3K-Akt signaling pathway, focal adhesion, tuberculosis pathways, and cell adhesion molecules (Figure 5C).

Univariate Cox regression analysis identified univariately significant genes. ccRCC patients were classified into two gene subtypes (geneCluster A and geneCluster B) based on the highest correlation coefficient K-value (Figures 5D,E). Survival analysis revealed that ccRCC patients in geneCluster A exhibited better survival outcomes compared to those in geneCluster B (Figure 5F). The boxplot showed that the expression of most PRGs, with the exceptions of AIM2, ZBP1, and PYCARD, were higher in geneCluster A than in geneCluster B (Figure 5G). A heatmap illustrated the relationships among gene clusters, PRGclusters, and select clinical characteristics (Figure 5H).

Considering the complexity and heterogeneity of individual PANoptosis patterns, a risk scoring system incorporating prognosis-related PANoptosis DEGs was constructed to predict PANoptosis patterns in individual ccRCC patients. ccRCC patients were randomly divided into training (n = 286) and testing groups (n = 286) using the “caret” R package. LASSO regression was applied to reduce gene overfitting, resulting in the selection of six significant genes (Figures 6A,B). Multivariate COX regression analysis identified the three top-performing genes (WDR72, ANLN, and SLC16A12) The PANoptosis risk scoring formula is as follows: Risk score = Exp (WDR72) × (−0.243) + Exp (ANLN) × (0.132) + Exp (SLC16A12) × (−0.167). All sets consisted of training and testing group files. ccRCC patients were divided into low- and high-risk groups according to their median risk score.

The alluvial diagram depicts changes in the distribution of ccRCC samples (Figure 6C). Figures 6D,E show that PRGcluster C and gene cluster B exhibited higher risk scores. The differential gene analysis of PANoptosis expression between the low- and high-risk groups is presented in Figure 6F. The boxplot results indicated that AIM2, ZBP1, and PYCARD were the only genes expressed in the high-risk group. The Kaplan‒Meier curve demonstrated that the high-risk group had poorer overall survival probabilities across all sets, including the training and testing sets (Figures 6G–I). ROC curves showed that the AUCs at 1-, 3-, and 5-year survival rates were 0.700, 0.692, and 0.722, respectively, in all sets. In the training set and testing set groups, the AUCs consistently exceeded 0.70 (Figures 6J–L). These findings demonstrated satisfactory predictive accuracy for ccRCC prognosis.

Figures 7A–C depict the differences in WDR72, ANLN, and SLC16A12 expression levels across all datasets, including the training and testing sets. The heatmap results indicated that WDR72 and SLC16A12 were highly expressed in the low-risk group across all three datasets. Across all datasets, including the training and testing sets, the risk score was negatively associated with survival times, indicating that ccRCC patients with lower scores experienced better survival outcomes (Figures 7D–J). A nomogram incorporating the PRG risk score, age, and clinical characteristics was constructed to predict 1-, 3-, and 5-year survival outcomes (Figure 7J). The calibration curve indicated a strong correlation between predicted and observed OS (Figure 7K).

The results of correlation analyses between the 22 immune cells and the risk score showed that the risk score was positively correlated with M0 macrophages, mast cells, plasma cells, CD4 memory-activated T cells, follicular helper T cells, and Tregs but negatively correlated with resting dendritic cells, M1 macrophages, resting mast cells, monocytes, resting NK cells, and CD4 memory resting T cells (Supplementary Figure S2). Subsequently, correlations between the three risk genes and immune cell populations were examined, revealing associations between certain genes and specific immune cells (Figure 8A). Results indicated that ANLN expression positively correlated with M2 macrophage infiltration, a hallmark of immune evasion and immunotherapy resistance. WDR72 may facilitate TME remodeling by negatively correlating with activated memory CD4⁺ T cells and positively with resting CD4 memory T cells, influencing immune checkpoint blockade efficacy and possibly leading to therapy resistance via altered T cell dynamics. SLC16A12, involved in metabolite transport, shows positive correlation with M2 macrophages and negative correlation with M0 macrophages, suggesting a role in macrophage phenotype modulation that could enhance TME remodeling toward a less aggressive state, potentially improving prognosis but also affecting sensitivity to certain therapies. Furthermore, ccRCC patients in the high-risk group exhibited significantly higher stromal, immune, and ESTIMATE score than those in the low-risk group (Figure 8B).

Correlation analysis between the risk score and cancer stem cells (CSCs) revealed a significant negative association (Figure 8C). Tumors with low TMB are typically linked to poor prognosis and limited responsiveness to immunotherapy. Our findings indicated that the low-risk group exhibited a lower TMB than the high-risk group, and the risk score was positively correlated with TMB (R = 0.14, p = 0.0076) (Figures 8D,E). Waterfall plots were generated to compare mutation frequencies between the high- and low-risk groups. In the high-risk group, TTN, MUC16, CSMD3, SYNE1, and HMCN1 were the five most frequently mutated genes. In the low-risk group, TTN, MUC16, ARID1A, PCLO, and HMCN1 represented the top five most frequently mutated genes. Overall, the high-risk group exhibited a higher mutation frequency than the low-risk group (Figures 8F,G).

Drug sensitivity was compared between the low- and high-risk groups to evaluate the signature’s predictive efficacy for ccRCC treatments (Supplementary Table S3). Our study found that low-risk patients exhibited greater sensitivity to afatinib, axitinib, bortezomib, and nilotinib. Conversely, high-risk patients displayed significantly lower IC50 values for 5-fluorouracil, cisplatin, irinotecan, and gemcitabine. These findings indicate that the PANoptosis signature could be useful in predicting the drug susceptibility of ccRCC patients (Figures 9A–H).

Additionally, the expression levels of the three marker genes and their ROC curve performance were verified using the GSE40435 dataset. The findings showed increased expression levels for ANLN, while WDR72 and SLC16A12 exhibited decreased expression levels (Figure 10A). As shown in Figure 10B, all three hub genes achieved ROC curve AUC values above 0.80 in the GSE40435 dataset (WDR72, AUC = 0.998; ANLN, AUC = 0.955; SLC16A12, AUC = 0.836). Additionally, the combined AUC for three hub genes was higher than that of any single gene, indicating improved predictive accuracy (AUC = 0.997) (Figure 10C). These findings suggest that the three marker genes may serve as potential diagnostic biomarkers for ccRCC.

Figures 10D–F validated the ccRCC diagnosis. RT‒qPCR was then conducted to compare mRNA expression levels of WDR72, ANLN, and SLC16A12 in adjacent normal and ccRCC tumor tissues. The results revealed ANLN upregulation and SLC16A12 and WDR72 downregulation in ccRCC tumor tissues compared to adjacent normal tissues (Figures 10G–I). We also examined the expression of the three marker genes in normal renal cell (Hk-2) and ccRCC cells (Caki-1, Caki-2, A498, 786-o). The result indicated that compared to HK-2 cells, ANLN expression was upregulated, whereas SLC16A12 and WDR72 were downregulated in ccRCC cell lines (Figures 10J–L).