The KIRC bulk RNA-seq dataset was available data from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and Clinical Proteomic Tumor Analysis Consortium (CPTAC, https://pdc.cancer.gov/pdc/browse) KIRC cohorts (26, 27). The previously published KIRC scRNA-seq data from normal kidney, tumor-adjacent, and tumor tissues was downloaded from CELLxGENE database (https://cellxgene.cziscience.com/collections/f7cecffa-00b4-4560-a29a-8ad626b8ee08) (28). Only patients diagnosed with clear cell renal cell carcinoma (ccRCC), also referred to as KIRC, were included in this study. All data were retrieved and analyzed using established bioinformatics pipelines; no novel algorithms or computational methods were developed in this study.

The scRNA-seq data were analyzed using the Seurat package in R (28). Cells were filtered based on quality control metrics, including the number of detected genes, total gene expression counts and the percentage of mitochondrial gene expression. Specifically, cells with less than 200 expressed genes, counts less than 1000 or more than 10000, or percentage of mitochondrial gene expression larger than 20%, were removed. Data were normalized using log normalization, and variable features were identified using the FindVariableFeatures function. Principal component analysis (PCA) was performed for dimensionality reduction[9], and the first 20 principal components were used for Uniform Manifold Approximation and Projection (UMAP) visualization. Cell clusters were identified using the FindClusters function with a resolution parameter of 0.5. Cell types were identified based on canonical gene markers as follows: CD3D, CD3E, CD4, LTB, IL7R, and MAL for CD4+ T cells; CD3D, CD8A, CD8B, GZMK, and NKG7 for CD8+ T cells; GNLY, GZMB, NKG7, KLRD1, KLRF1, and CD247 for natural killer (NK) cells; MS4A1, CD79A, CD79B, and CD19 for B cells; IGKC, IGHG1, JCHAIN, MZB1, and TNFRSF17 for plasma cells; LYZ, MRC1, CD68, CD163, and C1QC for macrophages; PLD4, ITM2C, LILRA4, IRF7, and IRF8 for plasmacytoid dendritic cells (pDCs); CD1C, CD86, HLA-DRA, and HLA-DRB5 for conventional dendritic cells (cDCs); TPSAB1 and TPSB2 for mast cells; PLVAP and TIMP3 for endothelial cells; ACTA2 and TAGLN for fibroblasts; GPX3 and GATM for proximal tubule (PT) epithelial cells; DEFB1 and UMOD for non-PT epithelial cells; and CA9 for renal cell carcinoma (RCC) cells.

The proportions of different cell types in normal kidney, tumor-adjacent, and tumor tissues were visualized using stacked bar plots. The expression levels of cell-type-specific markers were displayed using dot plots, where the size of each dot represented the percentage of cells expressing the marker, and the color represented the average expression level.

To evaluate the overrepresentation of specific macrophage subtypes across different tissue or clinical conditions, we calculated the odds ratio (OR) and significance levels using Fisher’s exact test based on a 2×2 contingency table framework (29). For each combination of cell type and condition, we computed the number of cells belonging to a specific cell type within a given condition, compared to the number of other cells across all other combinations. The following categories were considered: (1) target cell type in target condition, (2) target cell type in all other conditions, (3) all other cell types in target condition, and (4) all other cell types in other conditions. The contingency table was analyzed using the fisher.test() function in R. P-values were adjusted for multiple comparisons using the Benjamini-Hochberg (BH) method. Full code is available in github (xxx).

Macrophage subpopulations were identified using unsupervised clustering of scRNA-seq data. Five transcriptionally distinct subsets were delineated and visualized using UMAP. Subtype annotation was performed based on enriched marker genes, integrating both functional characteristics and previously reported molecular signatures (29, 30). Based on this comprehensive analysis, we identified and annotated five distinct macrophage subsets: (1) ALOX5AP+ lipid-associated macrophages (ALOX5AP+LA-Mac), characterized by high expression of ALOX5AP, APOE, and APOC1, indicative of lipid metabolism and arachidonic acid processing; (2) HERPUD1+LA-Mac, marked by HERPUD1, APOE, and APOC1, suggesting endoplasmic reticulum stress response within lipid-metabolic macrophages; (3) PRDX1+LA-Mac, defined by PRDX1, APOE, and APOC1, representing antioxidant-responsive lipid-associated macrophages; (4) FCN1+ inflammatory macrophages (FCN1+Inflam-Mac), characterized by FCN1, IL1B, and ISG15 expression, indicative of pro-inflammatory activation; and (5) oxidative phosphorylation macrophages (OxP-Mac), marked by mitochondrial genes ATP5ME, MT-ND3, and MT-ND2, reflecting enhanced metabolic activity. Heatmaps were generated to display the z-score normalized expression of subtype-specific markers. The proportions of each macrophage subpopulation in different tissue contexts were analyzed. The functional polarization analysis was performed to assess M1/M2 polarization scores by the AddModuleScore function from Seurat package. The M1/M2 gene signatures were obtained from previously published research (31). Pseudotemporal trajectory analysis was used to reveal the dynamics of macrophage subpopulation differentiation by R package Monocle3 (31). The OxP-Mac was set to be starting point in trajectory analysis. The trajectory regulation genes were determined by the graph_test function from Monocle3 package with the threshold of adjusted p value < 1e-26 and Morans’ I > 0.25. Pathway enrichment analysis and gene set enrichment analysis (GSEA) were performed to identify key biological processes associated with each subpopulation.

The hdWGCNA package was used to identify gene modules from macrophages in scRNA-seq data (32). A soft threshold of 10 was chosen to construct the network, achieving the best scale-free topology. Key modules were defined as those specifically enriched in macrophages and highly expressed in tumor tissues. Volcano plots were generated to highlight the most significantly differentially expressed modules between tumor and control samples (including both tumor-adjacent and healthy tissues). Core genes were identified based on high module membership (kME) and differential expression between macrophages from tumor and control samples.

To assess the prognostic value of macrophage-related core genes, we constructed 20 prognostic models by combining different machine learning algorithms and parameter settings. These included Random Survival Forest (RSF), Elastic Net (Enet), Stepwise Cox, CoxBoost, Partial Least Squares Regression for Cox (plsRcox), Supervised Principal Components (SuperPC), Generalized Boosted Regression Modeling (GBM), and Survival Support Vector Machine (survival-SVM). Genes that overlapped between macrophage-related core genes and those associated with poor prognosis in the TCGA cohort were selected for model training. The TCGA-KIRC cohort was used as the training set, while the CPTAC-KIRC cohort served as the validation set. Prior to machine learning analysis, the gene expression data were log2-transformed and standardized (zero mean, unit variance) to minimize heterogeneity and ensure comparability across features. The RSF algorithm was selected for its optimal performance. Survival analysis, including Cox proportional hazards and Kaplan–Meier (KM) analyses, was conducted using the survival R package. For KM analysis, the optimal cutoff value was determined using the surv_cutpoint function from the survminer package, which stratified samples into high- and low-risk groups. Kaplan–Meier survival curves were visualized using the ggsurvplot function from survminer to compare survival outcomes between the high-risk and low-risk patients, as defined by the RSF model. A time-dependent receiver operating characteristic (ROC) curve was used to evaluate the prognostic performance of the RSF model. Feature importance analysis was performed to identify key prognostic genes.

To quantify the proportions of macrophage subpopulations in bulk transcriptomic datasets, we applied the MuSiC (Multi-subject Single Cell deconvolution) algorithm, which uses cross-subject single-cell RNA-seq references to infer cellular composition in bulk samples (33). We used macrophage subpopulations of single-cell RNA-seq data as the reference. MuSiC deconvolution was performed on the TCGA and CPTAC KIRC bulk RNA-seq datasets. Prior to deconvolution, bulk expression matrices were log-normalized, and only genes present in both bulk and single-cell datasets were retained. We executed the analysis with the music_prop function from the MuSiC R package (v1.0.0) using default settings. The output provided relative abundance estimates of each macrophage subtype in every bulk sample.

The mutation profile, including missense, silent, nonsense, frameshift/in-frame insertions and deletions, and uninterrupted mutations, was analyzed, excluding synonymous mutations. Tumor mutation burden (TMB) was calculated based on the total number of somatic mutations. Pathway enrichment analysis, focusing on known oncogenic or user-defined pathways, was performed using the pathways function. Differences in mutational profiles between KIRC risk groups were assessed using the mafCompare function. All analyses in this section were conducted using the maftools R package (34).

The tumor microenvironment was estimated using the xCell R package (35). Differences in macrophage polarization between high-risk and low-risk groups were assessed using the Wilcoxon rank-sum test. For Kaplan–Meier (KM) analysis of macrophage polarization scores, the optimal cutoff value was determined using the surv_cutpoint function from the survminer package.

Differentially expressed genes (DEGs) between high-risk and low-risk patients were identified using R package limma with a threshold of |log2 fold change| > 1 and adjusted p-value < 0.05 (36). Gene Ontology (GO) enrichment analysis was performed on both upregulated and downregulated genes by R package clusterProfiler (37). Single-sample gene set enrichment analysis (ssGSEA) of the HALLMARK gene set was conducted to assess pathway activity in high-risk and low-risk tumors (38). The HALLMARK gene set was obtained from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). All statistical analyses were conducted using R version 4.4.3. A p-value < 0.05 and adjusted p-value < 0.05 were considered statistically significant.

