To determine the correlation between the immune-related risk score (IRS) and efficacy of cancer therapy, five immunotherapeutic cohorts with available RNA-seq data and clinicopathological parameters were included in this study: (1) IMmotion151 cohort, advanced patients with RCC treated with atezolizumab plus bevacizumab (22), (2) JAVELIN Renal 101 trial, advanced patients with RCC treated with the combination of avelumab (anti-PD-L1) + axitinib (TKI targeting VEGF receptors) vs. sunitinib (multitarget TKI) (23), (3) E-MTAB-3218 dataset, patients with mRCC treated with nivolumab (anti-PD-1 ICI) (24), (4) IMvigor210 cohort, patients with advanced urothelial cancer treated with atezolizumab (25), and (5) GSE78220 cohort, patients with metastatic melanoma treated with pembrolizumab (anti-PD-1 antibody) (26).

All data from the IMmotion151 cohort is deposited in the European Genome-Phenome Archive under the accession number EGAS00001004353, and we obtained it according to Hoffmann-La Roche policy. Clinical response information and normalized RNA-seq data were acquired from the supplementary material of Choueiri et al. (23). Original data of E-MTAB-3218 dataset were downloaded from https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-3218?accession=E-MTAB-3218# (24). Additionally, RNA-seq and clinical data for the IMvigor210 cohort were obtained from http://research-pub.gene.com/IMvigor210CoreBiologies. Raw data were normalized using the “DEseq2” R package and further transformed into TPM values. RNA-seq data (FPKM-normalized) and the clinical phenotypes of 28 melanoma patients in the GSE78220 cohort were downloaded from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE78220. Additionally, well-recognized immune-related genes were downloaded from http://www.gsea-msigdb.org/gsea/msigdb/index.jsp (27).

Above datasets merely contain anonymized and de-identified patient information. Secondary analysis of de-identified data was confirmed exempt from review by the medical ethics committee of Hunan Cancer Hospital and the Affiliated Cancer Hospital of Xiangya School of Medicine as it was classified as negligible risk research. Thus, our study was exempt from the ethical review or the patient consent.

WGCNA is a widespread systematic algorithm used to generate gene modules with similar expression patterns and determine the correlations between modules and clinical traits (28). In this study, we screened the immune-related gene expression profiles in the IMmotion151 cohort and further identified a correlation network including significant clinical characteristics and genes using the “WGCNA” R package. We also developed an adjacency matrix to characterize the correlation strength between the nodes, which was further changed to a topological overlap matrix. Subsequently, modules containing more than 30 genes were identified via hierarchical clustering. To compare the co-expression levels, modules were clustered based on their correlation with module eigengenes (MEs). When the correlation of MEs > 0.80, module merging was performed, indicating that the expression profiles of the modules were similar (29). Pearson’s correlation coefficient was used to assess the correlations between the modules and various clinicopathological parameters. Finally, gene significance (GS) and module membership (MM) were used to quantify the relationships between the genes and the clinicopathological characteristics in the module. Hub genes were considered as those with MM > 0.8 and GS > 0.2 (29).

Based on the results of WGCNA analysis, 269 immunotherapy-associated genes in the turquoise module were selected for further analysis. Meanwhile, DEGs between the patients with complete response (CR) and those with progressive disease (PD) in the IMmotion151 cohort were identified using the “ggplot2” R package (29). Nine immune-related genes significantly affecting the patient responsiveness to immunotherapy were ultimately identified using the intersection of the above two kinds of genes and used to construct an IRS model.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted using the “clusterProfiler” R package to identify the DEG-associated signaling pathways and biological processes (30, 31). Pathways with a nominal P < 0.05 and false-discovery rate (FDR) < 0.05 were considered to be statistically significant.

Given the individual heterogeneity and intricacy of the clinical outcomes of advanced RCC cases treated with ICB combined with anti-angiogenic drugs, we formulated a scoring system, termed as the IRS model, using the ssGSEA algorithm on the basis of the mRNA expression levels of the identified nine immunotherapy-related genes in a single sample to quantify the prognostic level of each patients with mRCC for in-depth analysis.

Optimal cut-off point identified by the “surv-cutpoint” function of the “survminer” R package stratified all advanced RCC cases in the IMmotion151 cohort into high- or low-risk subgroups. In this approach, different values are grouped as cut-off values for statistical testing, and the result with the lowest P value is considered as the optimal cut-off point that corresponds to the most significant association with the clinical outcome.

A heat map was constructed to visualize the IRS distribution and clinicopathological parameters. Survival analysis between high- and low-risk group were conducted using Kaplan–Meier curves with log-rank test and the “survival” R package (32). Hazard ratios (HRs) and the corresponding 95% confidence intervals (CIs) were estimated.

