The RNA-sequencing data (HTSeq-Counts) for ccRCC (n = 539) and normal (n = 72) tissue samples was downloaded from The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/) while the corresponding clinical data of ccRCC samples was downloaded from the cBioportal website (https://www.cbioportal.org/datasets). Raw counts of RNA-sequencing data were transformed into transcripts per million (TPM) values, and were further log2- transformed (log2TPM) for subsequent analyses. Matrix files of gene expression profiles and clinical information of the E-MTAB-1980 cohort was downloaded from the ArrayExpress website (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-1980/).

Normalized transcriptomic and clinical data of CheckMate 025 (CM-025) cohorts of ccRCC patients treated with Nivolumab (anti-PD-1) therapy were obtained from the published article (26). The inclusion criteria for ccRCC samples were: i. Those with sequencing or array data; ii. Those with survival data; iii. Those with a follow-up of ≥ 30 days; and iv. Pathologically confirmed ccRCC cases. Samples from the CM-025 cohorts that did not have information regarding ICB therapeutic responses were excluded. Finally, 539 ccRCC and 72 normal tissue samples were used to identify differentially expressed genes (DEGs). A total of 509 ccRCC samples from TCGA, 101 from E-MTAB-1980 and 172 from CM-025 cohorts were enrolled in this study.

A total of 7399 immune-related genes were obtained from the InnateDB website (https://www.innated b.com/redirect.do?go=resourcesGeneLists). Then, the “edgeR” R package (version 3.30.0) was used to identify DEGs between ccRCC (n = 539) and normal (n = 72) tissue samples in TCGA. The threshold for statistically significant DEGs was set at |log2-fold change (FC)| > 1 and false discovery rate (FDR) < 0.05. Then, immune-related genes were intersected with significant DEGs to obtain immune-related differentially expressed genes (IRDEGs).

A total of 64 immune and non-immune cell types in the TME were downloaded from xCell (http://xcell.ucsf.edu/). xCell is an enrichment algorithm with 6573 gene signatures of 64 immune and non-immune cell types, including epithelial cells, hematopoietic progenitors, extracellular matrix cells, adaptive and innate immune cells. The xCell obtained the full cellular landscape of 64 human cell type from various sources and adopted an integrated approach combining the advantages of gene set enrichment and deconvolution approaches to identify cell types across data sources (28). Cox regression and Kaplan-Meier analyses were performed to identify prognostic immune and stroma cell types by “survival” R package. The “survminer” R package was used to determine the best cutoff.

WGCNA was performed using the “WGCNA” R package to identify IRDEGs associated with CD8⁺ T cells on the TCGA cohort. The expression profile of IRDEGs with log2TPM expression greater than 1 (n = 1339) derived from the above analysis was used as the input file for WGCNA. The WGCNA is a method for evaluating correlation patterns among genes across samples and for visualizing co-expression networks (29). Adjacency matrix was generated by Pearson’s correlation between all pair-wise genes. The soft threshold power of β = 6 was selected to achieve scale-free topology of the adjacency matrix. Then, the adjacency matrix was transformed into topological overlap matrix (TOM). According to TOM‐based dissimilarity measure with minimal module size as 30 and cut height as 0.25, IRDEGs with similar expression patterns were classified into the same gene module by average linkage hierarchical clustering. Then, we evaluated the correlation between module eigengenes (MEs) and clinical traits to identify clinically significant modules.

Kaplan-Meier and univariate Cox regression analyses were performed using the “survival” R package to select genes from the brown module, and genes with p < 0.01 were selected for unsupervised cluster analysis. Then, we performed the consensus non-negative matrix factorization (CNMF) algorithm based on above genes and selected default parameters for molecular classification of the TCGA cohort using the “CancerSubtype” R package (30). The CNMF algorithm is a classical dimension reduction method that extracts essential data from high-dimensional data. Since the CNMF function of the “CancerSubtype” R package is designed for multi-omics data analysis, we randomly divided gene expression profiles into two datasets, as two different omics. Silhouette coefficient, which ranged from -1 to 1, was used to estimate the result of the cluster. Silhouette coefficients near 1 indicate that the sample is distinguished from neighboring clusters and was used to determine the best number of clusters (31). The same classification procedure was performed for the E-MTAB-1980 cohort to validate molecular classification.

Cox regression analysis was performed to assess the prognostic significance of molecular clusters in the TCGA cohort. Moreover, the independence of molecular classification based on clinical characteristics (T stage, M stage, Stage and Grade) was also evaluated. Then, time-dependent receiver operating characteristic (ROC) analysis was compared with area under the curve (AUC), concordance index (C-index) and decision curve analysis (DCA) between molecular clusters and ClearCode34 (32). Brooks et al. developed a model based on 34 gene expressions to classify ccRCC patients and to predict patient survival outcomes (32). The “compare” function in “timeROC” R package and the “cindex.comp” function in “survcomp” R package were used to statistically test the difference in AUC and C-index between different signatures, respectively. Decision curve analysis (DCA) was performed to evaluate the net benefits derived from the use of molecular classification and ccA/ccB subtypes. Subsequently, the associations between molecular clusters, clinical characteristics and ClearCode34 were evaluated. The sankey diagram was used to visualize relationships among molecular clusters, stage, grade, survival status and ClearCode34.

Data for somatic mutations and CNAs were downloaded from the cBioportal website (https://www.cbioportal.org/datasets). Genes with mutation rates greater than 2% were included in this study. Copy-number alterations data were analyzed using the GISTIC 2.0 online version (https://cloud.genepattern.org/gp/pages/index.jsf). Threshold copy numbers at alteration peaks (kirc.all_lesions.conf_90.txt file) were obtained by GISTIC analysis. CNAs with occurrence rates greater than 2% were also included in this study. CNAs were divided into focal and arm levels by GISTIC according to the length of genome mutation fragments. Focal CNAs were defined as shorter than the chromosome arm and arm CNAs were defined as chromosome-arm length or longer (33). The load of copy number loss or gain was defined as the total number of CNAs at focal and arm levels.

Immune-related features, including tumor mutation burden (TMB), neoantigen load, homologous recombination defects (HRD), CTA scores and intratumor heterogeneity (ITH) were also retrieved (34). Immune scores for each TCGA sample were downloaded from the ESTIMATE website (https://bioinformatics.mdanderson.org/estimate/index.html). A total of 178 immunomodulators and chemokines were obtained from the TISIDB website (http://cis.hku.hk/TISIDB/index.php), while 44 immunomodulators and chemokines without corresponding expression profiles were excluded.

The immune suppression score signature (CD274, IDO1, FASLG, CTLA4, PDCD1, LAG3, HAVCR2, PDCD1LG2, IL10, TGFB1, PTGS2) (35, 36) and the signatures of 10 oncogenic pathways were obtained from a previously published study (37). ssGSEA algorithm of “GSVA” R package would rank-normalized gene expression values for each ccRCC sample, and a normalized enrichment scores (NES) of each ccRCC sample was generated using the Empirical Cumulative Distribution Functions of the genes and the remaining genes in each signature (38, 39). The NES represented the degree of absolute enrichment of a gene signature in each ccRCC sample, which can reflect the activity of a signature or a pathway. Subsequently, NES were compared among different clusters.

Annotated gene sets c2.cp.kegg.v7.2.entrez.gmt, c5.go.v7.2.entrez.gmt, c6.all.v7.2.entrez.gmt, c7.all.v7.2.entrez.gmt, c8.all.v7.2.entrez.gmt and h.all.v7.4.entrez.gmt were downloaded from the GSEA website (http://www.gsea-msigdb.org/gsea/downloads.jsp). Comprehensive functional enrichment analyses, including GO enrichment analysis, KEGG enrichment analysis, oncogenic signature, immune signature, cell type signature and “hallmark” signature were performed using “clusterprofiler” R package. Hallmark pools specific well-defined biological states or processes and demonstrates coherent expression, which was extensively utilized in medical studies (40–44). FDR < 0.05 was considered statistically significant.

GSVA, a non-parametric estimation method, is commonly used for variation analyses of biological pathways between different molecular clusters. By using the “c2.cp.kegg.v6.2.symbols” gene set as the reference gene set, we performed GSVA to identify differences in KEGG pathways among different clusters and selected default parameters using the “GSVA” R package. The FDR < 0.05 and t value > 2 were set as the cut‐off criteria.

Based on the top 30 up-regulated DEGs in each cluster, we constructed a simple gene classifier to predict the molecular clusters using nearest template prediction (NTP) function of “CMScaller” R package. The NTP algorithm provides a convenient method for cluster prediction in a testing dataset using a list of signature genes of each cluster, which can be used to single-patient, multi-cluster and cross-platform predictions (45). This method has been successfully applied to clinical classification and prognosis prediction based on gene expression (46, 47). Therefore, this method was suitable for this study to establish a convenient gene classifier for clinical application. Then, the gene classifier was validated in E-MTAB-1980 and CM-025 cohorts and was compared with previous molecular classifications based on the CNMF algorithm. Sankey diagram was used to visualize the relationships between molecular classification and the gene classifier.

R software (https://www.r-project.org/) was used for all computational and statistical analyses. Cluster quality measure- “in-group proportion (IGP)” analysis was performed to evaluate reproducibility of the molecular clusters, and IGP close to 100% meant credible reproducibility of the clusters (48). Expression data of different datasets were normalized by Z-scores prior to IGP analysis. Determination of statistically significant DEGs between two TCGA clusters using the “limma” package were defined as |log2FC| > 1 and FDR < 0.05. The Hazard Ratio (HR) and 95% confidence interval (CI) were generated by Kaplan-Meier and Cox regression analyses. Statistical significance of the comparisons between two groups for continuous variables and categorical variables was estimated by the Student’s T test or Mann-Whitney U test and Chi-square test or Fisher’s exact tests, respectively. Correlations between variables were assessed by Spearman correlation analysis. Two-tailed p ≤ 0.05 was considered statistically significant.

The flowchart of this study was shown in Figure S1A . We compared the RNA expression levels between ccRCC (n = 539) and normal tissues (n = 72) in the TCGA cohort (|log2FC| > 1, FDR < 0.05), and obtained a total of 5640 DEGs (3726 upregulated and 1914 downregulated) ( Supplement Table 1 ). The volcano plot showed the distribution of all DEGs in TCGA ( Figure S1B ). We obtained a total of 1399 IRDEGs with the expression greater than log2TPM 1 through the intersection of DEGs and InnateDB immune-related genes ( Figure S1C and Supplement Table 2 ).

Kaplan-Meier and Cox regression analyses showed that out of 64 immune and non-immune cells, 21 cell types were associated with overall survival of ccRCC. ( Figure S2 and Supplement Table 3 ). Multiple lymphoids were associated with poor prognosis, including CD8⁺ T cells, B-cells, Th1 cells, effector memory CD8⁺ T cell (CD8⁺ Tem), natural killer T cells (NKT), Plasma cells and Th2 cells ( Figure 1A and Figure S2 ). In myeloid cells, eosinophils and conventional dendritic cells (cDC) were associated with good prognosis whereas basophils, mesangial cells, M1 Macrophages and monocytes were associated with poor prognosis. High abundance of CD8⁺ T cells was associated with worse prognosis, which consistent with previous studies (22, 23).

The abundance of CD8⁺ T cells was extracted from the xCell website (http://xcell.ucsf.edu/). To identify key modules that were significantly correlated with the abundance of CD8⁺ T cells, we performed WGCNA on the TCGA-ccRCC dataset after incorporating the 1399 IRDEGs derived from the above analysis ( Figure 1 ). Clustering dendrograms of IRDEGs was shown in Figure 1B and color intensity varies were positively with abundance of CD8⁺ T cells. In terms of survival status, red means dead and white indicates live ( Figure 1B ). Analysis of scale-free fit index and mean connectivity for various soft-thresholding powers was shown in Figures 1C, D , respectively. By setting the cut height = 0.25, the brown and turquoise module eigengenes were combined ( Figures 1E, F ). Figure 1F also show the dendrogram of robust IRDEGs clustered based on a dissimilarity measure. By setting the cut height = 0.25 and β = 6 (scale-free R2 = 0.85), 1399 IRDEGs were divided into eight independent co-expression modules ( Figure 1G ). As shown in the relative diagram of module-trait relationship, the brown module including 648 IRDEGs was most significantly correlated with the abundance of CD8⁺ T cells ( Figures 1G, H and Supplement Table 4 ).

To evaluate the potential biological mechanisms of CD8⁺ T cells, we performed comprehensive functional enrichment analysis of CD8⁺ T cell related genes using the “cluster profiler” R package. Figure 2A and Supplement Table 5 showed that immune-related KEGG pathways, such as chemokine signaling pathway, natural killer cell mediated cytotoxicity, primary immunodeficiency, cytokine-cytokine receptor interaction and hematopoietic cell lineage were enriched. The GO term, including leukocyte proliferation, leukocyte cell-cell adhesion, regulation of T cell activation, positive regulation of cytokine production and T cell activation were mostly associated with CD8⁺ T cell related genes. “Hallmark” gene signatures analysis revealed that CD8⁺ T cells could be involved in multiple biological states or processes, such as complement, IL6 Janus kinase (JAK)-signal transducer and activator of transcription (STAT) signaling (signaling, interferon Gamma (INFγ) responses, inflammatory responses, and allograft rejection. Furthermore, based on the cutoff identified by the “survminer” R package, we divided the TCGA cohorts into two groups, high CD8⁺ T cell- infiltration and low CD8⁺ T cell- infiltration groups. The results from GSVA analysis revealed that the low CD8⁺ T cell- infiltration group was enriched with more metabolism-related processes (e.g., tyrosine metabolism, selenoamino acid metabolism, purine metabolism, butanoate metabolism and propanoate metabolism et al), as well as WNT signaling and transforming growth factor (TGF)-β signaling pathways ( Figure 2B and Supplement Table 6 ). High CD8⁺ T cell- infiltration group were related to many immune-related biological pathways, such as intestinal immune network for IGA production, antigen processing and presentation, T cell receptor signaling pathway, natural killer cell mediated cytotoxicity and cytokine receptor interactions et al, but more related to immunosuppressive pathways, such as primary immunodeficiency, vascular endothelial growth factor (VEGF) signaling, hedgehog signaling, p53 signaling, JAK-STAT signaling, and Toll like receptor signaling ( Figure 2B ). Meanwhile, the abundance of CD8⁺ T cells was positively correlated with immune suppression scores and expression of critical immune checkpoint genes (CTLA4, TIGIT, PDCD1, LAG3 and CD274) ( Figures S3B–F ). However, it was negatively correlated with the expression of CD107a (T cell degranulation related factors) ( Figure S3G ). The TMB or neoantigen load of ccRCC was not correlated with the abundance of CD8⁺ T cells ( Figures S3H, I ). The high CD8⁺ T cell- infiltration group exhibited elevated expression levels of critical immune checkpoint genes ( Figure 2C ). Previous studies reported that the TMB, neoantigen load and immune checkpoint gene expression might not be better indicators for ICB therapeutic efficacy in ccRCC, when compared to other solid tumor types (10, 49). Various subpopulations of CD8⁺ T cells were associated with ICB therapeutic responses in ccRCC (27, 50, 51). However, the above findings suggested that the CD8⁺ T cells and immune checkpoint gene expression, but not mutation or neoantigen loads, might be ideal predictors for ICB therapy. In addition, the high CD8⁺ T cell- infiltration group exhibited lower expression levels of CD107a and elevated immune suppression scores and abundance of immunosuppressive cells (M2 macrophages and Th2 cells) ( Figures 2C–E ), which may explain why ccRCC tumors progress despite high T cell infiltration. Overall, through functional enrichment analysis of CD8⁺ T cell- related genes, new insights can be provided for evaluating the heterogenous phenotypes of CD8⁺ T cells.

To retrieve key genes that were significantly associated with OS, univariate cox regression analysis was used to screen out the prognosis-related brown genes (p_Cox < 0.01). Then, the prognosis-related brown genes obtained from univariate cox regression analysis were further screen by Kaplan-Meier analysis (p_KM < 0.01, Supplement Table 7 ). Using the 84-gene panel, we performed molecular classification using the CNMF algorithm to characterize two distinct ccRCC molecular clusters; C1 cluster (N= 176 samples) and C2 cluster (N= 333 samples) ( Figure 3 ). The C2 cluster exhibited a better OS than the C1 cluster (Log rank p < 0.001, Figure 3A ). In addition, we attempted to divide ccRCC patients of TCGA cohort into 2, 3, 4 and 5 clusters, and found that the silhouette coefficients near 1were the largest when they were divided into two clusters ( Figure 3B and Figures S4A–C ). Therefore, we choose two clusters as the best clusters. Differential expression tested the expression difference between two clusters ( Figure 3C ). The C1 cluster was correlated a higher histological grade, T stage, N status, M stage and stage as shown in the heatmap ( Figure 3D and Table 1 ). The heatmap showed that 84 prognostic genes were correlated with the abundance of CD8⁺ T cells in the TCGA cohort ( Figure S4D ).

To evaluate the independence of molecular classification from clinical characteristics, a total of 475 ccRCC patients ( Supplement Table 8 ) with complete clinical data, including T stage, M stage, stage, and grade, were selected for Cox regression analysis. After the adjustment of clinical characteristics (T stage, M stage, Stage and Grade), molecular classification was found to be an independent prognostic factor for OS outcomes of ccRCC ( Figures 3E, F ). The 5-year AUC and 5-year C-index of molecular classification were somewhat higher than those of ClearCode34 (5-year AUC: 0.65 vs. 0.64; 5-year C-index: 0.64 vs. 0.61), but no statistical significance (p>0.05; Figures S4E, F ). Moreover, the DCA for 5-year OS prediction showed that the net benefit across 0% to 50% threshold probabilities of molecular classification was higher than that of ClearCode34 ( Figure 3G ). Sankey diagram was used to visualize the relationships among molecular clusters, ClearCode34 subtypes and clinical characteristics ( Figure 3H ). Overall, we successfully construct two CD8⁺ T cell-based molecular classifications to predict the OS for ccRCC. The molecular classification in this study was of greater significance for the prediction of clinical prognosis.

The same classification pipeline was performed on the E-MTAB-1980 cohort to validate the stability of molecular classification, and the result of Kaplan-Meier analysis was consistent with that of the TCGA cohort ( Figures S5A–C ). Furthermore, we attempted to divide ccRCC patients of E-MTAB-1980 cohort into 2, 3, 4 and 5 clusters, and found that the silhouette coefficients of two clusters were the largest. Therefore, we choose two clusters as the best clusters. The IGP values were 89.5% and 98.4% for C1 and C2 clusters, respectively, indicating that the two molecular clusters have a high degree of reproducibility in E-MTAB-1980 cohort.

To investigate differences in biological mechanisms between molecular clusters, GSVA was performed in the TCGA cohort. The results showed that both of two molecular clusters were notably enriched in metabolism related pathways ( Figure S6A and Supplement Table 9 ). Strikingly, gene sets related to primary immunodeficiency, nucleotide metabolism, JAK-STAT signaling, Toll like receptor signaling, vascular endothelial growth factor (VEGF) signaling, hedgehog signaling, p53 signaling and cell cycle/apoptosis related pathways were significantly elevated in the C1 cluster. In terms of the C2 cluster, it was dramatically enriched in amino acid metabolism (e.g tyrosine and glutamate metabolism), mammalian target of rapamycin (mTOR) signaling, transforming growth factor (TGF)-β signaling, insulin signaling, peroxisome proliferators-activated receptors (PPARs) signaling and tight junction pathways. It has been reported that these pathways might modify the cancer-related immune microenvironment. For example, previous studies have revealed that VEGF had direct or indirect effects on components of the immune system, including suppressing DC maturation and CD8⁺ T cell proliferation and regulating ICAM1 to suppress NK cell and T cell trafficking, resulting in immunosuppressive outcomes (52). Therefore, it could be deduced that the C1 cluster might develop cancer immune escape through multiple immunosuppressive pathways.

To investigate the mechanisms associated with clinical phenotypic heterogeneity between molecular clusters, we analyzed differences in extrinsic immune escape mechanisms between them. Extrinsic immune escape mechanisms involve four aspects: presence of immunosuppressive cells, lack of immune cells, more fibrosis, and high concentrations of immunosuppressive cytokines, which imply that non-tumor cell components in the TME lead to immune escape of tumor cells (37, 53, 54). Figures 4A, B showed that the C1 cluster had higher immune infiltration and immune suppression scores. Expression level of TMEM173 (STING), a signature of spontaneous initiation of innate immunity (37), was higher in the C1 cluster than in the C2 cluster ( Figure 4C ), which were consistent in the E-MTAB-1980 cohort ( Figure S6B ). The C1 cluster exhibited higher infiltrations of immune cells, such as CD8⁺ T-cells, CD4⁺ T-cells, Th1 cells, dendritic cells (DC), and B cells as well as higher infiltrations of immunosuppressive cells, such as Th2 cells, and stromal cells, including fibroblasts and endothelial cells ( Figure 4D and Supplement Table 10 ). We found that the abundance of Th2 cells of C1 cluster was also higher than that of C2 cluster in the E-MTAB-1980 cohort ( Figure S6C ). Interestingly, the C1 cluster exhibited higher expression levels of chemokines, such as CXCL10, CXCL9 and CCL4, which could attract CD8⁺ T cells and DC (37) ( Supplement Table 11 ). Furthermore, immunosuppressive cytokines were upregulated in the C1 cluster, including IL-10 pathway- and TGF-β signaling pathway-related genes ( Figure 4C ), which were basically consistent in the E-MTAB-1980 cohort ( Figure S6B ) (35). In contrast, the GZMB/CD8A ratio, a signature of anti-tumor immunity efficacy (55), was significantly higher in the C2 cluster ( Figure 4E ). Thorsson et al. identified six immune subtypes (C1–C6) encompassing 30 cancer types based on the pan-cancer immunogenomics analyses and reveal that the immune subtype C3 was enriched in most ccRCC patients, which characterized by immune equilibrium and a better prognosis than the other immune subtypes (34). Notably, they identified two distinct CD8+ T cell-related molecular clusters, which had distinct proportions of the immune subtype C3, of which C1 cluster-immune subtype C3 accounted for 75.0% and C2 cluster-immune subtype C3 accounted for 92.3% (p<0.001) ( Figure 4F ). Consist with previously reports, our results demonstrated that although less immune cell infiltration, anti-tumor immune components and immune escape components might achieve a balance in the C2 cluster. Overall, although the C1 cluster has a higher immune infiltration level, there are multiple extrinsic immune escape mechanisms that result in worse prognostic outcomes. C2 cluster may represent immune equilibrium and a better prognosis.

We further evaluated intrinsic immune escape mechanisms between molecular clusters. It has been reported that there are at least two factors that mediate intrinsic immune escape, including immunomodulators and tumor immunogenicity (37, 54). Potential signatures determining tumor immunogenicity were compared between the two clusters: TMB ( Figure 5A ), neoantigen load ( Figure 5B ), HRD ( Figure 5C ), CTA score ( Figure 5D ), ITH ( Figure 5E ), and MHC-related antigen-presenting capability ( Figure 5G ). The C1 cluster had a higher TMB, neoantigen load, HRD, CTA score, and ITH, indicating that it had more sources of tumor antigens. However, when compared to the C2 cluster, the C1 cluster had low expression levels of MHC I- related antigen-presenting genes ( Figure 5G and Supplement Table 11 ), which contribute to suppressed antigen presentation on tumor surfaces and lower immunogenicity (56, 57). We also found that the expression levels of MHC I- related antigen-presenting genes of C1 cluster was lower than that of C2 cluster in the E-MTAB-1980 cohort ( Figure S6D ). Immunomodulators are involved in other intrinsic immune escape mechanisms and play an essential role in cancer ICB therapy (58). The C1 cluster exhibited elevated expression levels of CD8⁺ T cell exhaustion markers (59), including CTLA4, TIGIT, PDCD1 and LAG3, and suppressed levels of degranulation markers (CD107a) (50) when compared to the C2 cluster ( Figure 5F and Supplement Table 11 ), which were basically consistent in the E-MTAB-1980 cohort ( Figure S6E ). In addition, we confirmed a closely correlation between the expression levels of most immune checkpoint genes and tumor immunogenicity in ccRCC ( Figure S7 ). Overall, the C1 cluster exhibits intrinsic immune escape mechanisms that lead to poor prognostic outcomes.

Tumoral genomic alterations play a pivotal role in cancer initiation, promotion, progression, and therapy. We assessed differences in genomic alterations between two molecular clusters in the TCGA cohort. First, we compared the NES of 10 oncogenic pathways generated by ssGSEA among the molecular clusters. Cell cycle and TP53- related pathways were highly enriched in the C1 cluster while the Hippo, NRF2, PI3K, RAS and TGF-b- related pathways were enriched in the C2 cluster ( Figure S8A and Supplement Table 12 ). Subsequently, we identified different frequencies of somatic mutations between the C1 and C2 clusters. Mutations of BAP1, SETD2, TACC2, MTUS2, POLR2A, ZZEF1, APOB, MEGF10, RELN, NOTCH2, USH2A, and PTEN in the C1 cluster were more frequent than in the C2 cluster ( Figure 6A and Supplement Table 13 ). Loss of PTEN, a tumor suppressor and a member of PI3K-AKT pathway, upregulates the expression of immunosuppressive cytokines and downregulates the expression of IFNG to inhibit T-cell mediated infiltration (56, 60). In addition, PTEN mutation was associated with acquired ICB therapeutic resistance (56, 60). Low expression level of SETD2 is associated with resistance to tyrosine kinase inhibitors (TKIs) in ccRCC patients (61). Our cluster-specific CNAs analysis revealed that chromosomal deletions and amplifications, for example, 9p21.3 deletions were more frequent in the C1 cluster ( Figure 6A and Supplement Table 14 ). Deletion of 9p21.3 is more enriched in tumors with high immune infiltrations and is associated with worse prognostic outcomes in ccRCC after anti-PD-1 therapy (26), indicating that ICB therapeutic efficacy in the C2 cluster might be higher than in the C1 cluster. Furthermore, we found significantly higher gain and loss of CNAs load in the C1 cluster than in the C2 cluster ( Figure S8B ). High CNAs load has been associated with a more aggressive phenotype of ccRCC (62) and ICB therapeutic resistance (63). Based on these above analyses of genomic alterations, we postulate that poor prognostic outcome of the C1 cluster might be correlated with genomic alterations, and that the C1 cluster might be resistant to ICB therapy.

ICB therapy (e.g., anti-PD-1/PD-L1 therapy) has highly enhanced ccRCC treatment. In the CM-025 cohort the same classification procedure as TCGA and E-MTAB-1980 cohorts was performed to validate molecular classification ( Figures 6B–D ). In addition, IGP values were 84.5% and 82.7% for C1 and C2 clusters, respectively, indicating the high reproducibility of molecular classification in the CM-025 cohort. Kaplan-Meier analysis revealed that the C2 cluster exhibited better OS outcomes ( Figure 6B ) and tended to have higher ICB therapeutic response rates than the C1 cluster ( Figure 6E ), indicating that the C2 cluster is more likely to benefit from anti-PD-1 therapy. The genomic differences between the two clusters in the CM-025 cohort were analyzed ( Figure S9 and Supplement Tables 15 , 16 ). We found that the difference results of 9p21.3 deletion, 9q34.3 deletion and 14q32.33 deletion were consistent with the results of TCGA cohort ( Figure 6F and Figure S9 ). The VHL and PBRM1 mutations of the C2 cluster were more frequent than those of the C1 cluster ( Figure 6G and Figure S9 ). It has been proven that the PBRM1 mutations is associated with ICB therapeutic efficacy in ccRCC (26). Overall, we confirm that molecular classification in the CM-025 cohort and the C2 cluster could be used to identify ccRCC patients who are more suitable for anti-PD-1 therapy.

To facilitate the clinical application, the NTP algorithm was used to generate the PIE classifier (PIE classifier: prognosis and ICB therapeutic efficacy of the ccRCC classifier) and to predict molecular classification ( Figure 7A and Supplement Table 17 ) based on the top 30 up-regulated genes in each molecular cluster of TCGA. The PIE classifier was used to predict molecular clusters in TCGA, E-MTAB-1980 and CM-025 cohorts ( Figures 7B–D ). There was a high degree of concordance between prediction results of the PIE classifier and molecular classification, implying that the PIE classifier could simply and reproducibly characterize molecular classification. We also performed Kaplan-Meier analysis of the gene classifier and the results were consistent with the original classification ( Figure S10 ).