We first employed a KIRC scRNA-seq to explore the cellular heterogeneity across the tissue regions. After the quality control, normalization and cell annotation, a total of 187,290 cells were divided into 14 celltypes, including NK cells, B cells, plasma cells, classical dendritic cells (cDCs), plasma dendritic cells (pDCs), CD4+ T cells, CD8+ T cells, macrophages (Macs), mast cells, endothelial cells (ECs), fibroblasts (FIBs), non-proximal tubule epithelial cells (nonPT-Epi), proximal tubule epithelial cells (PT-Epi), and renal cell carcinoma (RCC) cells (Figure 1A). The tissue origins of cells were classified as normal kidney, tumor-adjacent, or tumor tissue (Figure 1B). The cell types were determined by the canonical marker genes, including CD3D, CD3E, CD4, LTB and IL7R for CD4+ T cells; CD3D, CD8A, CD8B, GZMK and NKG7 for CD8+ T cells; GNLY, NKG7, KLRD1 and CD247 for NK cells; MS4A1, CD79A, CD79B and CD19 for B cells; IGKC, IGHG1, JCHAIN and MZB1 for plasma cells; LYZ, MRC1, CD68, CD163 and C1QC for Macs; PLD4, LILRA4, IRF7, IRF8 and SERPINF1 for pDCs; CD86, HLA-DRA, HLA-DRB5 and HLA-DRB1 for cDCs; TPSAB1 and TPSB2 for mast cells; PLVAP and TIMP3 for ECs; DEFB1 and UMOD for nonPT-Epi; GPX3 and GATM for proximal tubule cells; ACTA2 and TAGLN for FIBs; and CA9 for RCCs (Figure 1C). The stacked bar chart displays the proportions of different cell types in normal kidney, tumor-adjacent, and tumor tissues, while the proportions of CD8+T cells and cDCs are significantly increased in tumor (Figure 1D). Moreover, the OR analysis quantifies the enrichment of different cell types in tumor tissues relative to normal tissues (Figure 1E). CD8+T cells, cDCs and RCCs are significantly enriched in tumor tissues (OR > 1), while the ECs and epithelial cells demonstrated the opposite pattern. Interestingly, Macs demonstrated the limited difference in OR (Figure 1E), suggesting that the functional alternation might be more important in KIRC progression. These findings revealed the significant heterogeneity of immune and stromal cells in KIRC.