Area under the curve (AUC) values of the receiver operating characteristic (ROC) curves established using the “survival ROC” R package were used to evaluate the predictive efficiency of the IRS model (33). In addition, to test the robustness of our IRS model, we verified its predictive capability using other external independent datasets: JAVELIN Renal 101, E-MTAB-3218, IMvigor210, and GSE78220 cohorts.

An identical median value based on the IMmotion151 cohort was also applied to the validation groups, effectively categorizing all patients into high- and low-risk groups.

We further determined the correlations between IRS and subsequent biological processes. Mariathasan et al. (25) curated multiple gene sets associated with specific biological pathways, including (1) CD8⁺ T-effector signature (34), (2) antigen processing machinery (35), (3) epithelial–mesenchymal transition (EMT) biomarkers (26, 36, 37), (4) immune checkpoints (25).

RNA-seq data of approximately 1000 tumor cell lines, AUC values for evaluating the efficacy of antineoplastic drugs in tumor cell lines, and targets or pathways of drugs were downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) (38). Spearman correlation coefficient was used to evaluate the correlation between drug sensitivity and the IRS model, and |Rs| > 0.2 and P-value < 0.05 were considered to be significant.

All immune-related genes from GSEA are listed in Supplementary Table 1 . After we intersected the above immune-related genes with the RNA-seq data of 407 patients with advanced RCC from the IMmotion151 cohort, 942 immune-related genes were identified and further subjected to WGCNA analysis ( Supplementary Table 2 ). In line with the standard scale-free network distribution, the soft threshold power value was determined to be three ( Supplementary Figure 1 ). Based on this dissimilarity, a dendrogram of all the gene clusters was formulated, which displayed 10 different modules ( Figure 1A ). Correlations among all clinical modules are illustrated in Figures 1B, C . We further assessed the correlations between MEs and clinical traits, including PD-L1 IHC, MSKCC risk score, PFS, objective response, metastatic status, and sarcomatoid histology. Turquoise module was most significantly associated with the objective response of patients with mRCC to atezolizumab plus bevacizumab (r = 0.45, P < 0.0001) ( Figure 1D ), indicating that genes in the turquoise module potentially exert a crucial effect on the clinical outcome of atezolizumab plus bevacizumab interventions. This turquoise module was further analyzed, and the genes in this module were found to be significantly correlated to the efficacy of atezolizumab plus bevacizumab therapy (r = 0.76, P < 0.0001) ( Figure 1E ).

We analyzed the gene expression profiles of 23 responders (patients with CR) and 75 non-responders (patients PD) in the IMmotion151 cohort, and identified 72 DEGs associated with the effects of atezolizumab plus bevacizumab (│log₂FC│ > 1, P < 0.05, FDR < 0.05) ( Supplementary Table 3 ). Moreover, we cross-referenced the above 72 DEGs and all genes in the turquoise module to select a total of nine DEGs, including seven downregulated DEGs in responders (serine protease inhibitor Kazal type 5 [SPINK5], semaphorin 3E [SEMA3E], roundabout guidance receptor 2 [ROBO2], bone morphogenetic protein 5 [BMP5], orosomucoid 1 [ORM1], C-reactive protein [CRP], and cathepsin E [CTSE]; log₂FC < 1; P < 0.05) and two up-regulated DEGs in responders (promelanin-concentrating hormone (PMCH) and C-C motif chemokine ligand 3-like 1 [CCL3L1]; log₂FC > 1; P < 0.05) ( Figures 1F, G ).

We further identified the mRNA expression profiles of the above nine genes in mRCC and found that, compared with the responders to atezolizumab plus bevacizumab therapy in the IMmotion151 cohort, the expression levels of seven genes (SPINK5, SEMA3E, ROBO2, BMP5, ORM1, CRP, and CTSE) were significantly increased and those of two genes (PMCH and CCL3L1) were decreased in the non-responders (P < 0.05) ( Figure 2A ). GO analysis of the nine DEGs linked them to neutrophil-mediated immunity, leukocyte migration, acute inflammatory response, humoral immune response, and cell chemotaxis ( Figure 2B ; Supplementary Table 4 ), most of which were associated with the modulation of immunity and immunotherapy. Similarly, KEGG pathway analysis revealed that these DEGs were correlated with cytokine–cytokine receptor interactions, complement and coagulation cascades, antigen processing and presentation, and neutrophil extracellular trap formation ( Figure 2C ; Supplementary Table 5 ), indicating their significance and conferring the basis to investigate a potential association between these genes and immunophenotypes.

