We retrieved data on gene expression, clinicopathology, and survival of patients with KIRC and KIRP from The Cancer Genome Atlas (TGCA) database and the University of California Santa Cruz (UCSC)database. The KIRC samples with prognostic data, comprising 172 normal samples and 539 tumor samples, were integrated to form a training set. A verification set containing prognostic data on 256 KIRP patients was also established.

In total, the expression profiles of 2,483 immunity-related genes were retrieved, with the complete gene names obtained from the Immunology Database and Analysis Portal (ImmPort) database’s Gene List module (Supplementary Table S1).

The HPA contains extensive data on the transcriptome and proteome of various human specimens, including tissue, cell, and pathology atlases. Moreover, the database also provides data on protein immunohistochemistry in tumors and normal human tissue samples.

Differential gene expression analysis was used to identify the differentially expressed genes (DEGs) between the tumor and normal groups. We used the ‘limma’ software package in R (v. 3.6.3), With a false discovery rate (FDR) < 0.25 and Log2| Fold Change | > 2. The DEGs were visualized through volcano plots using the “ggplot2,” “Cairo,” and “ggrepel” packages in R (v. 3.6.3).

We used the ‘clusterProfiler’ package (v. 3.14.3) to perform GO and KEGG analyses to discover the prospective biological functions and enrichment pathways of the identified DEGs. The KEGG database incorporates information on chemical, genomic, and systematic functional data. The molecular functions, biological processes, and cell composition are gene features that are primarily described by the GO analysis. An adjusted p-value of <0.05 indicated significance. We used the “ggplot2” package (version 3.3.3) for visualization.

To identify the differentially expressed immunity-related genes (DEIGs), a total of 1734 genes associated with immunity were identified in KIRC tissues and compared with the DEGs obtained in the previous analysis. Univariate Cox regression analysis was then used to identify immunity-related genes that had significant relationships with prognosis. Following this, we used LASSO (least absolute shrinkage and selection operator) regression analysis to identify a subset of genes that could be used as independent prognostic indicators. The LASSO regression allowed us to improve our model’s accuracy and interpretability, while also addressing the issue of collinearity among the independent variables (Alhamzawi and Ali, 2018). Independent prognostic factor-related regression coefficients were then obtained from multivariate Cox regression analysis. Lastly, the formula for an IRSS, which was established using the multivariate Cox regression coefficient beta value as: I R S S = E X P g e n e 1 * β 1 + E X P g e n e 2 * β 2 + E X P g e n e 3 * β 3 + .   .   . + E X P g e n e n * β n wherein this formula, EXP denotes the expression level of the gene and β signifies the regression coefficient from the multivariate Cox regression analysis (Zeng et al., 2017).

Groups with high-IRSS and low-IRSS were established where the risk score for each TCGA-KIRC sample was computed using the above formula; the median risk score was used as the cut-off value for the establishment of these groups. Additionally, the log-rank test was used for comparing overall survival (OS) rates between the two groups, which was visualized using Kaplan-Meier (KM) curves. The AUC (area under the curve) for the receiver operating characteristic (ROC) was used to analyze the prognostic value of the six immunity-related genes identified in this study and the IRSS in KIRC patients. The horizontal and vertical axes of the ROC curve denote specificity and sensitivity, respectively (DeLong et al., 1988). The AUC provides a probability value that can be used to assess the model’s predictive accuracy and varies in value from 0.5 to 1, where the larger the value, the better the model’s accuracy. In our study, larger AUC values denote higher levels of agreement between the actual OS and the predicted OS. For these analyses, we used the “glmnet” package (version 4.1-2) and the ‘survival’ package (version 3.2-10).