We next sub clustered the Macs to explore the alternative activation in KIRC. According to the enriched marker genes integrated both functional characteristics and distinctive molecular signatures, macrophages (Mac) were identified five transcriptionally distinct subsets: (1) ALOX5AP+ lipid-associated macrophages (ALOX5AP+LA-Mac), characterized by high expression of ALOX5AP, APOE, and APOC1, indicative of lipid metabolism and arachidonic acid processing; (2) HERPUD1+LA-Mac, marked by HERPUD1, APOE, and APOC1, suggesting endoplasmic reticulum stress response within lipid-metabolic macrophages; (3) PRDX1+LA-Mac, defined by PRDX1, APOE, and APOC1, representing antioxidant-responsive lipid-associated macrophages; (4) FCN1+ inflammatory macrophages (FCN1+Inflam-Mac), characterized by FCN1, IL1B, and ISG15 expression, indicative of pro-inflammatory activation; and (5) oxidative phosphorylation macrophages (OxP-Mac), marked by mitochondrial genes ATP5ME, MT-ND3, and MT-ND2, reflecting enhanced metabolic activity (Figures 2A, B). Proportions of macrophage subpopulations shifted markedly across tissue contexts: OxP-Mac was enriched in normal tissues (up to 60%), while the other subtypes dominated tumor regions (Figure 2B). Subtype-specific markers validated by the z-score-normalized expression dotplot (Figure 2B). Functional polarization analysis showed these macrophage sub types were not canonical M1/M2 features, while tumor-enriched PRDX1+ and HERPUD1+LA-Mac exhibited mixed M1/M2 features (Figure 2C). Although the Mac substypes exhibited mix M1/M2 signatures, normal kidney enriched OxP-Mac demonstrated the lowest M1 as well as M2 scores compared to the other subtypes (Figure 2C), indicating a quiescent M0 state. Additionally, KEGG GSEA results highlighted an alternative functional activation among the Mac subtypes (Figure 2D). Specifically, HERPUD+ and PRDX1+LA-Macs demonstrated a similar phenotype characterized by hyper-inflammation and immune response while normal kidney enriched OxP-Macs were activated in ribosome and oxidative phyphorylation (Figure 2D).

