Figure S1 illustrated the workflow for this study. We searched and downloaded ccRCC expression datasets with complete clinical annotation and mutations from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Two datasets (TCGA-KIRC and GSE36895 datasets) were analysed in this work. For TCGA-KIRC dataset, we obtained gene expression data from UCSC website (https://xenabrowser.net/datapages/) and converted them to kilobase per million values. GSE36895 dataset were downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo/). “Sva” package was performed for correcting batch effects in two datasets (10). Samples lacking complete clinical information and mutation data were excluded. Clinical information was summarised in Table S1 .

Human renal tubular epithelial cells (HK-2) and ccRCC cell lines (786-O and Caki-1) were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). These cells were cultured in DMEM or RPMI-1640 medium containing 10% fetal bovine serum and 1% streptomycin-penicillin. All cells were incubated in a sterile incubator at 5% CO2 and 37°C.

TRIzol reagent (Invitrogen, USA) was applied to isolate and extract total RNA from the cells. NanoDrop 2000 spectrophotometer (Thermo Scientific, USA) was applied for evaluating of RNA quantity control and concentration. Reverse Transcription Kit (Takara, China) was applied to reverse transcribe total cellular RNA into cDNA. ABI 7500 real-time fluorescence quantitative PCR instrument was designed for carrying out qRT-PCR process. The cycling threshold (Ct) for each gene was recorded and 2-ΔΔCt method was applied to calculate gene mRNA expression. All experiments were repeated 3 times and procedures were carried out according to reagent instructions. Primer sequences were listed in Table S2 .

Unsupervised clustering analysis were applied to identify IRF-related molecular subtypes. Consensus clustering algorithm was performed for determining the number of clusters. “ConsensuClusterPlus” package was employed to perform consistency clustering analysis (11). The process was repeated a thousand times to ensure consistency of classification.

GSVA is a non-parametric unsupervised analysis method that transforms gene expression matrices into gene set expression matrices to evaluate gene set enrichment results of transcriptome (12). Based on the “c2.cp.kegg.v6.2.symbols” gene set obtained from MSigDB database, GSVA analysis was conducted using “GSVA” package.