Marker genes for each of the 24 types of immune cells were obtained from a research report by Bindea et al. (2013). We used the single sample Gene Set Enrichment Analysis (ssGSEA) method to identify the different immune cells that infiltrated KIRC tumors. We used Wilcoxon’s signed-rank sum test and Spearman’s correlation to analyze the associations between the gene expression profiles of the six DIEGs and immune infiltration patterns. The analysis was performed using the “Xiantao tool” module (https://www.xiantao.love/products) and the ‘GSVA’ package (version 1.34.0).

We screened for KIRC-related prognostic predictors, specifically clinical characteristics, and examined their association with the IRSS that we have developed. We paid special attention to the relationship between the OS and the IRSS, which was examined with the univariate Cox proportional risk model. Multivariate Cox regression analysis was used to establish the potential of the IRSS as an independent prognostic predictor. A nomogram was constructed by incorporating clinicopathological data from KIRC patients, including age, gender, histologic grade of the tumor, pathologic stage of the tumor, and IRSS; this was used to comprehensively evaluate patients' survival data. We used the “rms” package [version 6.2-0] and the ‘survival’ package [version 3.2-10] for these analyses.

In addition, decision curve analysis (DCA) were also generated across 1, 3, and 5 years to assess the suitability of nomograms for clinical use. The threshold probability percentage and net income are displayed in a graph (Figure 8) on the x and y-axes, respectively.

The TCGA-KIRP dataset (n = 256) was chosen to verify our IRSS model and assess the generalization of the signature since KIRC and KIRP tumors are anatomically and histologically quite similar. Every patient’s risk score in the KIRP cohort was computed and ranked using the IRSS formula that we developed using the TCGA-KIRC dataset. The TCGA-KIRP cohort’s cutoff values were used to classify the KIRP samples into groups with high- and low-risk as per the scores. The KM curve was used to compare survival across the two groups. The ROC curve was used to evaluate the precision of the signature prediction. In addition, the nomogram was used to evaluate the survival probabilities of the KIRP patients using clinicopathological data, including age, gender, histologic grade of the tumor, pathologic stage of the tumor, and the IRSS.

We first created a flow chart (Figure 1) to illustrate the steps followed in this study. The training set had samples from 539 KIRC patients along with 172 normal samples with data on gene expression, prognosis, and survival. A risk score model was developed based on the expression patterns of the immunity-related genes to predict the prognoses of the KIRC patients to examine if these genes can be used as biomarkers. The differential gene expression analysis identified 6,533 DEGs (|log2(FC)|>2 and p <0.05) as depicted in the volcano plot (Figure 2A) Supplementary Table S2. Of these, 5,055 were upregulated (logFC positive) and 1,478 were downregulated (logFC negative).

GO and KEGG enrichment analyses were conducted to investigate the prospective relationship between immune activity and gene expression in the tumor and normal tissues in the TCGA-KIRC cohort. The DEGs were enriched in a range of processes in both tissue types, with immune-related pathways making up the majority of these processes. As per the KEGG enrichment analysis (Figure 3), natural killer cell-mediated cytotoxicity, cytokine-cytokine receptor interaction, BP modules implicated in receptor-ligand activity, cytokine receptor binding, and cytokine activity, in GO enrichment analysis. Collectively, these results demonstrate that certain immune genes are differentially expressed in KIRC tumors. Supplementary Table S3, S4 provide more details on the enrichment analysis results.

According to the Venn diagram showing the overlay of the immunity-related genes and DEGs (Figure 2B); 487 DEIGs were identified (Supplementary Table S5; Figure 2B). Of these, six genes that had a strong link to prognosis were found during the LASSO regression analysis. Merging the results depicted in Figures 2D, E, the six immunity-related genes that were selected and included in the model were PAEP, PI3, SAA2, SAA1, IL20RB, and IFI30 (Figure 2C; Supplementary Table S6). The associated regression coefficients for these genes (β1–β6) were 0.067, 0.027, −0.174, 0.182, 0.038, and 0.346, respectively (Supplementary Table S7). The IRSS was developed using the following formula: I R S S = EXP   P A E P * 0.067 + EXP   P I 3   * 0.027 + EXP   S A A 2   * − 0.174 + EXP   S A A 1 * 0.182 + EXP   I L 20 R B   *   0.038 + EXP   I F I 30   * 0.346

Subsequently, an OS heatmap of the six immune-related genes was generated (Figure 4A) and the relationship between the IRSS and prognosis in KIRC patients was examined using the Kaplan-Meier (K-M) plotter. As shown in Figures 4B–G, high expression levels of PAEP, PI3, SAA2, SAA1, IL20RB, and IFI30 is correlated with poor prognosis in patients with KIRC. We also analyzed the relationships between the expression levels of the six immunity-related genes and overall survival outcomes in KIRP. The results showed that high expression levels of SAA2, SAA1, and IL20RB are linked to poor prognosis (Supplementary Figure S1). Additionally, when the KIRC samples were split into low- and high-risk groups based on the IRSS, we found that patients in the high-risk group had poorer prognoses than those in the low-risk group (Figure 5A, log-rank p < 0.001; HR = 3.12, 95% CI = 2.23–4.37). Subsequently, the single-indicator ROC curve analysis indicated that AUC (Area under curve) values were 0.764, 0.521, 0.474, 0.568, 0.957, and 0.891, respectively. (Figure 6). Following this, the effectiveness of our model in predicting OS for KIRC patients was evaluated using the ROC curves. As shown in Figure 5B, the 1-, 3-, and 5-year AUC values were 0.735, 0.708, and 0.720, respectively, which demonstrates the model’s reliability and accuracy in predicting patient prognosis.

An investigation of the relationships between the clinicopathological characteristics and the IRSS revealed that the progression of KIRC was closely associated with the six immunity-related genes identified in this study. Patients with pathological stage III–IV and histologic grade G three to four tumors, and those who died because of KIRC were found to have greater risk scores as compared to those with pathological stage I–II (p <0.001) and histologic grade G 1–2 (p <0.001) tumors and patients who survived (p <0.001; Figures 5C, E, F). No significant association between age and the IRSS were observed (p >0.05; Figure 5D). Overall, our findings demonstrated that the six immunity-related genes identified in this study are strongly linked to KIRC progression.

As shown in Figure 7, PAEP expression was significantly positively correlated with infiltration by regulatory T cells (Tregs), Th2 cells, and Th1 cells (p <0.001); PI3 expression was significantly positively correlated with infiltration by Tregs, macrophages, and Th2 cells (p <0.001); SAA2 expression was significantly positively correlated with infiltration by macrophages, Tregs, and B cells (p <0.001); SAA1 expression was significantly positively correlated with infiltration by Tregs, macrophages, and B cells (p <0.001); IL20RB expression was significantly positively correlated with infiltration of Tregs, Th1 cells, and macrophages (p <0.001); IFI30 expression was significantly positively correlated with infiltration by Tregs, NK CD56 bright cells, and T helper cells (p <0.001). Interestingly, we found that Tregs were positively associated with all the six genes identified in this study, with correlation coefficients of 0.385, 0.415, 0.399, 0.451, 0.485, and 0.333, respectively. We also analyzed the relationships between the six immunity-related genes identified here and 24 tumor-infiltrating immune cells in KIRP tissues. The results showed that the six immunity-related genes were also associated with multiple immune cell infiltration in KIRP (Supplementary Figure S2).

We used univariate and multivariate Cox regression analyses to evaluate the prognostic significance of the IRSS that we have developed. According to univariate Cox regression analysis, in the TGGA-KIRC dataset, age, histologic grade, pathologic stage, and risk score, were all prognostic predictors, although gender was not (Figure 8A). According to the multivariate Cox regression analysis (Figure 8B, log-rank p < 0.001; HR = 3.334, 95% CI = 1.785–6.229; Table 1), the risk score was an independent prognostic predictor. According to the DCA (Figure 8C), the IRSS in combination with various clinicopathological characteristics has greater clinical applicability. According to the aforementioned findings, our IRSS can be employed as a robust and innovative prognostic biomarker for KIRC.

Nomograms are useful because they simplify statistical prediction models to single numerical estimates of event possibilities suited to specific patient characteristics (Iasonos et al., 2008). Since numerous clinical characteristics have prognostic significance in clinical practice, we constructed a nomogram using numerous clinicopathological characteristics and the IRSS to predict patient prognosis in KIRC. As depicted in Figure 8D, it is possible to compute and combine the scores of each variable to generate a prognosis. According to the DCA (Figure 8E), the nomogram’s clinical application value is also higher when in combination with a variety of clinical variables.

Data from the TCGA-KIRP cohort were used as a verification dataset to validate the usability of the IRSS that we have developed for RCCs. After the individual patients’ risk scores were computed as per the IRSS equation mentioned before (Supplementary Table S8), the KIRP samples were classified into low- and high-risk groups using the median IRSS value. We found that the risk score was negatively correlated with the OS rates (Figure 9A). In addition, the KM analysis (Figure 9B) indicated that the patients in the low-risk group had significantly better prognoses than those in the high-risk group (log-rank p = 0.049; HR = 1.84, 95% CI = 1.02–3.32. The model performed well in predicting OS as demonstrated by the AUCs of the survival ROC curves which were 0.601, 0.571, and 0.518 over 1, 3, and 5 years, respectively (Figure 9C). Several clinical features including age, gender, and pathologic stage of the tumor together with IRSS were utilized to construct a nomogram for predicting prognostic survival over 1, 3, and 5 years (Figure 9D).

Collectively, the IRSS that we have developed has potential as a prognostic predictive biomarker for both KIRC and KIRP and may be useful in other RCCs also.