We further employted pseudotime trajectory analysis to identify the development trajectory as well as regulation genes. Pseudotime trajectory analysis positioned OxP-Mac as an early-state population branching into tumor-specific subsets (ALOX5AP+, HERPUD1+ and PRDX1+LA-Mac), suggesting differentiation plasticity (Figure 2E). The regulation genes were divided into 5 clusters, among which cluster1 and cluster2 promoted this differentiation while cluster4 and cluster5 exhibited the opposite (Figure 2F). GO enrichment revealed that cluster1 regulation genes were functionally enriched in antigen presentation (MHC class II assembly, peptide antigen processing) (Figure 2G), while cluster2 genes were involved in leukocyte apoptosis and immune suppression (Figure 2H). These findings underscore functional diversification of tumor-associated macrophages and their potential roles in KIRC progression.

To delineate macrophage-associated transcriptional networks, we applied high-dimensional weighted gene co-expression network analysis (hdWGCNA) to single-cell RNA-seq data. The optimal power of 10 was employed to construct a scale-free network which identified 10 gene co-expression modules (Figures 3A, B). The macrophage module 4 (Mac-M4) and Mac-M8 demonstrated a high correlation among these modules (Figure 3C). As Figure 3D shown, the top 5 genes were labelled in the UMAP. Among the 10 modules, Mac-M2 was highly enriched in macrophages compared to other cell types (p < 0.0001, Figures 3E, F), highlighting the genes of M2 were specifically expressed in macrophages. Moreover, M2 demonstrated the most differentially expressed modules with elevated log2 fold changes (log2FC > 3) and statistical significance (p < 1 × 10^-175) (Figure 3G). Thus, Mac-M2 was selected for advanced analysis. To further clarify the hub genes in the Mac-M2 cluster, we employed two criteria: module membership (kME) and differential expression levels between tumor and normal macrophages (Figure 3H). Interestingly, these two values showed a strong correlation (r = 0.53, p < 2 × 10^-16), suggesting that the associated genes may play a key role in macrophage differentiation (Figure 3H). Genes with kME > 0.5 and log2FC > 0.3 were selected for further analysis.