Single sample gene set enrichment analysis (ssGSEA) was performed to assess immune infiltration levels based on immune cell-specific gene expression. The immune gene set file is derived from Charoentong et al (13, 14). ESTIMATE algorithm calculates immune and stroma score to estimate the amount of stroma and immune cells and compute tumour purity (15). CIBERSORT is designed to calculate the composition ratio of the 22 immune cells. 22 immune cell expression data are taken from CIBERSORT website (https://cibersort.stanford.edu/) (16). Considering that CD4 naive T cells was 0 in all ccRCC samples, CIBERSORT algorithm only analysed remaining 21 immune cells.

“limma” package is applied to filter differentially expressed genes (DEGs) between clusters (17). Genes with adjusted P-value<0.001 were recognized as DEGs. “ClusterProfiler” package is intended for GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis of DEGs (18).

Univariate COX regression screened for prognosis-related DEGs. Principal component analysis (PCA) was performed for constructing IRF gene signature. PC1 and PC2 were used as feature scores to calculate IRFscore for individual samples (19). IRFscore = ∑ (PC1i + PC2i), where i represented DEGs’ expression.

The TCGA-KIRP and TCGA-KICH cohorts were used to validate the clinical performance of the IRFscore. Information on both queues can be downloaded from the online website (https://portal.gdc.cancer.gov/).

The four different immunophenotypic scores (antigen-presenting, effector, suppressor, checkpoint) are calculated separately by immunophenoscore (IPS), IPS z-score is the integration of the four, and the higher the IPS z-score, the more immunogenic the sample (20). IPS was obtained from The Cancer Immunome Atlas (https://tcia.at/home).

GDSC (https://www.cancerrxgene.org/) database contains massive genomic data on tumour therapeutics and drug sensitivity data (21). We predicted the response of ccRCC patients to five chemotherapeutic agents, including sunitinib, sorafenib, nilotinib, temsirolimus and pazopanib. “pRRophetic” package was performed for quantifying the half maximal inhibitory concentration (IC50).

Protein-protein interaction (PPI) network maps between IRFs was obtained from STRING database (22). Wilcoxon rank sum test was designed to comparative analysis of two groups, Kruskal-Wallis and one-way ANOVA was designed to calculate differences between three and more groups. Spearman correlation analysis was designed to determine correlation coefficient. Kaplan-Meier and log-rank test were performed for plotting survival curves and calculating statistical differences. Multivariate COX regression analysis was conducted to detect independent prognostic factors. “maftools” package was conducted to describe mutations. Statistical analyses were all two-sided and P<0.05 was considered statistically different. All data were analysed by R software (version 4.1.1).

IRF1-9 were included in this work. First, we analysed mRNA expression levels of IRFs in TCGA and GSE36895 cohort, respectively. IRFs were severely imbalanced in expression and the results of both databases remained largely consistent ( Figures 1A, B ). All genes were up-regulated in ccRCC except IRF6. ROC and PCA analysis indicated that IRFs can distinguish well between ccRCC and normal samples ( Figures 1C-E ). We then used two databases, CTPAC and HPA, to compare differential protein expression. CTPAC database results were consistent with the above database ( Figure S2A ). Figure S2B illustrated that in HPA database, IRF1, IRF3, IRF7-9 were upregulated in tumour, while the opposite is true for IRF6. IRF2 was highly expressed in both tissues. IRF4 and IRF5 were low or undetectable in both tissues. Furthermore, we observed that IRFs were highly correlated in expression ( Figure 1F ) and interacted with each other in PPI network ( Figure 1G ).

To validate IRFs mRNA expression, we performed qRT-PCR analysis in three cell lines. Most IRFs were more highly expressed in tumour cells ( Figure 2 ), which is generally consistent with the results above. Furthermore, we noted that IRFs were expressed with cell specificity in different cells ( Figure S2C ).

We then discussed clinical relevance of IRFs. We found that most IRFs were correlated with prognosis ( Figure 1F and S2D ). IRF6 exhibited a tumour suppressive profile and its expression was positively correlated with prognosis. In contrast, the higher the expression of other IRFs, the worse the prognosis of patients.

Using an unsupervised clustering approach, we classified ccRCC patients into different subtypes. We ultimately identified three IRF-associated molecular subtypes, termed IRF Cluster A-C ( Figure 3A and S3A-C ). Heat maps illustrated the relationship between three subtypes and clinicopathological features ( Figure S3D ). Prognostic analysis pointed to a much higher survival advantage for cluster B ( Figure 3B ). By analysing IRF expression profiles, we observed higher expression of protective factors (IRF6) in cluster B, while the opposite was true for risk factors (e.g. IRF3 and IRF7) ( Figure 3C ). This laterally explained why cluster B had a better prognosis.

GSVA analysis was performed to characterise different biological properties. Multiple immune activation-related pathways, including T and B cell receptor signalling pathways accumulated in cluster C ( Figures 3D, E and S3E ). Cluster B enriched for some matrix activation pathways, whereas cluster A was mainly associated with immunosuppression and base excision repair. We then proceeded to analyse TME immune infiltration. First, we evaluated 23 immune cell infiltrations using ssGSEA, and almost all immune cells were heavily infiltrated in cluster C ( Figure 3F ). We then ran ESTIMATE algorithm to calculate stromal and immune cell content. Apparently, cluster C had much higher immune and stromal scores, signifying that cluster C had lowest tumour purity ( Figure 3G ). However, no matching survival advantage was found for cluster C with this immune profile. Therefore, we counted the relative proportions of cell subpopulations via CIBERSORT. CD8+ T cells and M2 macrophages were more predominant ( Figures S3F-G ). Typically, the higher the expression of CD8+ T cells, the more positive the prognosis (23). Interestingly, we observed the greatest proportion of CD8+ T cells in cluster C and the lowest in cluster B, which is opposite to the prognosis. Researches have revealed that CD8+ T cells are exhausted in ccRCC and secrete numerous immune checkpoints, including PD-1 and CTLA-4. At this point, the higher the intensity of CD8+ T cell infiltration, the worse the prognosis of ccRCC (24). Here, we analysed T cell exhaustion-related immune checkpoint expression. Most checkpoints were highest in cluster C ( Figure 3H ). Combining with previous studies, we speculated IRFs may regulate T-cell exhaustion.

To further characterise biological functions of IRF-related subtypes, we filtered 547 DEGs from three subtypes and performed functional enrichment analysis ( Figure 4A ). These DEGs participated in many immune cell activation and proliferation-related pathways ( Figures 4B, C ). This implied that IRF-associated DEGs are actively engaged in immune processes and modulating immune infiltration. Subsequently, univariate COX regression analysis was performed to identify 426 prognosis-related DEGs ( Table S3 ). Similarly, we ran unsupervised cluster analysis on 426 DEGs and identified three gene clusters, termed IRF gene Cluster A-C ( Figures S4A-D ). Similarly, we compared clinicopathological characteristics and immune infiltration between different gene clusters and found that gene cluster A had superior prognostic prospects (p<0.001, Figures 4D and S4E ). Except for IRF6 and IRF8, the remaining risk genes were expressed in gene clusters in the order C, B and A ( Figure 4E ). CD8+ T cells and MDSC had lowest infiltration intensity in Cluster A ( Figure 4F ). This accounted for the greatest survival advantage of gene cluster A. Overall, the concordance of prognostic and immune infiltration characteristics among gene clusters justified this classification.

PCA analysis was conducted on 426 DEGs and IRFscore were calculated to accurately quantify individual IRF-related molecular subtypes. The samples were divided into high and low IRFscore groups following the threshold values determined by “survminer” package. Figures 5A-C exhibited the variation in attributes of individual patients in different clusters. Figure 5D demonstrated IRFs expression profiles in two groups. Prognostic analysis revealed that the higher the IRFscore, the worse the prognosis (p<0.001, Figure 5E ).

Next, we proceeded with a stratified prognostic analysis by different clinical characteristics. First, we observed a higher proportion of patients with advanced tumours were in high-IRFscore group (p<0.05, Figure S5A ). Patients with VHL, PBRM1 and BAP1 mutations also had higher IRFscore, although not statistically different ( Figure S5B ). Stratified prognostic analysis revealed that low IRFscore consistently showed marked survival advantages (p<0.05, Figure S5C ). Multivariate Cox regression analysis proved that IRFscore could be independent prognostic factor ( Table S4 ). ROC curves and nomograms demonstrated the performance of IRF scores in predicting patients’ rates at 1, 3, 5 and 7 years (AUC≥0.666, Figures 5F, G ).

To gain insight into IRFscore’s prognostic value, we further validated the effectiveness of IRFscore in predicting papillary renal cell carcinoma (KIRP) and kidney chromophobe (KICH) prognosis. Based on previous PCA results obtained from 426 DEGs, IRFscore was re-established and survival analyses were performed. In KICH, the prognosis was significantly better in low IRFscore group, while the opposite was true in KIRP (P<0.05, Figures S5D, E ). This suggested that IRFs are responsible for renal cancer progression, but for specific efficacy, it depended on cancer type.

To uncover how IRFscore works in regulating TME, we examined immune infiltration in two groups. High IRFscore group had a more significant immune infiltration ( Figures 6A, B ). Furthermore, we found that CD8+ T cells and M2 macrophages accounted for largest proportion in both groups ( Figure 6C ). Therefore, we speculated that these cells probably function primarily in ccRCC progression. Previous studies demonstrated that immune dysregulation occurs in advanced ccRCC (25), when massive exhausted T cells and M2 macrophages are simultaneously enriched in TME and substantial receptor-ligand interactions exist between two cells leading to worse prognosis (26). Table S5 listed receptors or ligands expressed by two cells. Expression analysis revealed that most co-stimulatory receptors, except for HAVCR2 and BTLA, were significantly overexpressed in high IRFscore group ( Figure 6D ). This suggested that CD8+ T cells in high IRFscore were mostly in exhausted state. However, a matching profile of M2 macrophages was not observed in high IRFscore group ( Figure S6A ). These results indicated that IRFs may not participate in interaction of exhausted T cells with M2 macrophages.

Tertiary lymphoid structures (TLS) are ectopic lymphoid tissues that surround the tumour. The higher the density of its presence, the better the patient’s prognosis (27). In ccRCC, TLS not only occurs significantly less frequently than other cancers, but also becomes dysfunctional (28). Interestingly, when TLS density and mature DCs are increased in ccRCC, a group of patients with high CD8+ T-cell infiltration and good prognosis emerges (29). This contradicted previous findings that CD8+ T cells cause worse prognosis in ccRCC (23). Therefore, scientists assumed that the emergence of TLS and mature DCs could be one reason for reduced T-cell exhaustion (30). We extracted TLS-related markers from published literatures, including three chemokines (CCL19, CCL21 and CXCL13) and two TLS-DC-related markers (HLA-DR and CD83). HLA-DR, CD83 and CCL13 were significantly upregulated in low IRFscore, while only CCL19 was downregulated ( Figure 6E ). Thus, we hypothesized that increased presence of TLS and mature DCs in low IRFscore may enhance ccRCC prognosis by reducing T-cell exhaustion.

Many studies proved that the more genetic mutations a tumour has, the more abnormal proteins it produces and the more likely the immune system is to be activated. This implied that tumour mutational burden (TMB) is somewhat predictive of immunotherapy effects (31). Furthermore, TMB can accurately predict multiple targeted and chemotherapeutic drug effects (32). Generally, the higher the TMB, the better the treatment effect. In this work, quantitative analysis and correlation analysis confirmed a positive correlation between IRFscore and TMB ( Figures S6B-C ). Survival analysis proved that lower TMB predicts a good prognosis (p<0.001, Figure 6F ). We further assessed the synergistic effect of these two scores in prognosis. Stratified survival analysis indicated that TMB and IRFscore did not interfere with each other, with IRFscore showing significant survival differences in two TMB subgroups (p<0.001, Figure 6G ). This meant that IRFscore could serve as a prognostic indicator independent of TMB.

Next, we discussed the performance of IRFscore in predicting targeted therapy efficacy. We compared estimated IC50 of five drugs ( Figures 7A-E ). Except for pazopanib, IC50 levels for remaining drugs were significantly higher in low IRFscore, meaning that low IRFscore was more sensitive to these drugs (p < 0.001). We then investigated the association between IRFscore and immune checkpoint inhibitor (ICI) therapy by IPS. Figures 7F-I depicted that four IPS scores were significantly higher in high IRFscore (p<0.001), signifying that higher IRFscore may have higher immunogenic phenotypes and be more sensitive to ICIs. Additionally, the higher the frequency of PBRM1 mutations, the better the outcome of anti-PD-1 treatment was found (26). Figure S5B demonstrated that PBRM1 mutations were more frequent in high IRFscore. Above results indicated that low IRFscore group may be more sensitive to targeted therapies, while high IRFscore subgroup were more sensitive to immunotherapy.

Numerous studies indicated that patient response to treatment options can be predicted by different tumour subtypes (33). We therefore sought to understand whether ccRCC-related phenotypes could explain why IRFs influence treatment outcome and analysed the extent to which IRFscore-related subgroups matched these tumour phenotypes. First, combining immune infiltration characteristics ( Figures 6A-C ), we hypothesized that high IRFscore group matched immune-inflamed phenotype, whereas low IRFscore group matched immune-desert phenotype. Generally, immune-inflamed phenotype was more responsive to anti-PD-L1/PD-1 therapies. In contrast, immune-desert phenotypes had no or the weakest response (34). This was consistent with our previous prediction that high IRFscore group was more sensitive to ICI therapies ( Figures 7F-I ).

Generally, targeted therapies are more effective in metastatic ccRCC (mccRCC) than other treatments (2, 3). To accurately predict the effectiveness of tyrosine kinase inhibitor (TKI) therapy in mccRCC, Benoit et al. identified four mccRCC molecular subtypes with different therapeutic effects on sunitinib based on tumour gene mutations, copy number variants (CNV) and methylation status (35). To determine whether this typing was applicable to our work, we collated the distribution of these features across two groups and summarised in Table S6 and Figures S6D-L . We considered that high IRFscore group may correspond to mccRCC 1/4 group, characterised by poor prognosis, low sunitinib sensitivity, increased methylation levels, slightly higher VHL and PBRM1 mutations, higher CNV, highly inflammatory immunosuppressive environment and low stem cell differentiation ( Figures S6D-L ). In contrast, low IRFscore group corresponded to mccRCC 2/3 group, which has the opposite characteristics. Although not all features match exactly, in general we assume that mccRCC subtypes can be applied to describe IRFscore grouping. These results pointed that IRFscore groupings can be well matched to ccRCC immunological and molecular typing, indicating that optimal treatment can be selected according to each patient’s tumour subtype.