Based on the strength of the optimal cut-off point (–0.46), 407 individuals with mRCC were stratified into high- and low-risk groups (high: IRS > –0.46 and low: IRS < –0.46) ( Supplementary Table 6 ). To investigate the potential mechanism of the impact of IRS on atezolizumab plus bevacizumab therapy for mRCC, a combined heat map was constructed to visually demonstrate the correlations between IRS and multiple clinicopathological parameters, including PD-L1 IHC, MSKCC risk score, objective response, metastatic status, and sarcomatoid histology, in the high and low IRS groups. We compared the high and low IRS groups in the IMmotion151 cohort and found that more patients in the IRS-low group had positive PD-L1 IHC results than those in the IRS-high group. Metastatic tumors were primarily distributed in the IRS-high cluster, whereas the proportion of sarcomatoid histological subtypes in the IRS-low group was greater than that in the IRS-high group. Additionally, patients with low IRS were primarily characterized by a higher expression of CD8⁺ T-effectors, antigen-processing machinery, and immune checkpoint signatures in the IRS-low group than in the IRS-high group. Conversely, patients in the IRS-high group showed relatively high expression levels of EMT-associated genes ( Figure 3A ). We further assessed the clinical outcomes of patients treated with atezolizumab plus bevacizumab in the IRS-high and -low groups. Survival analysis revealed that short PFS in patients with high IRS (HR = 1.91; 95% CI = 1.43–2.55; P < 0.0001) ( Figure 3B ).

Furthermore, IRS was significantly higher in the SD/PD group than in the CR/PR group ( Figure 4A ), indicating that IRS was negatively associated with the magnitude of response to atezolizumab plus bevacizumab in mRCC. Compared to tumors that were negative for PD-L1 IHC, tumors that were positive for PD-L1 IHC exhibited lower IRS ( Figure 4B ). We also observed enrichment of metastatic tumors ( Figure 4C ) and sarcomatoid histological subtypes ( Figure 4D ) in the IRS-high group. These findings suggest that the IRS can predict the efficacy of atezolizumab plus bevacizumab.

ROC curve was used to determine the ability of the IRS model to distinguish between immunotherapy responders and non-responders. The IRS model displayed a satisfactory performance to differentiate responders from non-responders, with an AUC of 0.822 (95% CI = 0.782–0.863) in the IMmotion151 cohort ( Figure 5A ). We also selected two external independent datasets: the JAVELIN Renal 101 cohort (patients with mRCC treated with a combination of avelumab and axitinib) and the E-MTAB-3218 cohort (patients with mRCC treated with nivolumab). When assessing survival prediction, we found that the AUC of our IRS model was 0.751 (95% CI = 0.699–0.803) in the JAVELIN Renal 101 cohort ( Figure 5B ) and 0.776 (95% CI = 0.684–0.868) in the E-MTAB-3218 cohort ( Figure 5C ). Additionally, we validated our model in the IMvigor210 cohort (patients with advanced urothelial cancer who received atezolizumab therapy) and GSE78220 (patients with metastatic melanoma who received pembrolizumab therapy), with AUC of 0.902 (95% CI = 0.868–0.936) ( Supplementary Figure 2A ) and 0.879 (95% CI = 0.7437–1) ( Supplementary Figure 2B ), respectively. Therefore, our results highlight that IRS has a favorable capability to stratify a subset of patients who will benefit from immunotherapy.

Patients with high IRS exhibited an inferior prognosis compared to those with low IRS in the JAVELIN Renal 101 cohort (HR = 1.77; 95% CI = 1.24–2.54; P = 0.002) ( Figure 5D ) and E-MTAB-3218 cohort (HR = 4.74; 95% CI = 1.31–17.2; P = 0.018) ( Figure 5E ). Similarly, patients in IRS-high group were characterized with a shorter OS than those in IRS-low group in the IMvigor210 cohort (HR = 2.59; 95% CI = 1.87–3.58; P < 0.001) ( Supplementary Figure 2C ) and GSE78220 cohort (HR = 4.22; 95% CI = 1.11–16.0; P = 0.034) ( Supplementary Figure 2D ).

A total of 26 correlated pairs between the IRS model and drug sensitivity in the GDSC database were analyzed using Spearman’s correlation analysis (38). There was significant correlation between drug sensitivity and IRS in 11 pairs, including CGP-082996, CGP-60474, and bicalutamide (Rs < –0.2, P < 0.05). In contrast, 15 pairs, including sunitinib, sorafenib, and temsirolimus (Rs > 0.2, P < 0.05), were characterized by significant correlation between drug resistance and the IRS model ( Figure 6A ). In addition, drugs whose sensitivity correlated with low IRS primarily targeted chromatin histone acetylation, p53 pathway, phosphoinositide 3-kinase (PI3K)/mammalian target of rapamycin (MTOR) signaling and protein stability, and degradation signaling pathways. However, drugs whose sensitivity was linked to high IRS mostly targeted the AKT2, IKK2, CDK2, and chromatin histone methylation signaling pathways ( Figure 6B ). Collectively, these results indicate that our IRS model is also associated with chemotherapy response in RCC.