We constructed a prognostic model based on TCGA-KIRC data and validated it in the CPTAC-KIRC cohort to assess the predictive value of macrophage-associated hub genes (Figures 4A–C). Shared genes from the Mac-M2 module hub and TCGA-identified survival-related genes were used to build survival models (Figures 4A–C). Seven genes—GPX1, LAIR1, FCGR1A, CD14, SERPING1, RGS1, and APOC1—were selected during model training (Figures 4B, C). The TCGA-KIRC transcriptome served as the training set, while CPTAC-KIRC data were reserved for validation. Twenty machine learning models, generated through combinations of algorithms and parameter settings, were evaluated (Figures 4A, D). Among them, the Random Survival Forest (RSF) model demonstrated the best performance in the training cohort (C-index = 0.94), with 1-, 3-, and 5-year AUCs of 0.97, 0.98, and 0.98, respectively (Figures 4D–F). Although the CPTAC-KIRC validation AUCs were lower, they remained clinically meaningful at approximately 0.70 for 1-, 3-, and 5-year predictions (Figure 4G). This performance reduction may reflect heterogeneity between cohorts, particularly in median overall survival (TCGA: 1,188 days vs. CPTAC: 781 days, Supplementary Figure S1A) and survival status distributions (32.9% vs. 19.4%, Supplementary Figure S1B), although these differences did not reach statistical significance (log-rank p = 0.30, Supplementary Figure S1C). Patients were stratified into high- and low-risk groups based on the RSF risk scores; both groups showed significantly different survival outcomes in TCGA and CPTAC cohorts (log-rank p < 0.0001), confirming the model’s prognostic utility (Figures 4H, I).

To bridge our single-cell macrophage subtyping analysis with the prognostic model, we first examined the expression specificity of RSF feature genes across different macrophage subtypes identified in our scRNA-seq analysis. Among the genes contribution to RSF survival mode, FCGR1A had the highest variable importance, followed by CD14 and LAIR1 (Figure 5A). Notably, RSF feature genes were predominantly enriched in the ALOX5AP+, HERPUD1+, and PRDX1+LA-Mac subtypes (Figure 5B). Similarly, the Mac-M2 score exhibited the same pattern (Figure 5C), indicating that LA-Macs represent the cellular origin of Mac-M2 hub genes.

We then integrated single-cell–derived macrophage subtypes with bulk transcriptomic data from the TCGA and CPTAC cohorts using MuSiC deconvolution and ssGSEA scoring. MuSiC deconvolution revealed that PRDX1+LA-Macs were significantly enriched in high-risk patients compared to low-risk patients in both cohorts, whereas other macrophage subtypes displayed variable enrichment (Figure 5D; Supplementary Figures S2A, B). Consistently, ssGSEA of the top 30 subtype marker genes confirmed higher PRDX1+LA-Mac scores in high-risk patients (Figure 5E; Supplementary Figures S2C, D), and ssGSEA of the Mac-M2 hub genes showed the same trend (Figure 5F; Supplementary Figures S2E, F). Finally, individual RSF feature genes were significantly upregulated in high-risk patients in both cohorts—except for SERPING1 in CPTAC (Figures 5G–H). Collectively, these results implicate PRDX1+LA-Macs and Mac-M2 hub genes in driving unfavorable prognosis in KIRC patients.

To further characterize molecular distinctions between high- and low-risk patient groups, we analyzed the mutation profile of the TCGA-KIRC cohort. In the high-risk group, 78.5% (84/107) of samples harbored at least one driver mutation; the most frequently mutated genes were VHL (46%), PBRM1 (38%), TTN (19%), and BAP1 (17%). Missense mutations and frameshift deletions predominated (Figure 6A). In the low-risk group, 80.8% (210/260) of samples carried driver mutations, with VHL (47%) and PBRM1 (40%) slightly more frequent; frameshift insertions were also more common (Figure 6B). Pathway enrichment analysis showed that mutations in the RTK-RAS (34.6%), Hippo (28.0%), and NOTCH (27.1%) pathways were significantly enriched in the high-risk group, whereas RTK-RAS (33.8%), NOTCH (23.5%), and PI3K (21.9%) pathways dominated in the low-risk group (Figures 6C, D). High-risk patients exhibited a significantly elevated tumor mutation burden (TMB) compared to low-risk patients (Figure 6E). Moreover, mutation frequencies of ZNF804B and FAM178A (5–6%) were significantly higher in high-risk patients, with odds ratios >3 (Figures 6F, G). These findings suggest that elevated TMB and specific gene mutations underlie the poorer prognosis observed in high-risk KIRC patients.

We further performed the analysis of immune cell infiltration in TCGA-KIRC cohorts to reveal critical associations between macrophage polarization and clinical prognosis (Figure 7A). Interestingly, the patients of high-risk group demonstrated a higher immune infiltration and lower stromal proportion compared to low-risk group (Figure 7B), indicative of a “hot” tumor microenvironment. Moreover, the patients in high-risk group had a higher score of overall macrophages and M1 subtypes while M2 subtype macrophages were not statistically significant (Figure 7C). As for M1/M2 ratio, high-risk group as well had a higher ratio compared to low-risk (Figure 7C). It was notable that Th2 cell score was dramatically enriched in high-risk group around 4.3 times (Supplementary Table S1). Patients with high M2 macrophage infiltration exhibited significantly better overall survival compared to those with low M2 levels (log-rank p = 0.018, Figure 7D). Conversely, a high M1 infiltration and M1/M2 macrophage ratio was strongly linked to improved survival (log-rank p < 0.01, Figure 7D), suggesting a harmful role of M1-polarized macrophages. Notably, total macrophage abundance alone showed no significant survival correlation (p=0.064, Figure 7D).

We further assess the clinical relevance of the newly identified macrophage subtype, PRDX1⁺LA-Mac. Patients with higher PRDX1⁺LA-Mac abundance exhibited significantly worse overall survival (log-rank p < 0.0001), with a more pronounced prognostic separation than observed for M1, M2, or total macrophages (Figure 7E). Spearman correlation analysis indicated only weak positive associations between PRDX1⁺LA-Mac and classical macrophage subsets (ρ = 0.11–0.18, all p < 0.05), suggesting that PRDX1⁺LA-Mac represents a transcriptionally and functionally distinct population (Figure 7F). Univariate Cox regression confirmed the strong prognostic value of PRDX1⁺LA-Mac, showing the highest hazard ratio and strongest statistical significance among all macrophage-related features (Figure 7G). Collectively, these results establish PRDX1⁺LA-Mac as a unique and highly informative prognostic macrophage subtype in KIRC.

Transcriptomic analysis of TCGA-KIRC data identified 394 differentially expressed genes (DEGs) between high- and low-risk patients (|log2FC| > 1, adj. p < 0.05). Key upregulated genes in high-risk tumors included SAA1 (log2FC = 3.2, p.adj = 8.4e-17), PAEP (log2FC = 2.5, p.adj = 2.3e-13), and SLPI (log2FC= 2.3, p.adj = 7.3e-14), associated with inflammatory responses and immune evasion (Figure 8A). Downregulated genes such as NPR3 (log2FC = −1.6, p.adj= 3.0e-17), G6PC1 (log2FC = −2.3, p.adj= 2.7e-12) and SLC16A12 (log2FC = −1.77, p.adj=2.9e-17) were linked to metabolic dysregulation (Figure 8A). Gene Ontology (GO) enrichment revealed distinct functional profiles: upregulated DEGs were enriched in immune-related processes (acute inflammatory response, neutrophil chemotaxis, acute-phase response), lipid metabolism (high-density lipoprotein particle) and cellular matrix components (collagen-containing ECM), while downregulated DEGs implicated metabolic pathways (urate metabolism, fatty acid transport) and epithelial polarization (apical plasma membrane, brush border) (Figures 8B, C). Subsequent evaluation of hallmark signaling heterogeneity between high- and low-risk groups by ssGSEA revealed higher levels of proliferation, inflammatory response, and epithelial-mesenchymal transition in high-risk patients (Figure 8D). These findings highlight a pro-tumorigenic microenvironment in high-risk KIRC, characterized by immune activation, metabolic reprogramming, and stromal remodeling.